Hans-Peter Bunge^*, Mark A. Richards, Carolina Lithgow-Bertelloni, John R. Baumgardner, Stephen P. Grand, Barbara A. Romanowicz

Computer models of mantle convection constrained by the history of Cenozoic and Mesozoic plate motions explain some deep-mantle structural heterogeneity imaged by seismic tomography, especially those related to subduction. They also reveal a 150-million-year time scale for generating thermal heterogeneity in the mantle, comparable to the record of plate motion reconstructions, so that the problem of unknown initial conditions can be overcome. The pattern of lowermost mantle structure at the core-mantle boundary is controlled by subduction history, although seismic tomography reveals intense large-scale hot (low-velocity) upwelling features not explicitly predicted by the models.

H.-P. Bunge, Institut de Physique du Globe de Paris, Laboratoire de Sismologie, 4 place Jussieu, 75252 Paris Cedex 05, France.
M. A. Richards, Department of Geology and Geophysics, University of California, Berkeley, CA 94720, USA.
C. Lithgow-Bertelloni, Department of Geological Sciences, University of Michigan, Ann Arbor, MI 48109, USA.
J. R. Baumgardner, Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87544, USA.
S. P. Grand, Department of Geological Sciences, University of Texas, Austin, TX 78713, USA.
B. A. Romanowicz, Berkeley Seismological Laboratory, and Department of Geology and Geophysics, University of California, Berkeley, CA 94720, USA.
^*   To whom correspondence should be addressed. E-mail: bunge@ipgp.jussieu.fr

Time-dependent Earth models are required to understand how the evolution of mantle flow affects Earth processes that occur on geologic time scales. For example, continental shelf and platform stratigraphy are controlled by vertical motions of the continental lithosphere in response to mantle convection (4). True polar wandering is caused by changes in the inertia tensor as a result of mantle convection (5), and the alternation between periods of rapid and slow magnetic field reversals is probably related to mantle-controlled changes at the core-mantle boundary (CMB).

The development of time-dependent Earth models has been delayed for several reasons: (i) Sufficient computer power to resolve the narrow thermal boundary layers in global mantle convection models has not been available; (ii) it is not obvious how the internal mantle density structure can be related to plate motion observations at the surface; and (iii) it is not known how time-dependent Earth models can be initialized at some starting point in the past, because the mantle density structure is known only for the present day (6).

Some of these difficulties have been overcome. (i) Advances in computer power allow three-dimensional (3D) spherical convection to be simulated at a resolution on the order of 50 to 100 km (7, 8). At the same time, large-scale mantle velocity heterogeneity structure has been mapped in greater detail (9, 10), and seismic tomography has imaged subducted slabs (11-13). (ii) The connection of internal mantle density structure to the history of subduction (14, 15) has allowed estimation of the internal buoyancy forces that drive plates (16). These developments allow convection models to be combined with plate motion reconstructions and such models to be tested with seismic data.

Figure 1B shows an Earth model obtained with the TERRA convection code (17, 18). More than 10 million finite elements provide an element resolution of about 50 km throughout the mantle, which allowed us to model convection at a Rayleigh number of 10⁸ (19). The history of plate motion is imposed as a time-dependent velocity boundary condition (20) starting in the mid-Mesozoic at 119 to 100 million years ago (Ma). We chose this starting time because well-constrained reconstructions exist only as far back in time as the 119 to 100 Ma period.

In computing the Earth model (Table 1) we assumed that (i) the mantle is of uniform chemical composition and (ii) is heated primarily from within by radioactivity, with a modest 20% component of bottom heating from an isothermal core. (iii) Viscosity increases by a factor of 40 in the lower mantle and in the lithosphere (21) relative to the asthenosphere. (iv) Mantle minerals undergo two phase transitions at a depth of 410 and 670 km. (v) The CMB supports no shear stresses. (vi) The motion of the upper boundary layer can be deduced from plate-motion reconstructions. (vii) The initial condition for thermal heterogeneity in the mantle is approximated by a quasi steady-state of convection derived by imposing mid-Mesozoic plate motions for all time before 119 Ma.

Table 1. Model parameters. Outer shell radius   6370 km Inner shell radius   3480 km T, surface   300 K T, CMB   2800 K  , surface   3500 kg m3  , CMB   5568 kg m3  , surface   4.011 × 105 K1  , CMB   1.256 × 105 K1   (upper mantle)   8.0 × 1021 Pa s   (lithosphere and lower mantle)   40 ×  upper mantle Thermal conductivity k   6.0 W m1 K1 Internal heating rate Qint   6.0 × 1012 W kg1 Heat capacity   1.134 × 103 J Kg1 K1  , 410 km   2 MPa K1  , 670 km   4 MPa K1

The first six assumptions roughly constitute a "standard model" for whole mantle convection (22). An internally heated mantle is consistent with cosmochemical and geochemical constraints on abundances of radioactive elements in the mantle (23), but a larger component of core heating is possible. Postglacial rebound and geoid modeling suggest an increase in viscosity of one to two orders of magnitude from the upper to the lower mantle (24, 25).

