Inferences of viscosity from the oceanic geoid: Indication of a low viscosity zone below the 660-km discontinuity

EPSL ELSEVIER

Science Letters 15 1 (I 997)

Earth and Planetary

I35- 137

Inferences of viscosity from the oceanic geoid: Indication of a low viscosity zone below the 660-km discontinuity Motoyuki Kido a, Ond?ej hdek a Ocean Research Institute, h Deptrrtment

of Geophysics. Faculty

Uniwrsity

qf Mathematics

Received

of Tok\o.

1-15-l.

and Physics, Charles

Minami-dai. Unirersit~.

b,*

Nakano-ku.

Tokyo 164. Japan

V Hole.Yo~~iFkrich 2. 18000 Prague.

Czech Republic

13 March 1997; revised 23 June 1997: accepted 23 June 1997

Abstract We have attempted to infer details of the viscosity structure in the top 1000 km of the mantle from the igeoid and tomographic structure beneath the oceans. In order to eliminate the gravity signal from problematic masses located below the subduction zones and the continents, we have considered only the intermediate degrees of the oceanic geoid (/= 12-25). A genetic algorithm has been used to determine the family of viscosity models which give the best correlatiofi with the observed geoid. Our inversion clearly identifies the asthenosphere just below the lithosphere and also confirms tHe viscosity increase in the lower mantle predicted by previous inferences, but suggests that the main viscosity jump occurs ad a depth of about 1000 km and not at the usually stated 660-km boundary. Somewhere in the depth range of 400-1000 lkm, a low viscosity zone may exist where the viscosity decreases to a value comparable with the asthenosphere. Existenc& of such a low viscosity zone is supported by recent analysis of deep mantle anisotropy which favours a flow pattern wifh a strong horizontal component in the top part of the lower mantle. Unfortunately, the resolution of the inversion as well aslthe quality of recent seismic tomographic models are not sufficient to localize the depth and to come up with a higher accudacy for the viscosity of this low viscosity channel. 0 1997 Elsevier Science B.V. Keywords:

geoid; seismology

; tomography;

mantle: convection

1. Introduction During the last decade, analyses of the relationship between the long-wavelength geoid and the seismic heterogeneity pattern of the mantle have led to significant progress in our understanding of the viscosity structure of the mantle [l-lo]. The common conclusion of these analyses is that the ob-

’ Comesponding author. ‘Tel.: + 42 2 8576 2544; fax: +42 8576 2555: e-mail: [email protected] 0012-821X/97/$17.00 P/I SOOl2-821X(97)001

2

