ENVIRONMENT AND ENERGY

Introduction

The idea of a thin, weak layer beneath the Earth’s rigid outer shell is over a century old, long predating the widespread acceptance of plate tectonics. This layer, termed the asthenosphere, was originally proposed in order to account for Earth’s internal support of topographic loads such as mountain ranges (isostasy). Our view of this layer expanded following the plate tectonic revolution of the late 1960’s, after which it came to be viewed as a lubricating layer to facilitate plate motion; passively flowing in response to the overlying plates. In the past two decades, it has become clear that the asthenosphere is a much more exciting and active part of the convecting mantle than previously imagined. Upwelling mantle plumes from deep within the Earth (found, for instance, beneath Hawaii) feed this layer with hot material. This creates high pressure regions within the asthenosphere, while sinking tectonic plates produce low pressure regions. The resulting pressure gradients can induce strong pressure-driven flow, in addition to the component of flow driven by the overlying plates. This fast asthenospheric flow can allow for the lateral trans-port of plume-sourced material over significant distances in relatively short geological timescales. As a result, surface uplift and volcanism may be induced in regions where we would not have previously expected it, which could even have implications for past and future changes in sea level. In addition, the flow can actually alter the motion of the overlying tectonic plates if pressure-driven flow velocities become high enough. The asthenosphere is therefore a crucial layer in the Earth system; providing a link between deep mantle dynamics and surface processes.

While the importance of the asthenosphere is now well understood, simulations which accurately include such a layer remain a computational challenge. Fluid dynamical simulations of past mantle convection have now become a ubiquitous tool in the field of global geodynamics. However, resolving the sharp viscosity contrasts and short flow wavelengths associated with the asthenosphere require a grid-point spacing of ~10 km at the very least. In a global model, this equates to ~670 million grid points. Solving the conservation equations for the velocity, temperature, and pressure at each grid point means solving a system of equations with several billion unknowns during each time step. A simulation covering a significant time span of geological history (say, 400 million years) may require ~10,000 such time steps. With these numbers in mind, the necessity of modern HPC systems in solving these problems becomes abundantly clear. Few geodynamics groups around the world have access to vast computational resources such as those found at the LRZ. Therefore, as it stands, the vast majority of published models remain at a grid spacing of ~25 km or coarser. It is due to this that influence of a thin, low viscosity channel in global models of past mantle convection has so far remained elusive.

Results and Methods

During this project, we have generated high-resolution simulations of past mantle convection which feature a low viscosity asthenospheric channel in the uppermost mantle. Simulations were carried out using the parallel finite element mantle convection code TERRA. Tectonic plate histories were assimilated at the surface of the simulation as boundary conditions on surface velocity, allowing us to study the resulting mantle flow. An Earth-like initial temperature field was first generated by simulating 1 billion years of convection with free-slip velocity boundary conditions at the surface. Following this, tectonic velocities were assimilated from 400 Ma (million years ago) to the present day. The generation of the initial state alone requires ~300,000 core hours on the SuperMUC-NG, while the remainder of each simulation requires an additional ~250,000 core hours. Even when parallelised amongst 2,048 cores, this is equivalent to 11 days of compute time for a single simulation. The simulation has a resolution of ~10 km, allowing us to create a thin (~100 km), low viscosity (~ 5^.1019 Pa.s) channel beneath a stiff outer shell, the tectonosphere. In addition to radial viscosity layering, viscosity was made temperature-dependent, meaning that cold, sinking tectonic plates can retain their strength while passing through the asthenosphere.

Figure 1 shows four transects of the African hemisphere, spanning from 60 Ma (shortly after the extinction of the dinosaurs at 66 Ma) to the present day (0 Ma). A contour of temperature deviation has been extracted at 400 Kelvin in order to highlight hot, upwelling mantle plumes. Old, sinking tectonic plates can be seen as cold (blue) regions in the background. For comparison, we also ran simulations at a more “standard” resolution of ~ 25 km grid-point spacing. Models such as this are far more limited in their ability to simulate asthenospheric dynamics, and instead can only feature a thicker (~ 200 km) and higher viscosity (~ 5^.1020 Pa.s) asthenosphere. A comparison between the two is shown in Figure 2. Temperature deviations within the asthenosphere of each model at 200 km depth are shown at the same time steps as Figure 1. The difference in the character of asthenospheric temperature anomalies is immediately clear. Note the long-wavelength character of the hot regions in the standard resolution model (left), as plume material cannot spread out effectively when it reaches the uppermost mantle. Meanwhile, the lower viscosity channel allowed by the high resolution model (right) allows for the fast transfer of plume material in response to pressure gradients and overlying plate velocities. Hot material is carried towards subduction zones (blue lines) where it builds up, an effect which is not reproduced in the standard resolution model. The differences between the asthenospheric flow patterns between these two models has broader implications for their prediction of surface uplift (dynamic topography) and volcanism. As such, high resolution models as the one shown in Figure 1 ought to become the new “standard” – something which can only become possible with computational resources such as the SUPERMUC-NG.

Ongoing Research / Outlook

While the results presented here have provided invaluable insight into asthenospheric dynamics, the grid spacing of ~10 km is still not fine enough to simulate viscosity structure of the real asthenosphere inferred from observations. For this, a grid point spacing of ~1–5 km would be required globally, leading to sparse linear systems with more than a trillion unknowns. This necessitates scalable codes which can take advantage of exascale-capable modern HPC systems. To this end, our TerraNeo project has focused on the development of a new mantle convection code, capable of extreme-scale simulations. The project is a collaboration between FAU Erlangen-Nürnberg, LMU Munich, TU Munich and the Leibniz Supercomputing Center (LRZ). The TerraNeo application is built upon the HyTeG [1] (Hybrid Tetrahedral Grids) framework, an open-source C++ framework aimed at efficient and scalable multigrid solvers for Finite Element simulations. The extreme scalability and performance of HyTeG’s matrix-free multigrid methods has already been demonstrated on the SuperMUC-NG. The solver was able to solve the Stokes equations with more than a trillion (10¹²) unknowns in less than a minute when parallelised amongst 147,456 processes [2]. The TerraNeo application has now been equipped with the essential physical features, such as tectonic plate velocities and temperature-dependent viscosity, in order to generate the first ever simulations of global mantle circulation with ~1 km resolution within the coming year.

References and Links

[1] i10git.cs.fau.de/hyteg/hyteg

[2] N. Kohl and U. Rüde, “Textbook efficiency: massively parallel matrix-free multigrid for the Stokes system”, SIAM Journal on Scientific Computing, Vol. 44, no. 2, C124–C155, 2022