Our seventh assumption is the most problematic, so we tested how mantle convection responds to a realistic change in plate movement. We first imposed 119 to 100 Ma plate motions for 2 billion years to produce an initial condition (Fig. 1A). We then imposed present-day plate motions for 500 million years (My). The global cross-correlation of the evolving mantle temperature field with the initial condition starts from exactly one [perfect correlation (Fig. 2)] and evolves toward a final stage, at which nearly all initial-condition information is lost. Most of the adjustment occurs during the first 150 My, and the correlation falls to about 0.3. It reaches a value of about 0.2 in the remaining 350 My, which reflects the correlation of the Mesozoic and the present-day plate configuration. The CMB correlation also declines rapidly to about 0.4 during the initial 150 My, accelerated by our inclusion of about 20% core heating (18, 26).

The modeled response time of 150 My is comparable to the vertical transit time in our convection model, that is, the time it takes for thermal disturbances to be advected from the surface plates to the CMB. Thus, our simulation demonstrates that most of the initial-condition "information" from past plate-motion regimes is lost after about 150 Ma, suggesting that the plate motion record is probably just adequate for modeling the genesis of present-day mantle heterogeneity (27).

We compared the Earth model with the seismic shear body-wave velocity structure of Grand (12), which is similar to the compressional body-wave study of van der Hilst et al. (13), and with the long-period horizontal shear-wave (SH) study (SAW 12) of Li and Romanowicz (10). Convection at 1100 and 2000 km (Figs. 3A and 4A) is dominated by narrow, linear cold downwellings under North America and the western Pacific caused by subduction of the Farallon, Izanagi, and Pacific plates. Temperature variations away from downwellings are small, as expected for mostly internally heated convection. The Farallon slab is displaced eastward relative to the simple model of Ricard et al. (15), because of return flow from sub-Pacific into sub-Atlantic mantle beneath America.

Smaller scale downwellings away from plate bounaries are the result of boundary-layer instabilities beneath slow-moving plates such as Africa. However, their role is minor as a result of the relatively stiff lithosphere and the radial mantle viscosity structure (28). There is a lack of active hot upwellings, owing to the small amount of core heating. Near the CMB at 2800 km (Fig. 5A), heterogeneity is dominated by large-scale structure as the cold downwellings spread laterally at the CMB.

The body-wave tomography shows narrow, sheetlike downwellings corresponding to old subducted slabs beneath the Americas (Farallon plate), the Tethys, and the western Pacific. There is a prominent low-velocity anomaly under Africa. Other mantle regions show minor heterogeneity. The CMB is characterized by broadscale fast-velocity anomalies in the circum-Pacific and a very large amplitude low-velocity anomaly beneath Africa. The long-period SH study, which gives a particularly good fit to the nonhydrostatic geoid, agrees in the overall location of these anomalies.

Some characteristics of GEMLAB1 are similar to those of Grand's S-wave study, but there are important differences. For example, the cold subduction-related CMB temperature pattern in GEMLAB1 correlates well with the seismic models, but there are no prominent hot regions with temperatures substantially above the mean. This result is expected inside a hot thermal boundary layer and suggests that the very low wave-speed anomalies may not represent purely thermal effects. They may, however, result from chemical heterogeneity, which may also help to explain why some of the low wave-speed anomalies are confined to the deeper mantle.

In the midmantle, the strong endothermic phase change delays the sinking of material through the transition zone, and GEMLAB1 provides a relatively poor match to Grand's image of subducted slabs beneath some weaker subduction systems under South America and the Tethys. We compared GEMLAB1 to an Earth model GEMLAB2 (29) without the strong phase transitions, and found that GEMLAB2 provides a better fit to Grand's study, as evidenced by the spectral heterogeneity plots (Fig. 6).

The midmantle differences between GEMLAB1 and GEMLAB2 illustrate the effect of model parameter changes (phase transitions in this case), but uncertainties in the plate motion input data may be equally important. For example, our reconstructions indicate considerably more subduction beneath the northwestern Pacific than is evident in either seismic model, which suggest that reconstructions in GEMLAB models should be refined.

REFERENCES AND NOTES

