Final Report

Simulations of the Earth's Core and Mantle Dynamics
(NCCS5-147)

Peter Olson, Johns Hopkins University, Principal Investigator

Contents

Mantle

A Brief History of TERRA Code Improvements during Round Two

Self-Consistent Generation of Plate Tectonics

Models of Plate Boundary Formation

Core

Report on Dynamo

Geodynamo Models

A Brief History of TERRA Code Improvements during Round Two

John Baumgardner, Los Alamos National Laboratory
 

Following the Round Two site visit held at UCLA in April 1996, an effort was launched to improve the performance of TERRA on the expected Round Two initial hardware, the Cray T3D. TERRA had been originally developed on and highly optimized for the Cray vector XMP and YMP machines. This effort focused on optimizing instead for the cache-based superscalar architecture of the T3D Dec Alpha processors. It involved primarily increasing the reuse of data in the very small cache on these processors. A large fraction of the arithmetic operations in TERRA involve application of 21-point stencils on a logically regular mesh. This character of the numerical algorithms provides the potential of significant data reuse, since field values at a given grid point are used several times as a stencil operator is applied at neighboring grid points.

The main strategy for achieving improved data reuse involved breaking single loops that swept over two spatial dimensions into two shorter loops that each sweep over a single spatial dimension. Whereas the T3D processor cache generally was too small to hold the two-dimensional arrays required for high reuse, it was found to be large enough for the smaller one-dimensional arrays. This strategy improved the per processor performance on the main kernels in TERRA from about 12 Mflops to the range of 35-40 Mflops -- approximately a three-fold increase. It required, however, a rewrite of all the important kernels in the code.

John Baumgardner at Los Alamos accomplished this work mostly during the month of May 1996 in collaboration with Mark Dalton of Cray Research at Los Alamos and Maynard Brandt in the Benchmark Services group at Cray Research in Eagan, Minnesota. In addition, Maynard Brandt rewrote the communication module to use T3D shmem routines in place of PVM library routines. The latter resulted in a two-fold reduction in the time spent in inter-process communication, such that it constituted only about six percent of the total.

Maynard Brandt at Cray Research made a quantitative measure of TERRA performance after these improvements. The test problem was a non-trivial one involving a fully 3-D temperature-dependent viscosity model using 16 T3D processors and a grid subdomain size identical to that planned for the 512 processor benchmark. From the CPU time and the number of floating point operations as determined by Apprentice, Brandt in early June 1996 obtained a TERRA performance of 22.0 Mflops/processor. This performance projected to a total of 11.24 Gflops on 512 T3D processors and indicated that TERRA's raw performance was now sufficient to achieve the first performance milestone of 10 Gflops.

The next hurdle was to remove a limitation in TERRA that restricted the number of processes to a power of four (i.e., 1, 4, 16, 64, 256, ...) such that, instead, any power of two (in particular, 9, corresponding to 512 processors) was possible. This involved modifying the manner in which domain decomposition was done and therefore also many of the details of the interprocess communication. It required generalization of all the communications routines to accommodate a variable number (either five or ten) of blocks of grid points per processor. (There are ten logically square blocks of grid points that fit together to cover the spherical surface in the icosahedral grid used in TERRA.)

A second similar hurdle was to eliminate a piece of serial coding at the lowest levels of the multigrid solver. This task involved mapping the multigrid calculations at the lowest grid levels, rather than to a single processor as had been done previously, instead, to one-fourth of the remaining active processors. This strategy allows the work at these lowest grid levels to be done in parallel instead of serially. But it required writing new communication routines to handle the new communication patterns as well as major modification of the data structures at these lower grid levels.

Both these tasks were completed by mid-January, 1997. Shortly thereafter TERRA met the 10 Gflops milestone with a performance of 11.6 Gflops on the NASA Goddard 512 processor Cray T3D. Speed per processor was found to be constant within about one percent for configurations of 16 to 512 processors when the total problem size was scaled to give roughly the same per processor work. This implied excellent parallel scaling.

With the arrival of the Cray T3E at NASA Goddard in the spring of 1997, effort was focused on optimizing the code for that platform. The most significant issue was making use of the new T3E hardware stream buffers that allowed more efficient transfer of data from main memory to secondary cache. Unfortunately, there was a design error in the boards that used the 300 MHz EV5 Alpha chip that could cause a system crash when the stream buffers were utilized. To work around this problem it was necessary initially to use versions of the communication routines based on MPI rather than shmem. Spencer Swift of Cray Research/SGI at NASA Goddard provided considerable help working through these difficulties as well as in modifying the kernel routines to achieve improved utilization of the stream buffers. With stream buffers activated, TERRA performance exceeded what was required to meet the 50 Gflops milestone by a comfortable margin.