served non-hydrostatic geoid is consistent With a l-2 order viscosity increase in the lower mantlk which is in general agreement with the microphysic+l analysis of mantle minerals [I l] and some other g$ophysical observations [ 12-151. In spite of this bndoubted progress, details of the viscosity structuke remain uncertain. This uncertainty mainly arises from a generally non-unique interpretation of the gravity data which is additionally based on rather brecarious seismic information on the three-dimensiqnal structure of the mantle. The recent seismic topographic models are probably still contaminated b$ a significant error at higher degrees [ 16- 181 dhich may

0 1997 Elsevier Science B.V. All rights reserved. 13-l

126

M. Kido. 0. ddek/Earth

and Planemy Science Letters 151 (1997) 125-137

affect the results of geoid inversions. Another issue is the translation of seismic anomalies to densities. Although recent developments in this field are rather promising [ 19-2 I], uncertainties in determining the velocity-to-density scaling from mineral physics data are still large. Moreover, the seismic heterogeneity pattern is likely affected by chemical changes. The assumption that seismic anomalies have a thermal origin, generally adopted in the above inferences, is probably not satisfied in a large portion of the upper mantle [22], which may significantly bias the results of geoid inversions. From this point of view, seismically fast regions found beneath the continents are particularly problematic, and the magnitude of density anomalies in the subcontinental mantle as well as the nature of “continental roots” remain a controversial issue [8,22-241. Other problematic masses are located below subduction zones. Both the chemical composition and the rheology of subducted slabs are obviously different from the rest of the mantle [25]. As suggested by Ricard and co-workers [15], it is particularly the subducted lithosphere that dominates inferences of viscosity from the lowermost-degree geoid. In this paper we attempt to infer details of the viscosity structure in the top 1000 km of the mantle from the oceanic geoid at intermediate degrees (/= 12-25). We will demonstrate later that such a geoid is predominantly sensitive to the mass anomalies located in the sub-oceanic mantle and is not affected by the problematic masses associated with the subducted slabs and the continents. The use of the intermediate-degree range is, however, accompanied by increasing uncertainties in seismic tomography. To estimate the effect of these uncertainties on the resultant viscosity profiles, we compare results obtained for three different tomographic models. Seismic anomalies are recomputed to densities by using recent mineral physics data, and the robustness of the viscosity profiles with respect to the velocity-to-density scaling is tested.

2. Seismic input information The upper limit of reliable lateral resolution of contemporary global seismic tomography is usually set somewhere between degree 12 and 16. Informa-

tion at higher degrees is considered to be limited although tomographic models are often presented up to degrees 18-36 [26-301. The main drawback of models obtained with seismic body waves is deteriorated resolution at places with poor ray coverage. The resolution of these models is generally worse in a shallow mantle below oceans than below continents. More uniformly resolved images of the upper mantle are provided by tomographic models obtained with seismic surface waves. However, good lateral resolution in these models may be accompanied with a problematic depth resolution due to a spurious effect of the inversion [3 11. Laske [18] has recently demonstrated that global maps of the phase velocity only agree up to harmonic degree /‘= 6, which raises doubts as to the lateral resolution of these models at higher degrees. The question arises whether the use of intermediate-degree (/= 12-25) tomographic pattern under the oceans can lead to robust estimates of the viscosity structure. It is evident that the advantage of presumably simple seismic velocity-to-density scaling relationship in the mantle below oceans may be devaluated by insufficient resolution of seismic tomography in this part of the mantle. To assess at least partially uncertainties associated with the use of seismic anomalies in the intermediate-degree range, we will carry out our inversion for three different tomographic models. Comparison of the viscosity profiles, obtained from different tomographic models, can provide an approximate idea about possible impacts of errors in seismic models on the solution of the inverse problem. The first model, derived by Fukao and co-workers [27] and hereinafter denoted as FK, is based on a whole mantle inversion of P-wave travel times. A similar data set but a more sophisticated inversion technique was used by Zhou [29] for deriving his Pl200 model (denoted as ZH in this paper). This model only shows the top 1200 km of the mantle but with a higher depth resolution (= 50 km) than model FK. Since the resolution of models FK and ZH, constructed from tomographic inversions employing body wave travel times, may show deterioration at places with poor ray coverage. we also employ the tomographic model obtained by Zhang and Tanimoto [26] from an analysis of the surface waves (model ZT). This model images the top 500 km of the mantle. Three models of the

M. Kido. 0. Cadek / Earth and Planeta? Science Letters 151 C19971 125-137

mantle structure, used as inputs in our inversion, have been derived from the tomographic models FK. ZH and ZT discussed above. The first input model, denoted as ZT + FK, is a combination of model ZT. imposed in the top 500 km, with model FK used in the rest of the mantle. Model ZH + FK has been constructed from model ZH and FK in an analogous way, only the boundary between the models lies deeper, at a depth of 1200 km. The last input model (FK + FK) is identical with tomographic model FK throughout the mantle. The seismic inputs of our inversion can be summarized as follows: *ZT+FK(ZT:O-500km;FK:500-2900km) *ZH+FK

(ZH: 0-

1200 km; FK: 1200-2900

*FK+FK

(FK: O-2900

km)

anomalies have been considered in the top 300 km). In the asthenosphere, the velocity-to-density scaling may be strongly influenced by anelasticity [21]. To quantify the effect of this layer on the solution of our inversion. we test two extreme cases: a very low value of aln p/am L’(models B and D), anld a value comparable with the rest of the upper mantle (models A and C). The weak increase of aIn p/iUn I’ imposed in the mantle below the asthenosphere roughly corresponds to the estimate of alnl~~/U’ by Karat0 [21] recomputed to irln p/am 1‘ by using the thermal expansivity taken from [32]. ln the lower mantle we use MgO data [ 191 corrected for the effect of anelasticity [2 I]. In model E a constant conversion parameter is prescribed throughout the mantle.

km)

In all three cases, the load distribution in a deeper mantle is computed from model FK. It should be mentioned. however, that the effect of density anomalies located below a depth of 1200 km on the intermediate-degree geoid is limited.

3. Velocity-to-density

127

4. Inversion

technique

We maximized the correlation observed and predicted geoid,

c between

the

scaling

Some of the previous inferences, e.g. 161, have indicated that the preferred viscosity profile is rather insensitive to the velocity-to-density scaling. However, a suitable choice of aln p/aln L’ can improve the fit of the geoid significantly. The velocity-to-density scaling factor can be either included into the set of unknown model parameters, together with viscos ity. or considered as known from mineral physics. In this paper, models of aln p/am ~1 are assumed a priori. Five profiles of the velocity-to-density scaling factor are tested (Fig. la). The profiles significantly differ in the most problematic regions, namely the lithosphere and the asthenosphere, which enables the role of these layers in the geoid inversion to be estimated. In the lithosphere we prescribe aln p/aln I’ = 0.6 for P-waves and 0.45 for S-waves (models A and B), corresponding to the values measured in olivine [20]. As the nature of the seismic anomalies in the uppermost mantle may be rather complex [22], we also test the case when the loads in the lithosphere are completely omitted, aln p/aln I’ = 0. in models C and D (cf. [9] where no density

where &No,,, is the observed geoid and Sn(pred is the predicted dynamical geoid, both considered at degrees 12-25. In a dynamic mantle. the value of 6 Npred depends both on the distribution of density heterogeneities 6p and on the viscosity structure n of the mantle. The convective currents induced by anomalous densities deform the boundaries and internal surfaces of the mantle and these defformations also contribute to geoid anomalies. Therefore. to determine 6 Npred, we have to solve the system of equations governing viscous flow in the mantle [I]. For simplicity. we neglect the effect of compressibility, and we assume that viscosity, 7. only depends on the radius, 7 = n(r). The scalar product in the numerator of Eq. ( I ) as well as the L,-norms in the denominator are computed over the area of one of the three oceans under consideration (Adantic, Indian and Pacific). The computation of the correlation is carried out in the space domain in a 2 X 2 degree grid. The borders of the oceanic areas have been chosen far enough from continental edges and subduction zones for the geoid not to be affected by

M. Kido, 0. tadek/Earth

28

and Planeta? Science Letters 151 (1997) 125-137

..,......., S_wave

0.5 6logp/6logv

0.0

-4

-3

-2

4

0

0.5 610gp1610gv

0.0

1

2

3

log,,ll

-2

-1

0

0.5 610gp/610gv

0.0

1

2

-3

log,$l 5

1,

-

-2

-1

log,,ll *I

-

0

0.0

0.0

*1ogp,610~~

1

2

-3

-2

I

log,$l * I

*

0

1

I.

610gp,610~~

2

-3

*

I.

-2

-1

0

1

2

logdl s,,

i

0:3739,

,

-4 .3

.2

-1

‘0gdl

0

1

2

d

.2

0;3699

-1

‘ogdl

0

1

2

-3

.2

.I

‘ogd

0

1

2

-3

.2

I

‘ogdl

0

1

2

.3

-2

-1

‘ogdl

0

1

2

M. Kido, 0. &dek/Earth

and Planetap Science Letters 151 (1997) 125-137

subducted slabs and the masses beneath the continents. The harmonic expansions of the geoid and the seismic heterogeneity are bounded from below by degree 12, corresponding to the minimum characteristic length of the oceanic areas, and truncated at degree 25, which is a cut-off degree of model FK, used in the lower part of the mantle. Since both the geoid spectrum and the spectra of seismic anomalies decay with a power-law-like dependence on the degree, our inversion is mainly sensitive to lower-degree harmonics (/= 12- 18) which are probably less contaminated by errors than the high-degree terms. The viscosity model is parameterized by seven layers having boundaries at 0, 100, 200, 400, 660, 820, and 1000 km. Our intermediate-wavelength inversion is not sensitive to deep mantle viscosity which accounts for the very thick bottom layer. To find a maximum value of the correlation coefficient we use a genetic algorithm 1331 with a penalty against wildly oscillating viscosity models. The main advantage of the genetic algorithm in comparison to other inversion techniques is its ability to find the set of solutions that fit acceptably well the data, and, thus, to show a range of possible implications for geodynamics. King [6] was the first to employ a genetic algorithm for determining viscosity profiles which best fit the low-degree geoid. Our modification of the genetic algorithm differs from [6] by larger random inputs at each iterative step. For each density model we have effected 100 inversion runs started from different model populations. Together with an extensive investigation of the morphology of the fitting function in terms of two-dimensional projections, this number of runs should guarantee that the main correlation extremes in the model space have been detected. Our formulation of the inverse problem, Eq. (1). differs from the traditional minimization of the

squared misfit between data [34],

the observed

IlaN,,h, - ~NprcdllL, = min,

129

and predicted

(2)

which, in contrast with Eq. (11, constrain); not only the pattern but also the amplitudes of the; predicted geoid. However, we believe that our formulation, based on a comparison of geoid patterns, id appropriate if mineral physics information on the velocityto-density scaling is to be employed. This information is mainly qualitative, indicating the niain trends and tendencies but without giving accurate: numerical values. We have tested whether the viscosity profiles obtained by minimizing the squared misfit, Eq. (2). are significantly different from the profilles giving the maximum correlation. In spite of tie formal dissimilarity, the results of both inversionsi are found to be rather similar, showing qualitativeli the same viscosity changes. The results of inversidns mainly differ in amplitudes of the predicted geoLds. While the maximization of the correlation, Eq. I;I>, prefers - lI~N,,Y,h,llL1 if an ppropriate amplitudes II~~N~~~~IIL~ velocity-to-density scaling is imposed, th inversion based on minimization of the misfit. Eq. (921, usually tends to underestimate the amplitudes oif the predicted geoid so that II~N,,~~~II~, - ~116N~~,ll~, (see Appendix for details).

5. Results In Fig. lb we show the viscosity profiles which best satisfy constraint (1) imposed in ~the whole oceanic area. Three viscosity profiles pave been plotted for each of the tomographic mod 1s and the velocity-to-density scaling profiles. The 1ppropriate correlation coefficient between the observed intermediate-degree geoid and the predicted geo!d is speci-

Fig. 1. (a) Velocity-to-density scaling profiles employed in this study. The solid lines are for P-velocities (tomographic modelb ZH and FK) while the S-velocity anomalies (model ZT) have been scaled by using the dotted lines. The profiles have been derived from dineral physics information discussed in Section 3. (b) Viscosity profiles which give the best correlation with the observed intermediatetdegree geoid (f= 12-25) in the oceanic areas. The profiles are depicted separately for each of the tomographic models (ZT t FK. bH + FK and FK + FK) and the velocity-to-density scalings (A-E). Three different viscosity profiles that fit acceptably well the data are shown for each input density model. The appropriate correlation coefficients, E$ (l), are specified inside each panel. The horizontal thin 1)ine marks the position of the 660-km boundary. Since the geoid is insensitive to the absolute value of viscosity. only relative variations of vikcosity can be inferred from these models.

M. Kido, 0. cadek/

Earth and Planetary Science Letters IS1 (1997) 125-137

Degree 12-25

.

0.1

..

(a)

All

_________

I&_

_______...._______..___

pm.

0.2

0.3

0.4

0.5

correlation Fig. 2. Relationship between the correlation and the confidence level, computed by a Monte Carlo test for each of the oceanic regions and the functions showing the same spectral decay at degrees 12-25 as the observed geoid. A great number (100.000) of random functions with modulated power spectra have been correlated with the real geoid at degrees 12-25 and the resultant set of the correlation coefficients have been processed statistically (for details see e.g. [49] where an analogous approach is used). The confidence level corresponding to a given correlation coefficient obviously depends on the size of the area under consideration. The confidence level corresponding the whole oceanic area is depicted by a dotted line.

fied inside each block. The correlation coefficient ranges between 0.33 and 0.44. Although these values indicate that only a small part of the geoidal pattern can be explained by our inversion, the correlation is always found to be significant (Fig. 2). The viscosity profiles in Fig. lb are chosen such that they illustrate a diversity in the set of acceptable viscosity models. The differences between individual models are, however, surprisingly small. The basic characteristic of all of the viscosity models is the existence of a

Fig. 3. Comparison between the intermediate-degree geoids predicted from the density anomalies considered (a) throughout the mantle and (b) below the oceans only. The two geoid patterns are almost identical in the oceanic regions which illustrates that the prediction of intermediate-degree oceanic geoid is not affected by problematic masses associated with the subducted slabs and the continental roots. The geoids are computed for tomographic model ZT+FK, velocity-to-density scaling A and the viscosity profile giving the largest value of the correlation coefficient (see Fig. lb). Only the intermediate-degree part (/ = 12-25) of the geoid is shown. The shading marks regions of negative values. Isoline interval is 2 m.

low-viscosity zone in the vicinity of the 660-km discontinuity. The minimum viscosity is either found just below the 660~km discontinuity or somewhat deeper, at a depth of 820-1000 km. In the majority of cases the viscosity model also shows a lowviscosity zone in the upper mantle which can be identified with the traditional asthenosphere. The

Fig. 4. Viscosity profiles which give the best correlation with the observed intermediate-degree geoid (e= 12-25) in three different oceanic regions. The profiles are depicted separately for each of the tomographic models (ZT + FK, ZH + FK and FK + FK), the velocity-to-density scalings (A-E) and the oceans (Atlantic, Indian and Pacific). The appropriate correlation coefficients, Eq. (I), are specified inside each panel. The horizontal thin line marks the position of the 660-km boundary. Since the geoid is insensitive to the absolute value of viscosity, only relative variations of viscosity can be inferred from these models.

M. Kido. 0. Cadek/Earth

0.4477

u= -4

-3

-2

-1

1 0

0.4442 1

2

-3

.2

13

and Planetary Science Letters 151 (1997) 125-137

I

1 0

1

1 0.4427

1

2

0

-3

-2

-1

’

”

I 0.5277

”

I

”

I

rr

n

’

’

I

0.5320 1

0.4795

0.5230 ~ I 0.5235 ~

’

’

L

132

M. Kido, 0. Cadek/Earth

and Planetary Science Letters 151 (1997) 125-137

asthenosphere is especially pronounced in the viscosity profiles derived from seismic model ZT + FK, which is expected to provide the most reliable resolution in the top part of the mantle. The two lowviscosity zones are usually separated by a highviscosity hill peaking somewhere between 400 and 820 km. All the profiles are characterized by larger values of viscosity below a depth of 1000 km than in the upper mantle. In Fig. 3 we show the geoids computed with (a) and without (b) density anomalies below the continents and the subduction zones. The two geoid patterns are almost identical in the oceanic regions which illustrates that the prediction of intermediatedegree oceanic geoid is not affected by problematic masses associated with the subducted slabs and the continental roots. In Fig. 4 we show the viscosity profiles which best satisfy the constraint given by Eq. (1) in the individual oceanic regions. The profiles have been plotted separately for each of the tomographic models, the velocity-to-density scaling models and the oceans. Only the viscosity profile giving the best correlation is shown. It should be mentioned, however, that the other solutions with comparable correlation coefficients show basically the same features as the profiles in Fig. 4. The correlation coefficient ranges between 0.35 and 0.58. Model ZH + FK successfully predicts the geoid patterns in the Atlantic and Indian oceans (c = 0.5) but gives rather small correlation coefficients in the Pacific (c = 0.4). The best correlation in the Pacific (c = 0.48) is obtained for model FK + FK. This result is surprising if we consider that almost one third of the Pacific region is not well constrained in tomographic models ZH and FK in the depth interval of O-500 km (see plate 1 in [29]). Relatively large values of the correlation coefficient have been found for tomographic model ZT + FK in all three oceans. This suggests that the shallow part of the mantle is better resolved in tomographic model ZT than in models FK and ZH. The success of our prediction can visually be assessed from Fig. 5 where the observed (a) and predicted (b) geoids are compared. The low-viscosity zone close to the 660~km discontinuity has been found in all the three oceans. The relative decrease of viscosity in this second asthenosphere is, however, rather uncertain because

la) Observed Geoid (12-25)

I

I 1

(b) Predicted Geoid ZT+FK/A

Fig. 5. Comparison between the observed geoid (a) and the geoid predicted for tomographic model ZT + FK and velocity-to-density scaling A (b). The predicted geoid is calculated from three different viscosity profiles (Fig. 4) which give the best correlation with observations in each of the oceanic regions, Only the intermediate-degree part (! = 12-25) of the geoid is shown. The shading marks regions of negative values. Isoline interval is 2 m.

of a trade-off between the layer thickness and the viscosity contrast. Moreover, the inversion algorithm often tends to prefer extreme viscosity changes (2-3 orders of magnitude), although these changes lead to only slightly better correlation than obtained for less oscillating profiles. Thus, the viscosity decrease by 2-3 orders of magnitude, usually found below a depth of 660 or 820 km, is probably only the upper limit and the real viscosity contrast will be lower. The resolution power of the inversion is not sufficient to localize the exact depth of the low viscosity channel either. Our results indicate that the lowviscosity zone may be located either between 660 and 1000 km, or somewhat deeper, between 820 and 1000 km. A shallower location of the second asthenosphere, between 400 and 1000 km, is only obtained in the Indian ocean for model FK + FK. Such a broad low viscosity zone seems to be inconsistent with the hard transition zone, recently pre-

M. Kido. 0. &dek/Earth

133

and Planetary Science Letters 151 (1997) 125-137

dieted by mineral physics [35]. In the Atlantic Ocean. the zone is always found in a depth interval of 660- 1000 km, while the deeper location of the second asthenosphere is typical of the viscosity profiles obtained in the Pacific. The situation is ambiguous in the Indian Ocean. In this region, each of the seismic models gives a different viscosity profile. The rather different viscosity profiles which have been obtained for different oceans but the same tomographic models may be a reflection of the inherent uncertainties associated with the tomographic data. The solution of our inversion is non-unique as far as the viscosity change at a depth of 400 km is concerned: Both an increase and a decrease of viscosity at this depth can lead to a satisfactory prediction of the geoid (cf. [6]). In the case of the models with a viscosity decrease at 400 km, the viscosity hill separating the low viscosity zones is located somewhat shallower, between 200 and 400 km, or two viscosity peaks are found, one in the upper mantle and the other located below the 660-km discontinuity. King [6] has showed that the viscosity profiles inferred from the lowest-degree geoid (/= 2-8) are basically insensitive to variations in the velocity-todensity scaling. Our results confirm this conclusion also for higher degrees. Although the correlation coefficients obtained for scaling models A-E may differ significantly, the characteristic features of the resultant viscosity profiles are similar. The best correlation is usually found for scaling models A, B and E, and only rarely for models C and D, in which the density anomalies in the lithosphere have been omitted. If we only take into account model ZT + FK, which obviously provides the best resolution of the top 200 km of the mantle, then the correlation coefficients obtained for scaling models C and D are always smaller than those obtained for models A, B and E. This suggests that the omission of the density anomalies in the top part of the mantle, which is justified at low degrees 191, may not be appropriate if an intermediate-degree part of the geoid is to be predicted. A simple analysis in terms of the geoid kernels [II shows that the lowest-degree geoid ($:= 2-6) is mainly generated by mass anomalies located in a deep mantle, while the gravity effect of large scale anomalies located close to the surface is almost

perfectly compensated by a topographic response to mantle flow. Thus, the omission of the lithospheric and asthenospheric loads is not accompanied with a big change of the low-degree geoid but it can help to reduce the amplitudes of the long-waveltngth dynamic topography to values predicted for the oceanic lithosphere [36,37]. In contrast, the geoid at degrees C= 12-25 is mostly generated by loads located in the top 1000 km, and the importance of the shallow loads increases with degree. Our results indicate that even the omission of the very shallow loa&. located in the top 100 km. can lead to a deterioration of the geoid prediction. The robustness of the solution with ;respect to different parameterizations in the model and data

1 O-25 0, g

EsI,

i

Z@e

’

12-25

14-25

’

0.7140

! 0.4042

0.4849 -1

f Es. SC Y*

’ ‘I *

iy

1 zig. (0 $7 : 0.7317

0.5317

1 :

: : 0.5343 ::,:,,:1

0.5747

T-4 -3

-2 -1 0 'Wdl

1

.4 -3 -2 .l 0 bl,,~

1

-I -3 .1 -1 0

1

+

2

h,o~

Fig. 6. Testing the robustness of the solution with respect to different parameterizations of the model and diffeient ranges of the geoid spectrum. The panels in individual rows dhow solutions obtained for the number of viscosity layers specified on the left. The solutions in individual columns differ by the

ing A and the geoid in the Atlantic ocean.

M. Kido, 0. Cadek/ Earth and Planetary Science Letters 151 (1997) 125-137

ZT+FK/Atl/A

____________

8

20

.I.,~,~~.,~,~,’ ...... ._....__........._. 24 T

’ -4

-3

-2

-1 xl

Fig. 7. Geoid 4) computed normalized so due to a unit The horizontal discontinuity.

0

1

2

3

4

O-s[m/(kg/m-2)]

kernels for viscosity model ZT+ FK/Atl/A (Fig. for degrees 12, 16, 20 and 24. The kernels are that they directly represent deflections of the geoid surface density anomaly located at a given depth. dotted line marks the position of the 660-km

spaces have been tested. The results of these tests are shown in Fig. 6 where the viscosity profiles obtained for different numbers of viscosity layers and different degree ranges are shown. If only a 5-layer model is considered (top panels), the resultant viscosity profiles are rather similar to the profiles inferred from the long-wavelength geoid (see e.g. El]). In this case, the viscosity model is characterized by a stiff lithosphere and a gradual increase of viscosity with depth. If the number of layers is changed from 5 to 6 or 7, a low viscosity zone appears below the 660~km discontinuity and the value of the correlation increases (see panels in the second and the third row). The increase of the correlation is especially significant for /= 12-25 but it is relatively small if lower degrees of the geoid spectrum are also considered. Further increase of the number of layers from 7 to 8 (bottom panels) does not lead to any significant change in the viscosity profile and is accompanied by only negligible increase of the correlation. Note that the viscosity profiles obtained for a number of layers greater than 5 and for any of the three degree

ranges are rather similar, showing basically the same viscosity changes with depth. This suggests that the results presented in Fig. lb and Fig. 4 are rather independent of the parameterization. It is of interest, however, that much higher correlation coefficients can be obtained for degree range of lo-25 (left column) than for degrees 12-25 or 14-25 (middle and right columns), which indicates a close relationship between the seismic anomalies and the geoid at degrees 10 and/or 11. In our inversion, degrees 10 and 11 have not been considered. The geoid at these degrees is associated with density anomalies of wavelengths which are comparable with characteristic dimensions of the oceans. Thus, the prediction of the oceanic geoid could be affected by subcontinental masses and subducted slabs. The question remains concerning the sensitivity of our inversion to the seismic tomographic structure in a deep mantle, which in all three cases is taken from tomographic model FK. An evaluation of the geoid kernels for the viscosity profiles given in Fig. 1b and Fig. 4 suggests that the role of density anomalies located below 1000 km depth is small for degrees 12- 16 and negligible for the higher-degrees (Fig. 7). Analysis in terms of the geoid kernels also indicates that the intermediate-degree geoid is rather insensitive to the viscosity changes in a deeper mantle.

6. Concluding remarks The idea of a low-viscosity zone below the 660-km discontinuity, discussed in this paper, is not new. The existence of a mid-mantle low-viscosity zone, but located above the 660&n discontinuity, was predicted from analyses of the low-degree geoid by King and Masters [3] and Forte and co-workers [381. Zhang and Christensen [4] have demonstrated that the long-wavelength geoid can also be fitted if the low-viscosity zone is located below the 660-km boundary. Their viscosity profile shows viscosity minima at depths of 150 and 800 km and a general increase of viscosity between upper and lower mantle, thus, the features which are characteristic of the viscosity models derived in the present paper. Using a genetic algorithm, King [61 has shown that the long-wavelength geoid alone may not be able to

M. Kido. 0. &dek/Earth

and PlanetarvScience

resolve low-viscosity transition zone models from models with a decrease of viscosity below the 660-km discontinuity. The results of the joint inversion of convection and post-glacial rebound observables [lo] have suggested the existence of a broad low-viscosity zone between 200 and 1000 km with the main viscosity minimum above the 660-km discontinuity. The model with such a broad low-viscosity zone, located between 400 and 1000 km, also appears among the models inferred from the oceanic geoid (Fig. 4) but it is rather rare. Most of our viscosity profiles are characterized by two narrower low-viscosity zones separated by a high-viscosity layer. This viscosity pattern is qualitatively rather similar to the viscosity models recently inferred from the postglacial rebound data by Peltier and Jiang [39] and Peltier 1401. It should be repeated, however, that the resolution power of our inversion as well as the quality of recent tomographic models are not sufficient for the details of the viscosity structure to be reliably distinguished. Recent developments in mineral physics have indicated that the viscosity may indeed increase in the depth range of 4 1O-660 km and decrease below the 660-km discontinuity [41]. New analysis of the deep mantle anisotropy [42] prefers a convection pattern with a strong horizontal flow in the top part of the lower mantle [43]. As shown by Cserepes and Yuen [44], such a flow can particularly develop in a low viscosity zone located below an endothermic phase transition. The low viscosity zone below the 660-km discontinuity followed by a large viscosity increase at a depth of 1000 km also provides a natural explanation for the significant correlation found between the positions of subduction zones in the past 120 Myr and the seismic structure at a depth of about 1000 km [45.46]. Detailed seismic studies in the West Pacific [47,48] have shown that the subducted lithosphere can penetrate through the phase transition boundary at 660 km to the lower mantle. However, no remnants of the old lithosphere have been detected below the depth of about 1100 km. This observation is consistent with a viscosity increase around a depth of 1000 km rather than with a viscosity jump at a depth of 660 km. We are aware that the results presented here are based on seismic tomographic models which may

135

Letters 151 (19971 125-137

not accurately reflect the density structure of the Earth at intermediate degrees. Nevertheless, we believe that it is useful to show that the preseht seismic information is consistent with the existence of a low viscosity zone below the 660-km discontinuity, especially if the decrease of viscosity in this rigion is in agreement with some other independent: observations.

Acknowledgements The authors thank Y. Fukao, T. Tanimdto and H. Zhou for providing their tomographic models. This work benefitted from discussions with S. Rarato. A. van den Berg, J. Trampert and D. Yuen. Critical comments and suggestions of two anonymous reviewers are gratefully acknowledged. This research has been supported by the Czech nationaligrant No. 205/96/0212 and the Charles University grant No. 228/96.[RV]

Appendix criterion

A. Correlation

and the lea$t-squares

In Section 4, we considered two different formulations of the inverse problem. The first has based on a maximization of the correlation o’ between observed data dobb and data dpred(m), pre icted for a set of model parameters m. In the spatial 1 omain the correlation formula reads:

where R is the data domain (0 is an oceanic area and dR = sin eded+ in our case). The other inverse problem formulation was based on a minimization of the squared misfit between the observed and predicted data: S(m)

= ,,,r d”h\ - dpre.( m)]‘dR.

t.42)

Although the use of the least-squares ct$erion, Eq. (A2), can be justified from the theory of ~probability [34], it can easily be demonstrated that this formula-

136

M. Kido, 0. ?adek / Earth and Planetary

tion may between using the S can be

lead to problematic results if the correlation prediction and observation is small. By definition of correlation c, Eq. (Al), misfit expressed as follows:

S(c,r)=(l

-2cr+r’)/d$,,dL?, 0

where r is the predicted/observed

(A3) amplitude

ratio,

Denote by C the maximum value of correlation c which can be reached for a given parameterization of the model space. The minimum misfit S is then obtained for amplitude ratio r which satisfies equation ;=2(r-C)j--d&dfl=O.

W)

This implies that the inversion based on the minimization of misfit S, Eq. (A2), tends to underestimate the amplitudes of the predicted geoid so that

The discrepancies between the predicted and observed amplitudes are not important if the lowest-degree part of the geoid is considered because C is close to 1 in that case. For the intermediate degrees (/= 12-25), however, correlation C is only = 0.5 and, thus, the application of the least-squares criterion leads to a prediction with significantly underestimated geoid amplitudes. The question remains whether the solutions of the two inverse problem formulations mentioned above can be identical. Although misfit S decreases with increasing value of correlation c for a fixed r (see Eq. A3), the solutions are not generally the same which immediately follows from comparing the expressions for as/(&n and aclam.

References [ll B.H. Hager, R.W. Clayton,

Constraints on the structure of mantle convection using seismic observations, flow models, and the geoid, in: W.R. Peltier fed.), Mantle Convection,

Science Letters I51 (1997) 125-137

Plate Tectonics and Global Dynamics, Gordon and Breach, 1989, pp. 657-763. [2] Y. Ricard, Bai Wuming, Inferring the viscosity and the 3-D density structure of the mantle from geoid, topography and plate velocities, Geophys. J. Int. 105 (1991) 561-571. 131 SD. King, G. Masters, An inversion for radial viscosity structure using seismic tomography, Geophys. Res. Lett. 19 (19921 1.551-1554. [4] S. Zhang, U.R. Christensen, Some effects of lateral viscosity variations on geoid and surface velocities induced by density anomalies in the mantle, Geophys. I. Int. 114 (1993) 531541. [5] A.M. Forte, R.L. Woodward. A.M. Dziewonski, Joint inversions of seismic and geodynamic data for models of three-dimensional mantle heterogeneity, J. Geophys. Res. 99 (1994) 21857-21877. [6] SD. King, Radial models of mantle viscosity: Results from a genetic algorithm, Geophys. J. bit. 122 (1995) 725-734. [7] G. Pari, W.R. Peltier, The heat flow constraint on mantle tomography-based convection models: Towards a geodynamically self-consistent inference of mantle viscosity, J. Geophys. Res. 100 (1995) 12731-12752. [8] A.M. Forte, A.M. Dziewonski, R.J. O’Connell, Thermal and chemical heterogeneity in the mantle: A seismic and geodynamic study of continental roots, Phys. Earth Planet. Inter. 92 (1995145-5s. [91 C. Thoraval, Ph. Machetel, A. Cazenave, Locally layered convection inferred from dynamic models of the Earth’s mantle, Nature 375 (1995) 777-780. [lo] A.M. Forte, J.X. Mitrovica. New inferences of mantle viscosity from joint inversion of long-wavelength mantle convection and postglacial rebound data, Geophys. Res. Lett. 23 (1996) 1147-1150. [Ill G. Ranalli, The microphysical approach to mantle rheology, in: R. Sabadini, K. Lambeck, E. Boschi feds.). Glacial Isostasy, Sea-level and Mantle Rheology. Kluwer, Dordrecht. 1991. pp. 343-378. 1121 M. Gurnis, B.H. Hager, Controls of the structure of subducted slabs. Nature 335 (1988) 317-321. 1131 M.A. Richards, Hotspots and the case for a high viscosity lower mantle, in: R. Sabadini, K. Lambeck, E. Boschi feds.1 Glacial Isostasy. Sea-level and Mantle Rheology, Kluwer, Dordrecht, 1991, pp. 571-588. [14] Y. Ricard, R. Sabadini, G. Spada, Isostatic deformations and polar wander induced by redistribution of mass within the Earth. J. Geophys. Res. 97 (1992) 14223-14236. 1151 Y. Ricard, M.A. Richards, C. Lithgow-Bertelloni. Y. Lc Stunff, A geodynamic model of mantle density heterogeneity, I. Geophys. Res. 98 (19931 21895-21909. 1161R.J. Pulliam. P.B. Stark, Confidence regions for mantle heterogeneity, I. Geophys. Res. 99 (19941 693 l-6943. seismic [I71M.H. Ritzwoller, E.M. Lavely. Three-dimensional models of the Earth’s mantle. Rev. Geophys. 33 (19951 l-66. [I81 G. Laske, Global observations of off-great-circle propagation of long-period surface waves, Geophys. J. Int. 123 (1995) 245-259.

M. Kido. 0. tad&/

[19]

Earth and Piartetary

A. Chopelas, Sound velocities of MgO to very high compression. Earth Planet. Sci. Lett. 114 (1992)

Letters 15/ (I0971

velocities in olivine at earth mantle pressures, Science 260

[36] P. Colin. L. Fleitout, Topography of the ocean floor: Thermal heterogeneities with the lithosphere. Geophys. as.

seismic tomography. Geophys. Res. Lett. 20 (19931

1623-

Kido.

residual

Geoid anomalies and

T.

Seno.

Dynamic

depth anomalies

age-depth [38] A.M.

1626.

Lett. 17

(19901 1961-1964. [37] M.

1487-1489.

[21] S. Karato. Importance of anelasticity in the interpretation of

[22] M.P. Doin. L. Fleitout, D. McKenzie,

137

125-137

evolution of the lithosphere and interaction of Beep mantle

185-192.

[20] J.M. Zaug, E.H. Abramson. J.M. Brown. L.J. Slutsky, Sound (19931

Sciencr

topography

compared

with

in oceans and implications

for

curves. Geophys. Res. Lett. 2 1 ( 19941 7 17-720.

Forte. A.M.

Dziewonski.

R.L. Woodwardi

Aspherical

structure of the mantle, tectonic plate motionsi nonhydro-

the structure of continental and oceanic lithospheres. J. Geo-

static geoid. and topography of the core-mantle

phys. Res. 101 (19961 16119-16135.

in: J.-L. Le Moue], D.E. Smylie. T. Herring (Eds.1, Dynam-

[23] T.H. Jordan, Structure and formation of the continental tectosphere, I. Petrol. Special Lithospheric Issue (19881 1 l-37. 96 abstract). Ann. Geophys.

DC. 1993, pp. 135

166.

[39] W.R. Peltier. X. Jiang. Glacial isostatic adjustmejn and Earth

(19961

rotation: Refined constraints on the viscosity of the deepest

in subducted

[40] W.R. Peltier. Mantle viscosity and ice-age ice slieet topogra-

14 (Supplement)

I I. [25] M.

ics of the Earth’s Deep Interior and Earth Rotatidn, Geophys. Monograph. 72. AGU, Washington

mantle

[34] G. Pari. W.R. Peltier. Dynamics of the s&continental (EGS

boundary,

mantle, J. Geophys. Res. 101 t 19961 3169-3290. Riedel.

S. Karato.

Grain-size

evolution

lithosphere associated with the olivine-spine1

transformation

and its effect on rheology. Earth Planet. Sci. Lett. 148 (1997) 27-33.

1359-136-1.

mantle minerals. in: D. Croaaley. A.M. Soward (Eds.1. Earth’s

[26] Y.-S. Zhang. T. Tanimoto. Global Love wave phase velocity variation and its significance to plate tectonics, Phys. Earth Planet. Inter. 66 (19911

160-202.

1271 Y. Fukao. M. Ohayashi.

H. moue, M. Nenbai,

Subducted

97 (199214809-4822. and Rayleigh

waves between 40 and 150 seconds,

Geophys. J. Int. 122 (I 995) 675-690. km of the mantle. J. Geophys. Res. 101 (1996127791-27810. Dziewonski.

locity structure: compatibility

Imaging upper mantle S veof different

Trans Am. Geophys. Union 77I46) (1996) H.-C.

B.L.N.

Kennett.

How to reconcile body-

Int. 125 (19961 279-248. [43] S. Karato. Seismic anisotropy in the deep mantle, boundary

Nataf,

Pageoph. 1997. [44] L. Cserepes. D.A.

Yuen,

Dynamical

consequences of mid-

mantle viscosity stratification on mantle flows iith

1291 H. Zhou, A high-resolution P-wave model for the top 1200

(311 Y. Ricard,

[42] J.-P. Montdgner.

layers and the geometry of mantle convection. submitted to

1281 J. Trampert, J.H. Woodhouse. Global phase velocity maps of

[30] G. Ekstom. A.M.

Deep Interior, Gordon and Breach. 1996. pp. 22t)-272. wave and normal-mode reference earth models. ‘Geophys. J.

slabs stagnant in the mantle transition zone. J. Geophys. Res

Love

phy. Science 273 (19961

[41] S. Karato. Phase transformation and rheoiogical properties of

J.-P. Montagner.

data sets. EOS, F483.

[45] L. Wen. D.L. Anderson. The fate of the slabs inferred from seismic tomogrdphy and 130 Ma subduction. &drth Planet. [46] H. Kyvalov6.

0.

Cadek,

D.A.

Yuen.

Correladion analysis

between subduction in the last 180 Myr and la/era1 seismic

tion with seismic data, J. Geophys. Res. 101 (19961 8457--

structure in the lower mantle: Geodynamical

8472.

Geophys. Res. Lett. 22 (19951 1281-128-t.

[32] A. Chopelas, R. Boehler. Thermal expansion measurements at very high pressure. systematics and case for a chemically 16 (19891

1347..

1350. Reading. MA,

1989.

[34] A. Tarantola. Inverse Problem Theory, Elsevier. Amsterdam. 1987. garnets: systematics and implications for the rheology of the mantle transition zone, Earth Planet Sci. Lett.

130 (19951

implications,

F. Niu, Seismic evidence for ai920-km

dis-

continuity in the mantle. Nature 371 (19941 3011-305. [48] T. Sakurai. Whole-mantle

P-wave tomography and differen-

time measurement, Master Thesis. Wniv. Tokyo.

1996. [49] T.W.

Ray.

D.L.

Anderson.

Spherical disharmpnics

Earth sciences and the spatial solution: Ridges. slabs. geochemistry

[35] S. Karate. Z. Wang, B. Liu, K. Fujino, Plastic deformation of

13-30.

[47] H. Kawakatsu.

tial PP-P

[33] D.E. Goldberg. Genetic Algorithms in Search, Optimization, and Machine Learning. Addison-Wesley,

( 1997)

181-184.

Sci. Lett. 133 (19951 185-198.

The tree-dimen-

sional seismological model a priori constrained: Confronta-

homogeneous mantle, Geophys. Res. Len.

the en-

dothermic phase transition, Geophys. Res. Lett. 24

and tomography

phys. Res. 99 (199419605-9614.

correlations.

in the

hotspots. J. Geo-