1. B. H. Hager and R. J. O'Connell, J. Geophys. Res. 84, 1031 (1979).
2. A. M. Dziewonski, ibid. 89, 5929 (1984); J. H. Woodhouse and A. M. Dziewonski, ibid., p. 5953; B. H. Hager and R. W. Clayton, in Mantle Convection: Plate Tectonics and Global Dynamics, W. R. Peltier, Ed. (Gordon Breach, New York, 1989), vol. 4, chap. 9.
3. M. A. Richards and B. H. Hager, J. Geophys. Res. 89, 5987 (1984); Y. Ricard, L. Fleitout, C. Froidevaux, Annal. Geophys. 2, 267 (1984).
4. J. X. Mitrovica, C. Beaumont, G. T. Jarvis, Tectonics 8, 1079 (1989); M. Gurnis, Nature 344, 754 (1990) .
5. Y. Ricard, G. Spada, R. Sabadini, Geophys. J. Int. 113, 284 (1993).
6. Because heat transport involves diffusion, mantle convection models cannot be run backward in time from a relatively well-known present-day state, even if a perfect forward model could be constructed.
7. P. J. Tackley, D. J. Stevenson, G. A. Glatzmaier, G. Schubert, Nature 361, 699 (1993) .
8. H.-P. Bunge, M. A. Richards, J. R. Baumgardner, ibid. 379, 436 (1996) .
9. W.-J. Su, R. L. Woodward, A. M. Dziewonski, J. Geophys. Res. 99, 6945 (1994); G. Masters, S. Johnson, G. Laske, H. Bolton, Philos. Trans. R. Soc. London Ser. A 354, 1385 (1996).
10. X. D. Li and B. Romanowicz, J. Geophys. Res. 101, 22245 (1996).
11. S. P. Grand, ibid. 99, 11591 (1994).
12. ___, R. D. van der Hilst, S. Widiyantoro, GSA Today 7, 1 (1997).
13. R. D. van der Hilst, S. Widiyantoro, E. R. Engdahl, Nature 386, 578 (1997) .
14. M. A. Richards and D. C. Engebretson, ibid. 355, 437 (1992) .
15. Y. Ricard, M. A. Richards, C. Lithgow-Bertelloni, Y. Le Stunff, J. Geophys. Res. 98, 21895 (1993).
16. C. Lithgow-Bertelloni and M. A. Richards, Geophys. Res. Lett. 22, 1317 (1995).
17. H.-P. Bunge and J. R. Baumgardner, Comput. Phys. 9, 207 (1995).
18. H.-P. Bunge, M. A. Richards, J. R. Baumgardner, J. Geophys. Res. 102, 11991 (1997).
19. The Rayleigh number based on internal heating approaches the dynamic regime of the mantle, which probably convects with an effective Rayleigh number on the order of 10⁸ to 10⁹.
20. The plate tectonic information consists of 11 plate-motion stages extending back to 119 Ma. During a stage, plate motions are relatively constant. Stage boundaries correspond to episodes of large or sudden changes in plate motions, that is, plate rearrangements. Plate stage information consists of the complete plate boundaries for all the major plates and the Euler rotation vectors of each plate. The first six plate stages (0 to 10, 10 to 25, 25 to 43, 43 to 48, 48 to 56, and 56 to 64 Ma) span the Cenozoic and are based on the global plate boundaries and rotation poles of Gordon and Jurdy (30); the remaining five Mesozoic plate stages (64 to 74, 74 to 84, 84 to 94, 94 to 100, and 100 to 119) are based on the global compilation of Lithgow-Bertelloni and Richards (31). All plate boundaries and rotation poles are in the hotspot reference frame. The inherent uncertainties associated with using plate tectonic information in geodynamic models are described in (31). Many of the Mesozoic reconstructions, particularly for plates that have disappeared (Izanagi, Phoenix, Kula) or are in the process of disappearing (Farallon), are only approximate, especially the positions of subduction zones and to a lesser extent the rotation poles.
21. In the Earth the lithospheric viscosity is many orders of magnitude larger than the upper mantle viscosity, but the numerical capability to resolve such large viscosity contrast is not available. Simulations (28) demonstrate that a 20-fold viscosity increase in the lithosphere is sufficient to suppress boundary-layer instabilities.
22. G. F. Davies and M. A. Richards, J. Geol. 100, 151 (1992).
23. G. J. Wasserburg, G. J. F. MacDonald, F. Hoyle, W. A. Fowler, Science 143, 465 (1964) .
24. J. X. Mitrovica, J. Geophys. Res. 101, 555 (1996).
25. B. H. Hager and M. A. Richards, Philos. Trans. R. Soc. London Ser. A 328, 309 (1989).
26. This CMB response time is accelerated by a factor of 2, relative to a model without core heating.
27. The correspondence of the mantle response time and the time period for reliable reconstructions is not coincidental. Reconstructions are largely dependent on the magnetic isochrons of sea-floor spreading, which is limited to the characteristic maximum age for oceanic lithosphere. Simple boundary-layer convection theory (22) predicts that this characteristic boundary-layer age should be similar to the vertical transit time for convection. Based on the prediction that advection velocities scale roughly with the natural logarithm of the viscosity contrast (18), the mantle response time should be accurate within a factor of 2.
28. H.-P. Bunge and M. A. Richards, Geophys. Res. Lett. 23, 2987 (1996).
29. H.-P. Bunge et al., data not shown.
30. R. G. Gordon and D. M. Jurdy, J. Geophys. Res. 91, 12389 (1986).
31. C. Lithgow-Bertelloni and M. A. Richards, Rev. Geophys. 36, 72 (1998).
32. We thank G. Davies and R. van der Hilst for constructive reviews, J. Painter for supporting the 3D graphics, and the Los Alamos Branch of the Institute of Geophysics and Planetary Physics for continuing support. Computing resources were provided by the Advanced Computing Laboratory of Los Alamos National Laboratory.

Volume 280, Number 5360 Issue of 3 Apr 1998, pp. 91 - 95
©1998 by The American Association for the Advancement of Science.