slab

How Strong are the Subducted Slabs ?

Shun-ichiro Karato1, Michael R. Riedel2 and David A. Yuen1,3

1: Department of Geology and Geophysics, University of Minnesota, Minneapolis, MN 55455, U.S.A.
2: GeoForschungsZentrum, Telegrafenberg, Potsdam, D-14473, Germany.
3: Department of Geology and Geophysics, Minnesota Supercomputer Institute, University of Minnesota, Minneapolis, MN 55415-1227 , U.S.A.

submitted to Science, June 1998

Abstract

Style of deformation of subducted oceanic lithosphere shows a good correlation with thermal parameter [(age of lithosphere) x (rate of subduction)]: slabs with large thermal parameters tend to deform more above the 660 km discontinuity. Causes for this correlation are investigated by taking into account the stress and grain-size sensitive rheology of slab materials. High stresses and large grain-size reduction, in addition to the latent heat release associated with the olivine-spinel transformation, can cause significant weakening of slabs. However, both the stress and grain-size effects are significant only for slabs with large thermal parameters, providing a possible explanation for the observed correlation between thermal parameters and style of deformation.

Introduction

The oceanic lithosphere subducts into deep mantle due to the negative buoyancy forces through fractures in the near surface layer and plastic flow in the central portions (1). Once subducted, lithospheric plates will encounter another resistant force near the 660 km discontinuity. Evidence for these opposing forces in the deep transition zone has been well documented through the analysis of focal mechanisms of deep and intermediate earthquakes (2), but the ultimate fate of subducted lithosphere has become detectable only recently with the development of seismic tomography (3). These studies have shown that, although most of the subducted materials eventually reach the bottom of the mantle (the D" layer), a wide range of style of slab deformation seems to occur near the 660 km discontinuity (3). In some areas (such as beneath the Americas), slabs penetrate into the lower mantle without much deformation around the 660 km discontinuity, whereas slabs are deformed and lie above the 660 km discontinuity in some western Pacific subduction zones (such as the Japan, the Izu-Bonin and the Tonga). In the Izu-Bonin subduction zone, for example, horizontally lying slabs on the 660 km discontinuity extend as much as ~ 2000 km toward the west of the subduction zone (3). Van der Hilst and Seno (4) suggested that the main cause for such a variety of slab behaviour is due to the differences in the trench migration velocity. However, a closer inspection of the relation between slab deformation and subduction parameters reveals that such a correlation is rather weak (see below in the next section). Furthermore, a simple physical consideration indicates that even if trench migration may affect slab geometry, it would be possible only for a weak slab (5). No significant deformation of slabs would occur, if slabs are significantly stronger than the surrounding mantle, as one might expect from the temperature-dependent rheology and the estimated temperature distribution in slabs. Contrary to this expectation, significant slab deformation appears to occur selectively in slabs where temperatures are likely to be low (old and fast subducting slabs in the western Pacific). The purpose of this paper is to provide a physical model to understand this paradox based on the mineral physics of the creep behaviour of Earth materials. We first review geodynamical observations of slab deformation and examine correlations between slab deformation and subduction parameters. The role of slab rheology will be highlighted along with trench migration in controlling slab deformation. Second, we examine slab rheology based on the available experimental data and theoretical models of deformation for minerals with the emphasis on the role of phase transformations. Third, we combine this mineral physics model with slab thermal models to estimate the deformability of slabs. Finally, the model is applied to explain the observed diversity of slab deformation in the transition zone.

Deformation and Rheology of Subducted Slabs

Geophysical observations

Deformation of slabs can be investigated through (i) the geometry of Wadati-Benioff zone that represents the zone of high seismic activities (6), (ii) the results of seismic tomography showing the regions of higher than average seismic wave velocities (3) and (iii) the analysis of the focal mechanisms of deep and/or intermediate depth earthquakes that provides a constraint on the geometry of deformation (7). Both (i) and (iii) are, however, applicable only to slabs shallower than ~660 km. Seismic tomography provides a valuable tool for inferring slab deformation. A serious problem in using the results of seismic tomography is its limited resolution. In the regional tomography of western Pacific, taken from Fukao and his co-workers, discrete blocks with a size ~110x110 (lateral size) x 200 (depth) km3 were used in the analysis. In contrast, much coarser blocks (400x400x200 km3) are used in van der Hilst et al's global tomography. Therefore, it is not straightforward to compare the results from two groups. Although some of these observations have large uncertainties, these results provide us with a means of classifying the style of slab deformation with emphasis on the deformation near the 660 km discontinuity. Using the results of high resolution seismic tomography (3) and the geometry of the Wadati-Benioff zones as defined by the seismicity (6), we classify the slabs into two categories (Table 1). [I] (type A) Slabs that penetrate directly into the lower mantle without much deflection in the transition zone. Slabs in the Americas (3), in the north Kurile (3), Java (3) and in the Kermadec (8) belong to this category. The slab in the Mariana (and perhaps in the Kermadec) subduction zone appears to penetrate deep into the lower mantle, but the slab image becomes broad in the lower mantle (3) suggesting significant deformation in the lower mantle. [II] (type B) Slabs that are significantly deflected above the 660 km discontinuity. The slab in the Izu-Bonin subduction zone is the best example where a horizontally lying slab above the 660 km discontinuity extends as much as ~2000 km toward the west. The slab in the Tonga is another example, although the horizontal extent of the slab deflection is not as large as that in Izu-Bonin. The New Guinea slab, the south Kurile and the Japan slab are also an example of this type (3). The degree of deformation of slabs in the Americas is difficult to estimate because of poor resolution. A regional tomography by Engdahl et al. (9) appears to suggest some deformation but the details are poorly constrained. It is, however, obvious that the degree of slab deformation near the 660 km discontinuity in these subduction zones is significantly less than those in the western Pacific. A wider variety of deformation behaviour has been observed based on laboratory experiments (10), numerical modeling (11) and on the Wadati-Benioff zone morphology (6). However, we adopt this simple classification because the main issue in this paper is the deformation of slabs at around the 660 km discontinuity and because the limited resolution of current seismic tomography does not seem to justify further sophistication. Such a variety of slab behaviour may be related to the various parameters of subduction. First, the age and velocity of subduction will control deformation behaviour through their effects on temperature which affects both slab strength and the buoyancy forces. Wortel and Vlaar (12) suggested a slab thermal parameter defined by

(t: age of subduction, vsub: velocity of subduction) to characterize the minimum temperature of slabs. Second, trench migration velocity (4), vm may influence the deformation behaviour of slab as suggested by laboratory studies (10) and by numerical modelling (11). We estimated the slab thermal parameters and trench migration velocities, using a compilation of current subduction parameters (6) and subduction history (13). Fig. 1 summarizes the correlation among slab deformation and subduction parameters. Note that in estimating the trench migration velocity and the slab thermal parameter, we choose the data corresponding to the slab properties from the surface to around 660 km depth. This is done because the geometry of slabs near the 660 km boundary in the transition zone is considered to be determined by the properties of slabs in that portion and the kinematic boundary conditions such as the trench migration velocity during the period when a slab interacts with the 660 km discontinuity. For example, for a slab in the Izu-Bonin subduction zone this corresponds to a period of ~10 my ago to the present. Several points can be noted from Fig. 1. First, despite an earlier claim (4), the local variation in deformation behaviour of slabs does not seem to show a clear relation with the trench migration velocity. A sharp change in slab geometry observed in the Izu-Bonin - Mariana, northern to southern Kurile and in the Tonga-Kermadec trench does not seem to correlate with the variation in the trench migration velocity. In the Izu-Bonin - Mariana subduction zone where the most dramatic change in the style of slab deformation occurs, the trench migration velocities are similar in each region from ~10 my ago to the present, although there was a significant difference during ~30 - 17 my year ago due to the opening of the Shikoku basin (4). Similarly, the change in the style of subduction from the northern to southern Kurile (3) does not have an obvious relation to the trench migration. Note also that in the Tonga-Kermadec trench, very significant trench migration is observed today back to ~50 my ago (13). Yet, no significant slab deflection is observed in the Kermadec. Furthermore, in the Americas, where trench migration similar to some of the western Pacific subduction zones occurs, evidence for slab deformation is less than that in the western Pacific. Thus the correlation between trench migration and slab deflection appears rather indirect and weak. Second, there is a broad correlation between the slab thermal parameters and slab deformation: all the slabs with significant deformation around the 660 km discontinuity are in the western Pacific and have large slab thermal parameters, i.e., low temperatures and/or high subduction velocities. The rather dramatic local variations observed in the slab geometry suggest that the nature of deformation of a slab is controlled to a large extent by the local force balance and/or local tectonic history (kinematic boundary conditions). Thus, factors controlling slab deformation may be sought in some local parameters rather than global circulation (14). Based on the rather weak correlation between the observed trench migration and slab deformation, we consider that the role of trench migration in slab deformation is rather secondary. In contrast, the difference in the style of deformation between slabs with large and small thermal parameters seems more fundamental. Other sets of observations on slab deformation or rheology include (i) the geoid anomalies around the subduction zones of western Pacific (15), (ii) the magnitude of stress drops of deep earthquakes (16, 17) and (iii) a sudden steepening of slab dip angle at the seismicity cutoff (18). The analysis of geoid anomalies in the western Pacific (15) suggests that these slabs are weak and not much stresses are transmitted along the slabs. In addition, stress drops associated with deep earthquakes in the Tonga subduction zone (16) are ~80 MPa, which are higher than those associated with the shallow ones (19), but not high as high as would be, if the slab stress is determined by viscosity that changes only with the temperature and pressure. Note also, however, the stress drop associated with the Bolivia deep earthquake of 1994 is significantly higher (~110 MPa) than those in the Tonga (17). All of these observations, taken together, strongly suggest that deep slabs, particularly those in the western Pacific, are rather weak. Similarly, the results of numerical modeling of slab deformation are consistent with the observations of slab morphology and focal mechanisms of deep earthquakes only when rather weak temperature-dependent rheology is assumed (20). These observations indicate that (i) a slab with a large thermal parameter (low temperature and/or high subduction rate) is easier to deform and that (ii) a high velocity of trench migration favors slab deformation (or lying over the 660 km discontinuity) but its effect is rather indirect. Therefore we suggest that variation in the slab rheology through the variation in the slab thermal parameters may play an important role in controlling the variety of style of slab deformation. However, the above observations imply a seemingly strange temperature dependence of slab strength: a cold slab appears to be weaker than a warm slab. To understand this paradox, we need to investigate the relation between slab rheology and its deformation and the dependence of slab rheology on temperature under conditions characteristic of the deep transition zone.

Slab deformation and slab rheology

The ability of subducted slabs to penetrate into the lower mantle or the nature of deformation of slabs near the 660 km discontinuity depends on various factors including the viscosity of slabs, the magnitude of the force acting on them and the kinematic boundary conditions such as the velocity of trench migration (4, 5, 11). Let us consider the bending of a viscous thin sheet embedded in a less viscous fluid which is inclined to the horizontal plane with a dip angle,

. We consider deformation of such a slab upon the action of a force, F, at its end (Fig. 2). For simplicity, we neglect viscous forces exerted by the neighbouring less viscous mantle. Then a simple analysis shows that the characteristic time for the deflection, t, is given by (21),

(1)

where

is the strain associated with bending and D is the flexural rigidity defined as

(2)

where h is the thickness of the slab, y is the coordinate normal to the slab surface measured from the center of a slab (Fig. 2) and h is viscosity which generally depends on strain-rate and M is the bending moment acting at the edge of the slab. When this time scale is comparable to or less than the time scale of subduction (~10 my), significant slab deformation becomes possible. The bending moment acting on a slab depends on the force and the geometry of a slab as,

(3)

where F is the vertical force and

is the dip angle. According to equation (1) through (3), the deformability of a slab depends on (i) rheology (through D), (ii) force and (iii) geometry such as the dip angle near the 660 km discontinuity. Resistance forces acting at the tip of a subducting slabs may result from the excess buoyancy forces associated with the deflection of the 660 km discontinuity and a resistant force due to the high viscosity of the deep transition zone or the lower mantle (22). These resistant forces increase with the velocity of subduction. For a reasonable range of parameters (L~700-900 km, F~1013 Nm), we obtain M~1023-1024 (Nm) (23). For these values of the bending moment, a critical flexural rigidity below which significant deformation occurs is Dc~1038 Nms. The flexural rigidity, D, is sensitive to the rheological properties of slab materials, as well as the temperature and stress distribution. We calculate effective viscosity using the available experimental data on the rheology of olivine and spinel (Table 2), a theoretical model of grain-size evolution associated with the olivine-spinel phase transformation (25) and on the stress-strain distribution in a slab obtained by solving the moment balance equation for a given force, F (26). We consider three mechanisms of deformation: the Peierls mechanism, the power- law creep and diffusion creep. Flow law parameters used are given in Table 2 (27). Note that these mechanisms have markedly different dependence of creep strength on temperature (and pressure), stress and grain-size. The Peierls mechanism and power-law creep are nonlinear in the stress-dependence. Temperatures in slabs are calculated using the analytical model by McKenzie (28) with additional corrections for adiabatic heating and latent heat release due to the olivine-spinel transformation. For the initial temperature profile, we assume a thermal model corresponding to 100 my old oceanic lithosphere (temperature profiles are similar to the model used by Riedel and Karato (25)). Initial grain-size of olivine is assumed to be 3 mm. Grain-size, however, will change during a phase transformation. We use a theoretical model developed by Riedel and Karato (25) which predicts that the grain-size of newly formed spinel is very small when transformation occurs at low temperatures. For example, in the center of a cold slab, newly formed spinel is predicted to have ~1 micron grain-size. This would lead to a reduction of the viscosity by as much as ~1010 times in comparison to the viscosity for ~3 mm grain-size. Therefore a weak core will occur in the central portions of a slab that is made of fine-grained spinels. The extent to which this zone develops depends on the thermal parameter (Fig. 3). Dominant deformation mechanisms in subducting slabs vary across a slab cross section as shown in Fig. 3. This pattern depends also on the velocity of subduction. Effective viscosity varies laterally and with the velocity of subduction (Fig. 4). Flexural rigidity, D, generally decreases with depth as a result of increase in temperature (Fig. 5). Flexural rigidity decreases abruptly at ~400 km deep due to the effects of latent heat release. Effect of grain-size reduction is large for cold fast slabs. This effect is important particularly in the regions deeper than ~500 km. Variation of the flexural rigidity with the subduction velocity is partly due to the assumed increase in bending moment with velocity (23). The bending moment can be modified by changing the force as well as the dip angle. A high dip angle near the 660 km discontinuity leads to a small bending moment.

Discussion

Effects of trench migration?

So far, the variation of trench migration velocity has been the only explanation for the diversity of slab deformation behaviour (4, 10, 11, 29). However, this kinematic explanation has two major problems. First, such a kinematic model would work only when slab deformation is dynamically possible. In fact, all the experimental (10) or theoretical modeling of slab deformation (11, 20, 29) assumed low slab viscosities for inducing its easy deformation. However, if deformability of slabs can vary from one slab to the other, then the variation in deformability must also be examined in addition to the influence of kinematic boundary conditions. Second, the observational basis to infer the importance of trench migration is not very strong (Fig. 1). There is no obvious reason to suggest a causal relation between trench migration and slab deflection from the seismological observation and plate kinematics as discussed before. The contrast between western Pacific slabs (slabs with large thermal parameters) and those in the Americas (slabs with small thermal parameters) appears more fundamental. We suggest that the rapid trench migration observed in some of the western Pacific subduction zones is a result of weak slabs and that the rapid trench migration is not a cause of deformation.

Deformation by earthquake faulting?

The rheological models considered here assume homogeneous deformation. However, there is evidence for localized deformation in the deep slabs. Several analyses suggest that a significant fraction of strain in subducted slabs is taken by seismic faulting (16). In this case, the assumption of homogeneous rheology does not work in its simplest form. However, deep earthquakes can explain only a part of total strain. This is true particularly for the deep portions of slabs beneath seismicity cutoff where the presence of soft slabs is inferred (18). Therefore, weak rheology in the ductile regime must play a significant role in many slabs. Weakening in the ductile regime, where more or less homogeneous deformation occurs, may also be important for the generation of deep earthquakes. Many mechanisms of deep earthquake faulting assume some form of rheological weakening (30). In addition, the complicated rheological structures of deep slabs predicted by the present model suggest that deep earthquakes may originate in the deep slabs, where high stress concentration occurs.

Weak zones?

The presence of zones of weakness such as fossil fault planes might influence the deformability of slabs (31). Such fossil fault zones may survive deep into the transition zone particularly in the cold portions of slabs. Fossil fault zones may be weaker than other regions due to small grain-size and/or possibly to higher content of hydrous minerals (or higher water content). Deformation along these zones may also contribute to slab deformation. However, it is not straightforward to explain observed variation in the style of slab deformation by this mechanism, which would predict positive correlation between seismicity at trenches and slab deformation in the transition zone. Such a correlation is not consistent with the observation (2).

Explanation of the variety of slab deformation behaviour

How could we explain the diversity of slab deformation behaviour by our model? First, our model provides a natural explanation for a puzzling observation that slabs with a large thermal parameter (slabs with low temperatures and/or large bending moment) tend to be more prone to deformation. In our model, this occurs because of the strong temperature dependence of grain-size after the olivine-spinel transformation in addition to the higher stresses. Obviously, if slab temperature is very warm, viscosity is so low and negative buoyancy forces are so small that it will not penetrate through the discontinuity. Likewise, if slab temperature is too cold, then no phase transformation will occur and hence slabs will remain strong and then no slab deformation will be possible. Within a reasonable range of slab temperatures, slabs with moderately high temperatures appear to have high enough strength to resist deformation near the 660 km discontinuity, whereas slabs with colder temperatures tend to become weak enough to cause significant deformation. Thus, there will be a temperature range (a "window") in which slabs will be strong enough to penetrate through the 660 km discontinuity without much deformation. The 660 km discontinuity will act as a filter for slabs for rheological reasons. This provides a possible explanation for the observed correlation in shown in Fig. 1. Among the slabs that are weak enough, the style of deformation appears also to depend on some other factors including trench migration velocity and other geometrical factors such as the dip angles.

Concluding remarks

In most of previous numerical modeling or laboratory experiments of convection or slab deformation, simple rheology of Earth materials such as constant viscosity or a very weakly temperature dependent viscosity was assumed. This is firstly due to computational difficulties in incorporating more "realistic" rheology but also because if one uses such "realistic" rheology (strongly temperature dependent rheology), one would not see any deformation of slabs (5). These studies have led to a notion that it is mainly the kinematic boundary conditions that control the variation in deformation behaviour of slabs and the importance of the rheology of slabs has not been well appreciated. We have shown from both theoretical and observational backgrounds that realistic rheology of slabs must be considered as an important parameter that plays a critical role in slab deformation. One of the key notions in our model is a strong dependence of the grain-size of spinel on temperature and the highly stress-dependent rheology such as the Peierls mechanism. We used a theoretical model to estimate grain-size (25), and rheological data on analogue materials of silicate spinel. The direct high-pressure deformation experiments of silicate spinel have recently been performed, which support our model (32). However, some of the rheological parameters remain poorly constrained and further experimental studies under high pressures are needed to improve our understanding of deep mantle dynamics. We have focused our attention to slab deformation above the 660 km discontinuity. This is firstly because the deformability of slabs in that portion controls their fate upon collision with the 660 km discontinuity, but secondly because much less is known about the rheology and the kinetics of phase transformation for the lower mantle minerals. However, most of the results of seismic tomography show significant thickening of slabs in the lower mantle (3), suggesting that slabs there are weak. Effects of grain-size reduction, similar to the one investigated in this study, may also operate in the lower mantle (33). Finally, if some slabs are indeed very weak in the deep transition zone as our model suggests, then we must conclude that the dynamics of plate motion must be largely decoupled with the dynamics of deep processes. The observed correlation between plate velocities and subduction zone parameters (34) should then be attributed to the control of plate velocities by forces acting on relatively shallow portions of the subduction slabs.

References and Notes

(1) C. Goetze and B. Evans, Geophys. J. R. astr. Soc., 59, 463 (1979).

(2) B. L. Isacks and P. Molner, Rev. Geophys. Space Phys. 9, 103 (1971); C. Frolich, Ann. Rev. Earth Planet. Sci. 17, 227 (1989).

(3) K. Okino, M. Ando, S. Kaneshima and K, Hirahara, Geophys. Res. Lett. 16, 1059 (1989); R. D. van der Hilst, E. R. Engdahl, W. Spakman, G. Nolet, Nature 353, 47 (1991); Y. Fukao, M. Obayashi, H. Inoue, M. Nenbai, J. Geophys. Res. 97, 4809 (1992); S. Grand, J. Geophys. Res. 99, 11591 (1994); X-Y. Ding and S. P. Grand, J. Geophys. Res. 99, 23767 (1994); R. D. van der Hilst, S. Widiyantoro and E. R. Engdahl, Nature 386, 578 (1997); Sakurai, T., "Whole mantle P wave tomography and differential PP-P time measurement", MSc Thesis, University of Tokyo, pp.31 (1996).

(4) R. D. van der Hilst and T. Seno, Earth Planet. Sci. Lett. 120, 395 (1993).

(5) U. R. Christensen and D. A. Yuen, J. Geophys. Res. 89, 4389 (1984), G. F. Davies, Earth Planet. Sci. Lett. 133, 507 (1995), S. King and J. Ita, J. Geophys. Res. 100, 20211 (1995).

(6) R. D. Jarrard, Rev. Geophys. 24, 217 (1986).

(7) D. Giardini and H. J. Woodhouse, Nature 307, 505 (1984).

(8) R. D. van der Hilst, Nature 374, 154 (1995).

(9) E. R. Engdahl, R. D. van der Hilst and J. Berrocal, Geophys. Res. Lett. 22, 2317 (1995).

(10) R. W. Griffths, R. I. Hackney and R. D. van der Hilst, Earth Planet Sci. Lett. 133, 1 (1996); L. Guillou-Frottier, J. Buttles and P. Olson, Earth Planet Sci. Lett. 133, 19 (1996).

(11) U. R. Christensen, Earth Planet. Sci. Lett. 140, 27 (1996).

(12) M. J. R. Wortel and N. J. Vlaar, PAGEOPH 128, 554 (1988).

(13) C. Lithgow-Bertelloni and M. A. Richards, Rev. Geophys. 36, 27 (1998).

(14) B. H. Hager and R. J. O'Connell [Tectonophysics 50, 111 (1978)] argued that subduction zone dip angles are controlled mainly by the global flow pattern. Our focus in this paper is the cause of more local variation in the nature of slab deformation and hence the local variation in subduction parameter could play more important roles. Furthermore, neglecting the high viscosity of slabs in Hager and O'Connell's model leads to an over estimate of the role of global circulation.

(15) L. Moresi and M. Gurnis, Earth Planet. Sci. Lett. 138, 15 (1996).

(16) N. Sugi, M. Kikuchi and Y. Fukao, Phys. Earth Planet Inter. 55, 106 (1989).

(17) M. Kikuchi and H. Kanamori, Geophys. Res. Lett. 21, 2341 (1994). The stress drop for the Bolivia earthquake was estimated using the same technique as the one used in (16). The large stress drop for the Bolivia earthquake is the result of a large magnitude (Mw = 8.3) and a small fault area (~40 x 40 km2).

(18) J. C. Castle and K. C. Creager, Earth, Planet, Space, in press.

(19) W. Chung and H. Kanamori , Phys. Earth Planet. Inter. 23, 134 (1980).

(20) P. Machetel and P. Weber, Nature 350, 55 (1991); W. R. Peltier and L. P. Solheim, Geophys. Res. Lett. 19, 321 (1992); S. Honda, S. Balachander, D. A. Yuen and D. M. Reuteler, Science 259, 1308 (1993); P. J. Tackley, D. J. Stevenson, G. A. Glatzmaier and G. Schubert, Nature 361, 699 (1993); W. C. Tao and R. J. O'Connell, Nature 361, 626 (1993); S. Zhong and M. Gurnis, Science 267, 838 (1995); G. F. Davies, Earth Planet. Sci. Lett. 133, 507 (1995); G. Houseman and D. Gubbins, Geophys. J. Int. (1997); U. R. Christensen, J. Geophys. Res. 102, 22435 (1997). In all of them, a much weaker temperature dependence of slab rheology than typical laboratory data on creep [e.g., S. Karato, Phys, Earth Planet. Inter. 55, 234 (1989)] was assumed to reproduce observed significant deformation of slabs in the transition zone.

(21) Biot's theory of deformation of a thin viscous plate [e.g., D. L. Turcotte and G. Schubert, Geodynamics, John Wiley & Sons, 1982] is used to estimate the time scale of deformation.

(22) One possible mechanism is the presence of a significant buoyancy force near the 660 km discontinuity due to the negative Clapeyron slope for the spinel to perovskite plus magnesiowustite phase transformation [E. Ito and E. Takahashi, J. Geophys. Res. 94, 10637 (1989)] or due to the presence of a buoyant garnet in the lower mantle [A. E. Ringwood and T. Irifune, Nature 331 , 131 (1988)]. Another possible mechanism for a high resistance for subduction is an increase in viscosity in the deep mantle. A significant increase in viscosity in the deep mantle was suggested by the analyses of geoid anomalies [B. H. Hager, J. Geophys. Res. 89, 6003 (1984); B. H. Hager and M. A. Richards, Phil. Trans. R. Soc. Lond. A328, 309 (1989)] and also by the analysis of post-glacial rebound [M. Nakada and K. Lambeck, Geophys. J. Int. 96, 495 (1989)]. However, the interpretation of these observations is not unique and much smaller increase in viscosity has also been proposed [e.g., J. X. Mitrovica and A. M. Forte, J. Geophys. Res. 101, 2751 (1997)]. Mineral physics study of plastic deformation suggests that the bottom portion of the transition zone may have a high viscosity [S. Karato, Z. Wang, B. Liu and K. Fujino, Earth Planet. Sci. Inter., 130, 13 (1995)].

(23) Forces acting on slabs have been estimated from buoyancy forces or viscous resistant forces to be ~1013 N/m [G. F. Davies, J. Geophys. Res. 85, 6304 (1980)]. Moment depends also on the dip angle and the wave length (curvature) of deformation. Assuming the length of slabs of ~700-900 km, one gets M~1023-1024 Nm. Because force F is likely to increase with the velocity of subduction, v, we consider M=1.6x1023-4x1023 Nm corresponding to v = 4-10 cm/y (in this range M is assumed to increase linearly with v). Since the moment depends on cosq, the moment also depends on the dip angle of a slab.

(24) Although the sequence of phase transformations in a slab could be complicated involving olivine to modified spinel and modified spinel to spinel transformation, we consider olivine to spinel transformation only in this paper. Rubie and Ross [D. C. Rubie and C. R. Ross, II., Phys. Earth Planet. Inter. 86, 223 (1994)] argued that kinetic parameters for the olivine to modified spinel are similar to those for the olivine to spinel transformation. If this is correct, then the assumption of including only the olivine to spinel transformation will provide a reasonable approximation to the consequence of phase transformations in subducting slabs.

(25) M. R. Riedel and S. Karato, Geophys. J. Int. 125, 397 (1996); M. R. Riedel and S. Karato, Earth Planet. Sci. Lett. 148, 27 (1997). We also estimate the effects of grain-growth by incorporating the grain-growth kinetics. Grain-growth kinetics in spinel or modified spinel have not been investigated in detail. We use the data on olivine [S. Karato, Tectonophysics 168 , 255 (1989)] and estimate the grain-growth kinetics in spinel by using homologous temperature scaling. The preliminary data on modified spinel and spinel phase are largely consistent with this estimate (S. Karato, unpublished data).

(26) We assume ductile flow law for the most part of the slab. The contribution to the bending moment from the central elastic portion is negligibly small for curvatures less than ~ 1000 km (1). C. Goetze and B. Evans (1) assumed a constant strain-rate for the whole slab. However, this is inappropriate in our case, because we ask whether the strain-rate of a slab for a given force is high enough for deformation. The moment balance equation reads (21),

(A-1)

where M is the bending moment and u is the displacement of slab materials along the y-direction. Equation (A-1) can be solved to get (21),

(A-2)

where we used the relations

and

(L: the length of the slab) are used. We assume, for simplicity,

(for the definition of coordinate, see Fig. 2). Rheology and D are related through equation (2). This equation determines the strain-rate (and hence stress) for a given force (F) in a self-consistent fashion. When several deformation mechanisms operate, we choose a mechanism that gives the smallest effective viscosity at a given condition as a dominant mechanism.

(27) We consider rheology of olivine and spinel only. Rheology of modified spinel is assumed to be similar to that of spinel. For olivine rheology we used the data compiled by S. Karato and P. Wu [Science 260, 771 (1993)] in addition to S. Karato and D. C. Rubie [J. Geophys. Res. 102, 20,111 (1997)]. Rheology of spinel phase of silicate is poorly constrained. We used the data on germanate analogues [P. J. Vaughan and R. S. Coe, J. Geophys. Res. 86, 389 (1981); C. Dupas, H. W. Green, II., N. Doukhan, J-C. Doukhan and T. N. Tingle, Phys. Chem. Mineral. in press; J. D. Lawlis, Y. Zhang and S. Karato, EOS, 79, S331 (1998)] and assumed a homologous temperature scaling to estimate the pressure effects. For the pressure effects of melting temperatures, we used the data compiled by E. Ohtani [Phys. Earth Planet. Inter. 33, 12 (1983)]. Effects of water is not considered, because oceanic lithosphere is considered to have been depleted with water due to partial melting at mid-ocean ridges [S. Karato, Nature 319, 309 (1986); G. Hirth and D. L. Kohlstedt, Earth Planet. Sci. Lett. 144, 93 (1996); S. Karato and H. Jung, Earth Planet. Sci. Lett. 157, 193 (1998)].

(28) D. P. McKenzie, Geophys. J. R. astr. Soc. 18, 1 (1969).

(29) D. Olbertz, M. J. R. Wortel and U. Hansen, Geophys. Res. Lett. 24, 221 (1997); D. Olbertz, The Long-Term Evolution of Subduction Zones: A Modelling Study, Ph. D Thesis, University of Utrecht, pp. 152 (1997).

(30) S. H. Kirby, J. Geophys. Res. 92, 13789 (1987); M. Ogawa, J. Geophys. Res. 92, 13801 (1987); H. W. Green and P. C. Burnley, Nature 341, 733 (1989); S. H. Kirby, W. B. Durham, L. Stern, Science 252, 216 (1991); Kirby, Stein, Okal and Rubie (1996); M. Liu, J. Geophys. Res. 102, 5295 (1997).

(31) P. G. Silver, S. L. Beck, T. C. Wallace, C. Meade, S. C. Myers, D. E. James and R. Kuehnel, Science 268, 69 (1995).

(32) S. Karato, C. Dupas-Bruzek and D. C. Rubie, Nature in press.

(33) D. Yamazaki, T. Kato, E. Ohtani and M. Toriumi, Science 274, 2052 (1996).

(34) D. W. Forsyth and S. Uyeda, Geophys. J. R. astr. Soc. 43, 163 (1975).

Table 1. Classification of slabs and subduction parameters

Name type Vm* (cm/y) Vs (cm/y) Age (my) F (km)
N. Kurile A 1 7.3 120 8760 S. Kurile B 1 8.7 130 11300 NE Japan B 0.7 9.6 130 12480 Izu-Bonin B 3.7 8.7 146 12700 Mariana A 1.2 9.5 155 14725 Java A 0.5 7.7 96-134 7400-10300 Banda B - 8.1 150 12150 New Ginnea B 1.0 7.0 - - Tonga B 4.2 9.9 113 11200 Kermadec A 1 8.9 98 8720 S. Chille A 2.8 4.4 26 1140 C. Chille A 3.1 4.4 48 2110 N. Chille A 3.4 4.5 82 3690 Peru A 3.1 3.7 45 1670 C. America A 1.4 5.7 23 1311 N. America A 1.3 - 0-5 - Vm: trench migration velocity Vs: velocity of subduction F: slab thermal parameter (=<age>xVs) : Data of Vm, Vs are from the compile by D. Olbertz (29) and C. Lithgow-Bertelloni and M. A. Richards (13). Classification of slab deformation is based on the results of seismic tomography (3).

Table 2. Creep law parameters for olivine and spinel.

A generic creep law formula of

is assumed (Tm(P): melting temperature). The Peierls mechanism is important in the cold and/or high stress regions of a slab. Power-law creep is important in the high stress and/or moderate temperature regions. Diffusion creep plays an important role in the deep portions of a cold slab after a phase transformation where the grain-size is small.

A (s-1) n m g sP (GPa) q
olivine
Peierls 5.7x1011 0 0 31 8.5 2 power-law 3.5x1022 3.5 0 31 - 0 diffusion 8.7x1015 1 2.5 17 - 0
spinel
Peierls 7.0x1011 0 0 31 10 2 power-law 4x1022 3.5 0 31 - 0 diffusion 2x1016 1 2.5 17 - 0

FIGURES

Figure 1

Correlation of style of slab deformation with the slab thermal parameters and trench migration velocity. Slab geometry is classified into two categories (type A and B). The data are summarized in Table 1. Ages and kinematic data are the average values from the surface to the 660 km discontinuity. Portions of slabs that are connected are indicated by a tie line. Note a general correlation that type A behaviour (penetration of slabs into the lower mantle without much deformation) occurs mostly for slabs with small thermal parameters and Type B behaviour (significant deformation above the 660 km discontinuity) occurs mostly for slabs with large thermal parameters. Correlation of deformation behavior with trench migration velocity emphasized by previous studies (4, 10) appears rather secondary. An exception is the slab at the Mariana trench, which has the largest thermal parameter, yet penetrates into the lower mantle without much deformation. The reason for this exception may be its very steep dip angle (~90o) which results in a small bending moment.

Figure 2

Geometry of a subducting slab. A slab encounters a barrier at the 660 km discontinuity. The force due to this barrier tends to bend a slab. The degree of slab deformation depends on its rheological properties which depend on, among others, temperature, the magnitude of force (F) and subduction geometry such as the dip angle (

Figure 3

Dominant deformation mechanisms in subducting slabs. Initial temperature distribution corresponds to that for a 100 my oceanic lithosphere. The cases for velocity of subduction of 4 cm/y and 10 cm/y are shown. Because stress, temperature (and pressure) and grain-size change significantly in space for a given slab, dominant mechanisms of deformation change in a complicated fashion. In high stress, low temperature regions, the Peierls mechanism dominates. In moderate stress, moderate to large grain-size regions, power-law creep dominates. Diffusion creep plays an important role in fine-grain regions in the center of slabs after a phase transformation. Note that such a pattern also depends on the velocity of subduction, which controls the temperature distribution and the magnitude of stress.

Figure 4

Distribution of effective viscosity in subducting slabs corresponding to the cases shown in Fig. 3. Note that the center of slab has high viscosities for 4 cm/y, but the central portions in the deep part of a slab for 10 cm/y have low viscosity because of small grain-size caused by a phase transformation.

Figure 5

(a) Depth-flexural rigidity (D) relations for slabs as a function of velocity of subduction. Initial temperature distribution corresponding to the oceanic lithosphere of 100 my old age is considered. A broken line for 10 cm/y shows the results for the case where the effects of latent heat release and grain-size reduction are neglected. Flexural rigidity generally decreases with depth due to the increase in temperature. In addition, there is an abrupt decrease in flexural rigidity around 440 km due to the effect of latent heat release. Also, the production of fine-grained spinel in the deep portions of fast slabs results in significant reduction in flexural rigidity. Velocity of subduction affects flexural rigidity through two different mechanisms. First, it affects flexural rigidity through the change in temperature (fast slabs are cold). Usually low temperatures result in high strength, but when the degree of grain-size reduction is large, strength will be low for low temperatures due to a large effect of temperature on grain-size (24). Second, velocity of subduction changes stress levels that also affect rheology. The second effect is important where stress sensitive Peierls mechanism dominates. The effects of grain-size reduction and of the stress magnitude on rheology result in a negative correlation between flexural rigidity and velocity of subduction.

(b) The flexural rigidity of a subducted slab at 600 km deep as a function of velocity of subduction. The rigidity increases with velocity in ~4-7 cm/y because of the decrease in temperature, but beyond ~8 cm/y, the flexural rigidity decreases with velocity of subduction mainly because of the effects of temperature on grain-size.