In early August 1997 a performance of 38.0 Gflops was measured on the 512 PE 300 MHz T3E (T3E-600) at NASA Goddard (considered by NASA to represent half the resources required for the 50 Gflops milestone). At about the same time a performance of 86.8 Gflops was measured on a 1024 PE 450 MHz T3E (T3E-900) at Eagan, Minnesota (considered by NASA to represent 1.5 times the resources required for the 50 Gflops milestone). In the first case the performance was 38.0/25.0 = 1.52 times the milestone requirement, and in the second it was 86.8/75.0 = 1.16 times the requirement. In regard to the latter case, it is to be noted that although the peak processor speed of the T3D-900 exceeded that of the T3E-600 by 50%, the bandwidth to main memory was identical for both systems. Since the TERRA performance is limited mostly by memory bandwidth, it is understandable why its performance did not scale simply with peak processor speed. In any case, these performance numbers were readily satisfied the requirements of the second performance milestone.

The final performance milestone required that the TERRA calculation include its particle-based treatment of lithospheric plates. Lithospheric plates in TERRA are modeled as perfectly rigid patches in the topmost layer of nodes in the spherical grid. The spatial identity of the plates is carried by sets of particles, one set for each plate, that are moved in a Lagrangian manner across the surface. A single Euler rotation vector describes the motion of each plate. In a dynamical calculation the set of plate rotation vectors is obtained via a Newton iteration procedure that minimizes the net torque on each plate. The resulting piecewise uniform surface velocity field is then applied as a velocity boundary condition in solving the momentum equation for the velocity field in the interior of the spherical shell domain. A set of rules is applied at plate boundaries to treat the removal of particles at zones of convergence and addition of particles and creation of new plate in zones of divergence.

The particles in this treatment are always stored in memory in locations corresponding to the nearest grid point. This permits the basic particle data structure to be identical to that of the other grid based variables and allows identical domain decomposition and parallelization. It therefore allows relatively high parallel performance for the particle algorithms. Indeed we found that adding the particle treatment reduced the overall performance by only a few percent compared with runs without this treatment.

We achieved the third performance milestone, that of 100 Gflops, in November 1998 during a window of opportunity when a large Cray T3E-1200 was available at Eagan, Minnesota. TERRA achieved 121.8 Gflops on 1024 processors with the plate treatment included. Table 1 below summarizes the TERRA performance on three versions of the Cray T3E hardware with processor clock rates respectively of 300, 450, and 600 Mhz.

Table 1 TERRA Performance in billions of floating point operations per second (gigaflops) on parallel Cray T3E systems as measured in 1998. Work per processing element (PE) is approximately equal in all cases.
 

Illustrative TERRA Calculation

To illustrate TERRA's present capabilities, we present a model calculation for the Earth's mantle that exploits TERRA's ability to treat strong local variations in the strength of silicate rock that arise from variations in temperature and stress. Including this realism in the silicate deformation law leads to a remarkably Earth-like solution, with plate-like surface velocities, an internal density structure with strong harmonic degree two components, and an amount of toroidal energy in the surface velocity field similar to that actually observed. This solution varies only slowly with time after a transient response from an initial state involving high spatial frequency random temperature perturbations.

Our computational mesh, as displayed in Fig. 1, has 65 layers, with 163,842 grid points each, for a total of 10,649,730 grid points. The horizontal resolution at the surface is about 60 km. For simplicity in the interpretation of our results we make the Boussinesq approximation. We assume the spherical shell is heated only from within. We use Earth-like values for all the primary physical parameters except for the specific heat, which we lower by a factor of ten to reduce the Rayleigh number so as to match the resolution constraints imposed by the available computer resources. We also reduce by about a factor of five the activation energy in the Newtonian viscosity because of limitations on the viscosity gradients imposed by our algorithms. The depth dependence of viscosity, as shown in Fig. 2, is relatively simple. The upper mantle has a reference value of 2 x 1020 Pa-s, and the top 65 km and depths below 1000 km have reference values 1000 times larger, with smooth variation in the depth regions 65-250 km and 450-1000 km. We apply a viscous yield stress of 140 MPa in the top 250 km. The resulting total viscosity variation is approximately a factor of 200,000.

The calculation is initialized with random temperature variations imposed on a radial temperature profile that is uniformly 2000 K in the interior and drops to 300 K through a boundary layer at the surface. The temperature field initially displays a large number of point like downwellings, but these tend to coalesce into linear downwelling features which in turn coalesce to form a single linear network as displayed in Fig. 3(a). The time scale for this development is about 109 years. Thereafter the pattern retains its topological character and changes only slowly with time as indicated in Fig. 3(b). Figures 4 and 5 show the viscosity field with large high viscosity patches at the surface characterized by low internal rates of strain. Weak zones that display high rates of deformation and represent zones of localized downwelling and upwelling surround these platelike patches. Numerous instances of oblique motion at downwelling zones lead to a toroidal component of velocity that is 30-35% of the poloidal component in magnitude at the surface. The geoid is dominated by its harmonic degree two components.

