Louis Moresi’s article in the December issue of SIAM News prompts readers to take a deeper look at models of Earth – deeper into Earth, that is, and more specifically into Earth’s mantle: the region between the solid crust and the metallic core.
The mantle makes up about 84% of Earth by volume and extends to a depth of 2,890 km. At its bottom, it is believed to reach temperatures of about 3,300–4,400 K. Yet because of pressure from overlying material, the mantle consists of solid rock (with few localized exceptions of partially molten material), as evidenced by its ability to transmit seismic shear waves. However, on sufficiently large time scales it behaves like a vigorously convecting fluid, cooled by the cold surface above and heated by both the hot core below and radioactive decay within (see Figure 1). Thus, one could view plate tectonics simply as the surface expression of convection cells. Within the mantle, average velocities are on the order of cm/year.
Figure 1. Schematic of convection in Earth’s mantle, including cold downwellings (blue), hot upwellings (red), and strong chemical heterogeneities (orange, green), all of which influence the material properties of the rocks.
At first glance there seem to be only a few visible expressions of mantle convection. However, most volcanism on Earth’s surface is closely related to the distribution and advection of thermal and chemical heterogeneities in the mantle and at mid-ocean ridges, subduction zones, and oceanic islands, for example. Furthermore, the importance of the mantle lies in its interactions with other parts of Earth:
• When hot material rising in the mantle approaches the surface, it causes massive melting that leads to large volcanic eruptions. The volume of erupted material can exceed 106 km3 in just a few million years, leading to the formation of large volcanic provinces with areal extensions greater than 106 km2. Due to the released gases, these events influence the global climate and may have led to mass extinction events.
• Cooling by the mantle drives convection in Earth’s core, which generates Earth’s magnetic field. A more viscous mantle would effectively insulate the core, and a less viscous mantle would have led to the core’s crystallization long ago; both scenarios would leave us without the magnetic field that offers protection from the harsh radiation of space and keeps solar wind from driving away the atmosphere.
• By ingesting water and sediment-laden subducting oceanic plates—and later releasing material again through volcanic outgassing—the mantle acts as a giant reservoir for water and carbon, thus creating a carbon cycle on time scales of tens of millions of years.
For these reasons, understanding the dynamics of convection and the conditions under which it happens have long been important in the geosciences. More recently, a fair number of computational and applied mathematicians have also entered the field, as mantle convection presents an attractive, accessible set of problems that require large-scale computations, are tough on linear and nonlinear solvers, yet still fall within a solvable range.
Many aspects of mantle convection can already be understood qualitatively by considering a set of equations in which fluid flow is slow and driven by density differences that result from temperature variation. The Stokes equations allow for adequate modeling of velocity and pressure:
\[-\bigtriangledown \: \cdot \: (2\eta \varepsilon (u)) + \bigtriangledown p = \rho g,\]
\[ \bigtriangledown \: \cdot \:u = 0,\]
with viscosity (ranging from 1018 to 1023 Pa·s; for comparison, the viscosity of water is about 10-3 Pa·s; that of air 10-5 Pa·s; and that of honey ~10 Pa·s), density \(\rho\)and gravity \(g\). An advection diffusion equation models temperature:
\[\frac{\partial T}{\partial t} + u \cdot \: \bigtriangledown T - \bigtriangledown \: \cdot \: \kappa \bigtriangledown T=Q.\]
This must be augmented with appropriate boundary conditions for the velocity and temperature, such as tangential or zero velocity around the boundary, and “hot” at the bottom and “cold” at the top of the domain. In the simplest description, the density depends linearly on the temperature, \(\rho (T) = \rho_{ref} - \alpha (T - T_{ref}\!)\), where \(\alpha\) is the thermal expansion coefficient. With realistically small values for the thermal diffusivity \(\kappa\), such a model already (correctly) predicts that heat is transported from the core to the surface primarily through “blobs,” or sheets of hot material that rise up from the core-mantle boundary and sheets of cold material that well down from the surface (see Figure 2). We can identify these hot blobs with mantle plumes believed to give rise to Earth’s hot spots, and the cold sheets with subducting oceanic plates.
Figure 2. Two hemispherical views of a global mantle convection model. Cold material (blue to grey colors) sinks towards the core-mantle boundary, while hot low-viscosity material (rainbow colors) rises towards the surface in focused upwellings
In this simple model, dynamics are easy to understand when considering the Rayleigh number \(Ra = \alpha \:\Delta Tg D^3/(\eta \: \kappa)\), where \(D\) is a length scale of the domain and \(\Delta T\) is the temperature difference between the surface and the core-mantle boundary; the larger the Rayleigh number, the smaller the features of the flow field. This already gives rise to a formidable computational challenge: with realistic values for the material parameters in the above equations, Earth’s flow features can be as small as a few kilometers across. With a volume of around 1012 km3, a finite element discretization using uniform meshes requires about a billion cells and a few 1010 unknowns to achieve appreciable accuracy. The current generation of codes can reach into this range, either through highly-tuned numerics or the use of adaptive mesh refinement. Some of the largest implicit finite element computations have indeed been performed on this problem [1, 2, 3, 5, 6, 7]. The most recent Gordon Bell prize was also awarded for the solution of this system on up to 500,000 cores [4].
A separate, equally-difficult challenge arises from the fact that realistic materials do not behave as outlined above. Rather than expand linearly with temperature, rocks undergo phase changes where density varies both continuously and discontinuously as a function of temperature and pressure. The same is true for viscosity, which may vary by many orders of magnitude even over small distances where hot and cold materials come together, such as when cold oceanic slabs subduct into the hot mantle. Viscosity also depends strongly and nonlinearly on stress, grain size, water content, and a number of other quantities. Finally, when considering the entire mantle, a mass conservation equation, \(\bigtriangledown \: \cdot \: (\rho u) = 0\), must replace the above incompressibility equation.
Cumulatively, compressibility and strong nonlinearities complicate the design of efficient and accurate solver strategies. The required nonlinear iterations also make the solution very expensive computationally. While the nonlinearity can be treated efficiently in a time-stepping scheme, solving the first time step self-consistently can become a significant challenge.
Equally difficult are the large jumps in viscosity that result from strong temperature gradients that often cannot be resolved adequately by the mesh. Such “essentially discontinuous” viscosity fields cause large discretization errors and pose enormous challenges to the design of linear solvers and preconditioners. Recent experiments show that appropriately averaging material parameters on each cell, without reducing the overall convergence order, can significantly reduce these discretization errors. This also vastly improves the efficiency of linear solvers, sometimes reducing the time to solution by a factor of ten on complex models.
Much more progress in all areas of computational mathematics—on discretizations, nonlinear and linear solvers, preconditioners, and parallel algorithms—is necessary to make the solution of complex models in mantle convection fast and routine. However, many of the most widely-used codes are open source and well-documented, such as Citcom, or our own contribution, ASPECT. Specifically, ASPECT is built as a modular platform that allows mathematical and computational scientists to test new discretizations and solvers on realistic problems, and geoscientists to develop and test new model descriptions. The Computational Infrastructure for Geodynamics has been collecting and curating these and other codes to facilitate both geodynamical research as well as experimentation with new numerical methods.
References
[1] Burstedde, C., Stadler, G., Alisic, L., Wilcox, L. C., Tan, E., Gurnis, M., & Ghattas, O. (2013). Large-scale adaptive mantle convection simulation. Geophysical Journal International, 192(3), 889-906.
[2] Gmeiner, B., Huber, M., John, L., Rüde, U., & Wohlmuth, B. (2015). A quantitative performance analysis for Stokes solvers at the extreme scale. Cornell University Library. Preprint, arXiv:1511.02134.
[3] Gmeiner, B., Rüde, U., Stengel, H., Waluga, C., & Wohlmuth, B. (2015). Performance and scalability of hierarchical hybrid multigrid solvers for Stokes systems. SIAM Journal on Scientific Computing, 37(2), C143-C168.
[4] Rudi, J., Malossi, A.C.I., Issac, T., Stadler, G., Gurnis, M., Staar, P.W.J.,…Ghattas, O. (2015). An extreme-scale implicit solver for complex PDEs: highly heterogeneous flow in earth’s mantle. Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis. Austin, TX: ACM.
[5] Stadler, G., Gurnis, M., Burstedde, C., Wilcox, L.C., Alisic, L., & Ghattas, O. (2010). The dynamics of plate tectonics and mantle flow: From local to global scales. Science, 329(5995), 1033-1038.
[6] Sundar, H., Biros, G., Burstedde, C., Rudi, J., Ghattas, O., & Stadler, G. (2012, November). Parallel geometric-algebraic multigrid on unstructured forests of octrees. Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis (p. 43). Salt Lake City, UT: IEEE Computer Society Press.
[7] Weismüller, J., Gmeiner, B., Ghelichkhan, S., Huber, M., John, L., Wohlmuth, B.,…Bunge, H.P. (2015). Fast asthenosphere motion in high-resolution global mantle flow models. Geophysical Research Letters, 42(18), 7429-7435.
Wolfgang Bangerth is a professor of mathematics at Texas A&M University. Juliane Dannberg is a graduate student at the GeoForschungsZentrum Potsdam in Germany and a research assistant at Texas A&M University. Rene Gassmöller is a postdoc at both the GeoForschungsZentrum Potsdam and Texas A&M University. Timo Heister is an assistant professor at Clemson University. All four are co-authors of the ASPECT code to simulate mantle convection. Bangerth and Heister are also among the principal authors of the deal.II finite element library.