From a scientific standpoint, this solution which displays such remarkable similarity to the Earth leads us to conclude much of the mantle's unique behavior and history can be traced directly to a few crucial but simple aspects of its rheology.

This case was run on 128 processors of the Cray T3E-600 `jsimpson' at NASA's Goddard Space Flight Center. A typical time step of length 40,000 years required 6.5 seconds of CPU time on 128 processors. This corresponds to 163 CPU seconds per Myr of problem time or 45 CPU hours per Gyr. Each time step typically involves four multigrid cycles to solve for 42.6 million unknowns (three velocity components plus pressure at each grid point) plus an update of the energy equation. The total execution rate for this case on the 128 processors of the T3E-600 is approximately 11 Gflops. TERRA also displays a similar level of performance on current Beowulf clusters of similar size.
 

Figure 1 Cutout view of a portion of the TERRA computational grid. Horizontal grid is constructed by successive refinements of the projection of the regular icosahedron onto the sphere. Replicating this horizontal grid at appropriate radial locations generates the 3D grid. The grid pictured here has 65 radial layers, each with 163,842 points, for a total of 10,649,730 grid points.

Figure 2 Reference radial viscosity profile. Lateral viscosity variation, arising from lateral temperature variations and yield stress weakening, multiplies this radial variation to yield the 3D viscosity field.

Figure 3a

 

Figure 3b

Figure 3 Temperature and velocity fields at 367 km depth and problem times of 1.0 Gyr (a) and 2.0 Gyr (b). Note the linear pattern of downwelling flow that resembles, at least in spatial frequency content, the pattern of subduction observed today for the Earth. Arrows denote horizontal velocity, triangles upward radial velocity, and squares downward radial velocity. Velocity magnitude is proportional to symbol size.

Figure 4a

Figure 4b

Figure 4 Viscosity and velocity fields at problem time of 2.0 Gyr. (a) Fields at surface. (b) Fields at 367 km depth. In (a) note regions of reduced viscosity at zones of flow convergence and divergence with regions of high viscosity elsewhere. In (b) the high viscosity in the narrow zones of downwelling flow is a consequence of the temperature dependence of the viscosity law.

Figure 5 Cutout view of 3D viscosity field with a high viscosity isosurface exposed in the cutout region. Top corner of cutout is at the north pole, and central longitude of cutout is approximately 115 degrees east relative to the central meridian of Figs. 3 and 4. Note the high viscosity platelike patches (in orange) at the outer surface, the low viscosity zone (in blue) corresponding to the Earth's upper mantle, and the contorted extensions of the high viscosity downwelling flow in the lower portion of the spherical shell corresponding to the Earth's lower mantle.
 

Related Publications

W.-S. Yang and J.R. Baumgardner, "Matrix-dependent transfer multigrid method for strongly variable viscosity infinite Prandtl number thermal convection," Geophys. and Astrophys. Fluid Dyn., in press, 2000. 

H. R. Wenk, J. R. Baumgardner, C. N. Tome, and R. Lebensohn, "A deformation model to explain anisotropy of the inner core," J. Geophys. Res., in press, 2000. 

M. A. Richards, H.-P. Bunge, Y. Ricard, and J.R. Baumgardner, "Polar wandering and inertial interchange events in mantle convection models," Geophys. Res. Lett., 26, 1777-1780, 1999. 

M.A. Richards, H.-P. Bunge, C. Lithgow-Bertelloni, and J.R. Baumgardner, "Mantle convection and plate motion history: Toward general circulation models," History and Dynamics of Global Plate Motions, AGU Monograph Series, 1999. 

C.C. Reese, V. S. Solomatov, J. Baumgardner and W. Yang, Stagnant lid convection  in a spherical shell, Phys. Earth Planet. Int., 116, 1-7, 1999.

H.P. Bunge, M.A. Richards, C. Lithgow-Bertelloni, J.R. Baumgardner, S.P. Grand, and B.A. Romanowicz, Time scales and heterogeneous structure in geodynamic Earth models,  Science, 280, 91-95, 1998.

H.P. Bunge, M.A. Richards, and J.R. Baumgardner, A sensitivity study of 3-D spherical mantle convection at 10^8 Rayleigh number, J. Geophys. Res., 102 (B6), 11991-12007, 1997. 

H.P. Bunge, M.A. Richards, The origin of large scale structure in mantle gonvection, Geophys. Res. Lett., 23, 2987-2990, 1996. 

H.P. Bunge, M.A. Richards, and J.R. Baumgardner, The effect of depth-dependent viscosity on the planform of mantle convection, Nature, 379, 436-438, 1996. 

Self-Consistent Generation of Plate Tectonics

Paul Tackley, UCLA

One of the fundamental goals of this project, and of mantle convection research in general, is to identify  a material description which will allow tectonic plates to form naturally in mantle convection simulations, rather than being imposed by the modeler as has been necessary in the past. This should help to answer fundamental  questions such as "Why does the Earth have plate tectonics, whereas Venus and Mars do not". This problem is particularly challenging in three-dimensional geometry due to the presence of transform plate boundaries, which do no arise naturally in a convecting system.

In summer 1997, Tackley obtained for the first time, self-consistently generated plates  in a three-dimensional mantle flow  calculation,  work which was subsequently published in Earth and Planetary Science Letters. They key component of the material properties  was strain-rate weakening,  i.e., decreasing stress with increasing strain rate. Strain-rate weakening had previously been investigated in a two-dimensional sheet by co-PI Bercovici. The material and numerical model used for these initial results was quite simple, and as a result the plates were not fully realistic. In particular, the calculations considered only the instantaneous flow for a given distribution of buoyancy  forces, rather than being time-dependent, and assumed a fixed-thickness lithosphere, rather than one which arose due to the temperature-dependence of viscosity in the cold upper boundar layer .

To overcome these limitations, Tackley developed a rheological parameterization which combines temperature- depth- stress- and history- dependence. The history-dependence is described through the time-dependent evolution of 'damage', which weakens the material. Damage arises where material is deformed at high stress, and heals with time at a rate which is temperature-dependent, allowing long-term memory at cold, lithospheric temperatures but rapid healing at hot mantle temperatures. Under steady-state conditions, the damage evolution reduces to the previously-used strain-rate weaking. This description was shown to lead to the generation of weak transform margins, suduction zones and spreading centers in steady-state models in a paper in press in an AGUGeophysical Monograph volume.

Fully time-dependent, self-consistent models followed in fall 1998, when Tackley developed the first fully time-dependent, self-consistent 3-D models which develop plate-like behavior evolving continuously in space and in time. These results were presented at several meetings, and will be submitted for publication shortly. An illustrative result shown in the Figure below shows passive rifting, passive oceanic spreading centers, single-sided subduction, and back-arc spreading , all characteristic features of Earth's plate tectonics.  In these results, the plateness is, to first order, caused by a pseudo-plastic yield stress rather than strain-weakening or strain-rate weakening. Having a low-viscosity asthenosphere under the lithosphere results in sharper plate boundaries but is not necessary for plates to exist. A systematic investigation of these models is under way.

Figure 6  Self-consistent plate generation, including continental plates (single-sided subduction)

An important component of Earth's plate tectonics are the compositionally-distinct continental plate regions. The interaction between buoyant continental material and mantle flowhas been studied by Tackley and student Li Ma using three-dimensional simulations with buoyant continental material at the top. The simulations are distinguished from previous works not only by dimensionality but also by treating the continents using a self-consistent rheological description allowing them to deform, break up etc., as opposed to the traditional 'rigid block' approach. Buoyant continental material has also been included in the fully self-consistent plate generation models, as shown above.

Elasticity and related brittle failure are undoubtedly important in the shallow lithosphere and may be responsible for some of the characteristic features of plate tectonics and instrumental in the formation of new plate boundaries. Progress was made by project postoc Moritz Heimpel on the inclusion of elasticity and brittle failure in the STAG convection code and it is expected that this will be completed in the near future.
 

Mantle-Core Coupling

a) Coupled mantle and core dynamics on Io. Jupiter's moon Io is a prime example of a body for which mantle-core coupling may play an important  role. Recent measurements by the Galileo spacecraft indicate that Io has an active component to its magnetic field. Circulation in the core is probably driven by differential heating and cooling from the mantle. Io's mantle is particularly interesting, being unique amonst planetary mantles so far considered because convection is driven by non-uniform tidal dissipation rather than uniform radiogenic heating or an isothermal hot lower boundary. Thus, Tackley in collaboration with Schubert and Glatzmaier have been simulating mantle convection in Io, using the spectral MAN code. Results to date indicate that the style of convection consistens of a mean flow driven by the distribution of tidal heating, with small-scale instabilities superimposed on the mean flow. As Rayleigh number is increased these small-scale instabilities increasingly distribute the heat, so that the surface distribution of heat flow is much more uniform than the tidal dissipation distribution would suggest. Patterns of heat flux across the core-mantle boundary are being used as a boundary condition for core dyamo simulations by Glatzmaier and Roberts.
 

b) The core-mantle-boundary region: thermo-chemical convection calculations with a dense D" layer. The D" region, lying at the base of the mantle directly above the core, has received a lot of interest and attention in recent years and is a key area for core-mantle interactions. One possibility for D" is that it consists of a layer of chemically dense material. Tackley has performed the first three-dimensional calculations of mantle convection with a dense layer at the base, in order to characterize the basic planforms and dynamics. Dense material is swept up into 'piles' underlying the upwellings, and slowly entrained. Entrainment is always in cylindrical forms, even in completely internally-heated convection. If depth-dependent properties (such as viscosity) result in very sluggish convection in the deep mantle (as indicated by a low 'local' Rayleigh number) then these piles can be very large, even if the chemical density contrast is small. These large thermo-chemical piles or 'megaplumes' may provide an explanation for features seen in seismic tomographic models under Africa and the Pacific.

Publications

Tackley, P. J., Self-consistent generation of tectonic plates in three-dimensional mantle convection,  Earth Planet. Sci. Lett., 157, 9-22, 1998.

Tackley, P. J., Three-dimensional simulations of mantle convection with a thermo-chemical basal boudary layer: D"?, in The Core-Mantle Boundary Region, AGU Geodynamics Series 28, Ed. Gurnis et al., pp231-253, 1998.

Tackley, P. J., Baumgardner, J. R., Glatzmaier, G. A., Olson, P., Clune, T., Three-dimensional spherical simulations of convection in Earth's mantle and core using massively-parallel computers, Proceedings of High Performance Computing 1999, pp. 95-100, Ed. A. Tentner, 1999.

Tackley, P.J., The quest for self-consistent generation of plate tectonics in mantle convection, in The History and Dynamics of Global Plate Motions, AGU Geodynamics Series, Ed. M.A. Richards et al., in press, 1999.

Tackley, P.J., Three-dimensional, time-dependent  models of mantle convection with self-consistent plate tectonics, in preparation for Geochemisty, Geophysics, Geosystems, 1999.

Tackley, P.J., G. Schubert, G.A. Glatzmaier, and J.T. Ratcliff, Three-dimensional simulations of mantle convection in Io, in  preparation for Icarus, 1999.
 

Meeting Presentations (reverse chronological order)

"The coupled system of mantle convection and lithospheric plate tectonics: towards a self-consistent model", P.J. Tackley, IUGG meeting, Birmingham, 1999

"Three-dimensional spherical simulations of mantle convection in Io", P.J. Tackley, G. Schubert, G.A. Glatzmaier, J.T. Ratcliff, IUGG meeting, Birmingham, 1999

"Dynamical constraints on the thermo-chemical megaplumes in the deep mantle beneath Africa and Hawaii", P.J. Tackley, Spring AGU 1999.

"Self-consistent generation of plate tectonics in time-dependent three-dimensional mantle convection", P.J. Tackley, Spring AGU 1999.

"Three-dimensional spherical simulations of convection in Earth's mantle and core using massively-parallel computers", P.J. Tackley, J. R. Baumgardner, G. A. Glatzmaier, P. Olson, T. Clune, At the High Performance Computing ASTC conference, San Diego 1999.

"Large-scale thermo-chemical structures in the deep mantle: Constraints from 3-D numerical simulations of mantle convection and seismology", P.J. Tackley, Fall AGU 1998.

"Self-consistent generation of plate tectonics from mantle convection, Part II: Memory and time-dependent evolution", P.J. Tackley, Fall AGU 1998.

"Three-dimensional numerical simulations of mantle convection in Io", G. Schubert, P.J. Tackley, J.-P. Matas, G.A. Glatzmaier, J.T. Ratcliff, Fall AGU 1998.

 "Towards self-consistent generation of plate tectonics in 3D mantle convection", P.J. Tackley, Gordon conference NH 6/98

 "Advances in three-dimensional mantle convection modeling: variable viscosity, plates, and geochemistry", P.J. Tackley, SEDI Tours, France 7/98

 "Three-dimensional simulations of mantle and core convection on Io", SEDI Tours, France 7/98

"Continental Drift, Aggregation, Breakup and Dispersal in a 3-D Mantle Convection Model", P.J. Tackley and L. Ma, AGU Chapman conference on The History and Dynamics of Global Plate Motions, 9/97.=20

 "Self-consistent Generation of Plate Tectonics in 3-D Mantle Convection",  Tackley, P.J., AGU Fall Meeting, 1997.

"Continental Drift, Aggregation, Breakup and Dispersal in a 3-D mantle convection model", Tackley, P.J., EOS Trans. AGU Fall Meeting, 1997.

Models of Plate Boundary Formation

David Bercovici, University of Hawaii
 

Figure 7 Approximate horizontal divergence (top) and vertical vorticity (middle) fields of the Earth (plate and continental outlines are shown for reference in the bottom frame).Horizontal divergence is a fine-scaled representation of poloidal flow at the Earth's surface and represents the so-called rate of creation (divergence) and destruction (convergence or negati e divergence) of the plates at ridges and subduction zones, respectively. Vertical vorticity is  a fine-scaled representation of toroidal flow and represents the rate of strike-slip shear between plates as well as plate spin (vorticity is essentially a measure of the local angular velocity, i.e., angular velocity of a point in the fluid). The divergence and vorticity fields shown are estimated directly from a model that uses finite-width plate margins based on seismicity distributions; for the standard plate model with infinitesimall thin margins, these fields are singularities and thus cannot be reliably estimated.

Figure 8 Generation of the Earth's toroidal motion via a source-sink model. The Earth's horizontal divergence (Fig. a -- top) (a simpler and smoother version of that shown in the previous figure)  is used as  a source-sink field to drive flow in a non-Newtonian thin (lithospheric) fluid shell.  The divergent zones (yellow and red) are used as sources of material, while the convergent zones (blue) are used as sinks. Ideally, the interaction of the source-sink  field and non-Newtonian rheology in the lithospheric fluid shell would recover the Earth's vertical vorticity field (Fig. (a) -- bottom) (again, simpler and smoother than that shown in \plateref{divvor}). Flow in the models with power-law rheologies do not generate Earth-like vorticity (Fig. (b) -- top and middle frames) for power-law indices that are mantle-like, n=3, and extreme, n=21.  Flow with the continuous stick-slip or self-lubricating rheologies (labelled as n=-1) (Fig. (b) -- bottom) generates a more Earth-like vorticity field.  See Bercovici [1995a] for further discussion.

Figure 9 A simple model of a thin lithospheric-type layer of fluid driven by a basic source-sink field S (top row; yellow and red represent a source of material, blues represent sinks). The viscosity of the fluid layer is variable. In one model the viscosity is a function of temperature Theta and heat is generated by viscous dissipation (middle row showing, left to right, dimensionless horizontal velocity v_h, vertical vorticity or strike-slip shear omega_z [see also Figs 7 and 8], and temperature Theta). In the second model viscosity is a function of volatile (e.g., water) content measured in terms of volume  raction of voids or porosity Phi; the voids are generated according to the amount that the fluid is stressed and damaged (bottom row, showing left to right, velocity, vorticity and porosity). The fields S, omega_z, Theta and Phi as shown are multiplied by a factor of 1000. The interaction of the source-sink field and the variable viscosity leads to toroidal flow. Toroidal motion (represented with vertical vorticity) and plate strength  or  plateness (measured by temperature or porosity distribution, where hot or porous is weak, cold and unporous is strong) remains diffuse and unplate-like in the viscous-heating based mechanism even with very temperature-dependent viscosity. However, toroidal motion in the void-volatile model can allow extremely focussed, nearly fault-like concentrations of strike-slip shear or toroidal motion (bottom middle frame) and plate- ike strength distributions (bottom right frame). (After Bercovici [1998].)

Relevant Publications

Bercovici, D., Y. Ricard and G. Schubert, A two-phase model for compaction and damage, Part 1: General theory, submitted to  J. Geophys. Res. 1999.

Ricard, Y. Bercovici, D., and G. Schubert, A two-phase model for compaction and damage,  Part 2: Applications to compaction, deformation and and the role of interfacial surface tension, submitted to  J. Geophys. Res. 1999.

Bercovici, D., Y. Ricard and G. Schubert, A two-phase model for compaction and damage,  Part 3: Applications to shear localization and plate boundary formation, submitted to  J. Geophys. Res. 1999.

Bercovici, D., Y. Ricard, and M.A. Richards, The relation between mantle dynamics and plate tectonics: A primer, in press AGU Monograph Series ``History and Dynamics of Global Plate Motions", M.A. Richards, R. Gordon and R. Van der Hilst, editors, 1999.

Bercovici, D., Generation of plate tectonics from  lithosphere-mantle flow and void-volatile self-lubrication,  Earth Planet. Sci. Lett., 154, 139-151, 1998.

Dumoulin, C., D. Bercovici,  and P. Wessel, A continuous  plate-tectonic model using geophysical data to estimate plate margin widths, with a seismicity based example, Geophys. J. Int., 133, 379--389, 1998.

Bercovici, D., Plate generation in  a simple model of lithosphere-mantle flow with dynamic self-lubrication, Earth Planet. Sci. Lett., 144, 41--51, 1996.

Bercovici, D., On the purpose of toroidal flow in a convecting mantle,  Geophys. Res. Lett.,  22, 3107--3110, 1995b.

Bercovici, D., A source-sink model of the generation of plate tectonics from non-Newtonian mantle flow,  J. Geophys. Res., 100, 2013--2030,  1995a.

Bercovici, D. and P. Wessel, A continuous kinematic model of plate tectonic motions, Geophys. J. Int.,  119, 595--610, 1994.

Bercovici, D., A simple model of plate generation from mantle flow,  Geophys. J. Int., 114, 635--650, 1993.

Gary A. Glatzmaier, Earth Sciences Department, University of California, Santa Cruz

Thomas T. Clune, Silicon Graphics Inc.

Paul H. Roberts, Institute of Geophysics and Planetary Physics, University of California, Los Angeles

A computer simulation of magnetic field generation by convection in the Earth's fluid outer core produces a time dependent field that, like the Earth's, is strongly dipole dominated outside the core and occasionally undergoes spontaneous dipole reversals.

The 3D spherical model solves the magnetohydrodynamic equations within the anelastic approximation. Because of limited spatial resolution, the viscous and thermal diffusivities crudely represent the transport of momentum and heat by unresolved eddies. An Earth-like heat flux imposed over the core-mantle boundary (CMB) controls the cooling rate of the core, which determines the vigor of the convection and intensity of the generated magnetic field. The Earth's rotation rate is essential for producing the sheared and helical fluid flows for the dynamo; it is also what makes the numerical method, which treats the Coriolis force implicitly, so challenging.

The DYNAMO code uses a spectral transform formulation based upon expansion in  spherical harmonics and Chebyshev polynomials.  The equations are integrated forward in time using a third-order Runge-Kutta scheme; nonlinear terms are treated explicitly and linear terms are treated implicitly.

The high-performance parallel code is based on a message passing paradigm. In the broadest terms, the parallelization strategy has been to perform the spectral transforms on data that is ``in processor'' while distributing the orthogonal directions among the processors in a load-balanced manner. Global memory transposes are then used to rearrange the data as necessary for each type of operation. The spectral transforms themselves are highly efficient operations (in terms of hardware utilization) on virtually all processors, but particularly so on cache-based hierarchical memory systems typical of modern parallel computers.

As part of NASA's Earth and Space Science High Performance Computing Cooperative Agreement, DYNAMO has been highly optimized for the Cray T3E. The T3E possesses an extremely low-latency, high-bandwidth network which is ideal for the global memory transposes used in spectral transform methods. Every effort has been made to maximize cache-reuse, reduce the number of communication packets, and to exploit the specialized  E-registers for local and remote memory transfers. Because the computational load of the Legendre transforms is dominant for all but very low resolution cases, the time spent in transforms also dominates the time spent communicating, thereby enabling excellent scalability. DYNAMO has run on as many as 1490 processors and attained performance in excess of 600 billion floating point operations per second (GFLOPS). For moderate grid resolutions suitable for our current research, DYNAMO  typically achieves between 100 and 200 GFLOPS on 512 processors of the T3E.

Several numerical simulations have been obtained with the Glatzmaier-Roberts model; they span over a million years using a numerical time step of 15 days or less.  Here we describe a medium resolution case with 49 radial levels (plus 17 in the inner core), 144 latitudinal levels, and 144 longitudinal levels. The magnetic energy, integrated through the core, is more than three orders of magnitude greater than the kinetic energy of the convection that maintains it.  Maximum magnetic field intensity, in the deep interior of the core, can be 20 mT and maximum fluid velocity is typically about a millimeter per second.  Although this is much less than, for example, atmospheric wind speeds, it is a million times greater than mantle convection.  Convection of the dense metallic core fluid is so efficient that the maximum temperature variation (on a sphere of constant radius) is only about 0.001K.

Figure 10 shows a typical magnetic energy spectrum (in spherical harmonic degrees) for the field at the CMB for the Earth (out to degree 12) and for the Glatzmaier-Roberts simulation (to degree 95).  (The Earth's core field is known only out to degree 12 because the crustal magnetic field significantly interferes for larger degrees.)  Both spectra show the dipole (l=1) component dominating.  Since the dipole part of the field decreases the least rapidly with distance from the core, the structure appears even more dipole dominated at the Earth's surface as i lustrated in Figure 2 where the radial component of the field is plotted in equal area projections.  One can also see in Figure 2 that although including degrees 13 to 95 produces no detectable difference in the surface field, it makes a significant difference in the structure of the field at the CMB.  This argues for much more small scale structure in the geomagnetic field at the CMB than is seen in the current degree 12 spectrum of the geomagnetic field. The intense concentrations of magnetic flux (core spots) seen in the simulated field at the CMB are somewhat similar to sun spots on the solar surface.

Figure 10 A snapshot of the magnetic energy density at the core-mantle boundary for today's Earth (solid circles, up to degree 12) and for the Glatzmaier-Roberts geodynamo simulation (open circles, up to degree 95).

Figure 11 The radial component of the magnetic field at the surface and at the core-mantle boundary for the Earth (up to degree 12) and for the Glatzmaier-Roberts simulation (up to degree 12 and up to degree 95). The surface field has been multiplied by 10 to give comparable color contrast.

Recent Publications

Glatzmaier, G.A., Coe, R.S., Hongre, L. and Roberts, P.H. The role of the Earth's mantle in controlling the frequency of geomagnetic reversals. Nature 401, 885-890 (1999).

Clune, T.C., Elliott, J.R., Miesch, M.S., Toomre, J. and Glatzmaier, G.A. Computational aspects of a code to study rotating turbulent convection in spherical shells.  Parallel Computing 25, 361-380 (1999).

Roberts, P.H. and Glatzmaier, G.A. Theory and simulations of the geodynamo. Rev. Mod. Phys., submitted (1999).

Roberts, P.H. and Glatzmaier, G.A. A test of the frozen flux approximation using geodynamo simulations. Phil. Trans. Roy. Soc., submitted (1999).

Glatzmaier, G.A. and Clune, T. Computational aspects of geodynamo simulations. Computing in Science and Engineering, submitted (1999).

Peter Olson, Johns Hopkins University

I. Technical Accomplishments:

In year 3 of the project we created and used a new workstation version of the DYNAMO code, along with user documentation and a graphical interface. The workstation version is based on the Boussinesq (incompressible fluid limit) version of the full code. It uses an improved explicit time step procedure, in which the Alfven waves in the boundary layers are numerically damped. This permits longer time steps, and allows for longer simulation times.  The code is readily portable to any workstation or PC running UNIX or LINUX.  The graphical interface is written in IDL (Interactive Data Language) and can be used on UNIX, Windows, or Mac platforms.

II. Scientific Accomplishments: A Systematic Parameter Study of Numerical Dynamos

We have analyzed about 50 3D numerical dynamos driven by thermal convection in a rotating spherical shell with the geometry of the Earth's core (Olson, Christensen and Glatzmaier, 1999; Christensen, Olson and Glatzmaier, 1999). We examined the influence of rigid versus stress-free boundaries and, by varying the Ekman, Rayleigh, and Prandtl numbers, different combination of magnetic, thermal and viscous diffusivities. We find mostly dipole-dominated dynamo solutions, except at the highest Rayleigh numbers, where the dynamos are sometimes octupole  dominated (see Figure 12 below).  When low-pass filtered to mimic the filtering effect of the mantle on the magnetic field, several of the models exhibit magnetic field strucutures that resemble the geomagnetic field at the core-mantle boundary.
 

  Figure 12 A comparison of numerical dynamos driven by thermal convection in rotating, electrically conducting fluid spheres. The left image in each panel shows contours of the radial component of the magnetic field on the outer boundary  (radius R). The right image in each panel shows contours of the radial component of the fluid velocity at radius 0.8R. Red colors indicate positive values and blue colors indicate negative values.  a: quasi-stationary dipolar dynamo with azimuthally periodic convection at Ekman number E=0.001 and Rayleigh number 1.8 times critical. b: nondipolar dynamo at E=0.001 and Rayleigh number 14 times critical. c: dipolar dynamo at E=0.0001 and Rayleigh number 3.5 times critical. d: dipolar dynamo at E=0.0001 and Rayleigh number 15 times critical. e: dipolar dynamo at E=0.0001, Rayleigh number 11 times critical and  magnetic Prandtl number 2. f: dipolar dynamo at E=0.0001, Rayleigh number 11 times critical and  magnetic Prandtl number 0.5. Among these, d and f are most similar to the geodynamo.

Publications Supported by NASA HPCC (Year 3 - 1999)

Olson, P., U. Christensen, and G.A. Glatzmaier, Numerical modeling of the geodynamo:  Mechanisms of field generation and equilibration, J. Geophys. Res., 104, 10.383-10,404, 1999.

Christensen, U., P. Olson and G.A. Glatzmaier, Numerical modelling of the geodynamo:  a systematic parameter study, Geophys. J. Int., 138, 393-409, 1999.

Tackely, P., J.R. Baumgardner, G.A. Glatzmaier, P. Olson and T. Clume, Three-dimensional spherical simulations of convection in Earth's mantle and core using massively-parallel  computers, Proc. High. Performance Computing Symposium-HPC'99 95-100, 1999.

Olson, P. and J. Aurnou, A polar vortex in the Earth's core, Nature, 402, 170-173, 1999.