Unravelling the interior Evolution of Rocky Planets Through Large-Scale Numerical Simulations

Principal Investigator:
Ana-Catalina Plesa

Institute of Planetary Research, Planetary Physics, German Aerospace Center/DLR, Berlin (Germany)

Local Project ID:

HPC Platform used:
Hornet of HLRS

Date published:

The large amount of data returned by several space missions to the terrestrial planets has greatly improved our understanding of the similarities and differences between the innermost planets of our Solar System. Nevertheless, their interior remains poorly known since most of the data is related to surface processes. In the absence of direct data of the interior evolution of terrestrial planets, numerical simulations of mantle convection are an important mean to reconstruct the thermal and chemical history of the interior of the Earth, Moon, Mercury, Venus and Mars. In this project, run on Hornet of HLRS, researchers used the mantle convection code Gaia to model the thermal evolution of terrestrial planets and in particular the early stage of their history. 

Rocky planets like the Earth, Venus, Mercury, and Mars have been built during the early stages of our Solar System through collisions of planetesimals and protoplanets. Such energetic collisions but also the decay of short-lived radiogenic elements like 26Al and 60Fe together with the energy released by core formation produced large amounts of heat that led to a global magma ocean. Upon planetary cooling, the molten material will crystallize, but the exact solidification sequence is a complex process that is only poorly understood. Nevertheless, this process is one of the most important episodes in a planet’s history as it sets the stage for the later evolution and may be responsible for the subsequent style of heat transport.

All terrestrial planets are essentially heat engines that cool over time by losing their interior heat through their surface. The slow creep of silicate material, a process also known as mantle convection, is ultimately responsible for the heat transport from the deep interior towards the surface and may be linked to surface expressions like volcanic provinces or tectonic plates.

On Earth, convection cells reach up to its surface, which is divided into seven major and other smaller tectonic plates that are permanently recycled back into the mantle at subduction zones (i.e., oceanic trenches) and created anew at spreading centers (i.e., mid-ocean ridges). On all other terrestrial bodies in our Solar System, instead, convection develops underneath an immobile layer, the so-called stagnant lid through which heat is transported solely by conduction. Why and when plate tectonics started to operate on Earth but never developed on other terrestrial planets like Mars, Mercury Venus or the Moon is still poorly known. Although an early and short episode of plate tectonics has been proposed for Mars and the surface of Venus is geologically very young (300 – 600 Myr), no indisputable evidence exists that plate tectonics has ever been active on such bodies. Moreover, the early history of the Earth was characterized most likely by crustal delamination, a type of surface mobilization that may be responsible also for the young crust observed on Venus, but that is different than plate tectonics observed today on Earth.

Because their crust has not been constantly renewed, old surfaces of terrestrial planets like Mars, Mercury or the Moon have been preserved and are an important piece of evidence of the earliest processes that may have been active also on the Earth.

The large amount of data returned by several space missions to the terrestrial planets, has greatly improved our understanding of the similarities and differences between the innermost planets of our Solar System. Nevertheless, their interior remains poorly known since most of the data is related to surface processes. Samples, which have been either collected during the Apollo missions to the Moon or delivered to Earth in form of meteorites from the Moon and Mars, are too few to unambiguously reconstruct the history and present-day state of the interior. Moreover, in-situ seismic and heat flow measurements, which can be seen as the most direct measurements of the interior since they can be used to constrain the interior structure and the amount of heat available, are solely available for the Earth and Moon but will be also performed on Mars by the upcoming InSight mission.

In the absence of direct data of the interior evolution of terrestrial planets, numerical simulations of mantle convection are an important mean to reconstruct the thermal and chemical history of the interior of the Earth, Moon, Mercury, Venus and Mars. Large-scale models that can accurately track the thermal anomalies and compositional heterogeneities over billions of years of evolution are typically run with thousands of computational cores in parallel on machines available at high-performance computing (HPC) centres like the Hornet System at HLRS in Stuttgart. Such models solve the conservation equations of mass, momentum and energy in various geometries and can tackle a complex rheology as inferred by deformation experiments performed under mantle temperature and pressure conditions.

In project MAntle THErmo-chemical COnvection Simulations (MATHECO), researchers of the Institute of Planetary Research of the German Aerospace Center, Berlin have used the mantle convection code Gaia to model the thermal evolution of terrestrial planets and in particular the early stage of their history. At present only a handful of mantle convection codes world-wide are able to tackle complex thermal evolution simulations in a 3D spherical geometry and track chemical heterogeneities in the interior of terrestrial planets.

Gaia is a fluid flow solver written in C++ that does not require any additional libraries and is thus easily portable on various HPC systems. Over the duration of the project the code has been improved and is now able to use arbitrary geometries (spherical, Cartesian or cylindrical) as long as the meshes are Voronoi meshes. Moreover, the code brings its own matrix solvers and, over the past years, various types of solvers have been included. The users can choose between TFQMR, GMRES and BiCGStab, with the latter being the best choice in most of the tested cases.

To accurately track a large number of chemical fields, Gaia employs the particle in cell method (PIC). This method uses massless tracer particles that move according to the velocity field computed on the mesh. The advantage of PIC method over classical grid based methods is that it can track the motion of a large number of chemical species by solving only one transport equation and is essentially free of numerical diffusion.

In addition, Gaia has its own routines to decompose the mesh into equal portions that are distributed on computational cores. Fig. 2a shows examples of domain decompositions for various meshes. When running the code in parallel, the so-called halo-cells are a first measure for the performance. Halo-cells or ghost-cells are computational cells used to exchange information between various computational domains. An optimal performance is achieved when their number is kept minimal. In the course of the MATHECO project, the code performance has been tested with up to 54,000 computational cores for a large 3D Cartesian grid with 275 million degrees of freedom (DOF) and with 12,000 computational cores for a 3D spherical grid with 60 million DOF and 240 million particles (Fig. 2b). A detailed description of the mantle convection code Gaia can be found in the following publications (Hüttig and Stemmer 2008a; Hüttig and Stemmer 2008b; Hüttig et al., 2013; Plesa et al., 2016).

Gaia writes periodically output and snapshot files. While output files contain data that is further post-processed to obtain the results of a specific model, snapshot files store additional information necessary to restart a simulation. During the accounting period of the MATHECO project, the output and snapshot writing routines have been improved. Gaia now uses the MPI-I/O facilities that in addition to a speed-up of a factor of 100 over the previously used method also reduce the number of files needed to be written during a specific simulation. The I/O performance of the recently implemented routines using the MPI-I/O facilities is shown in Fig. 2c. Further results obtained using the MPI-I/O routines are summarized in a Bachelor Thesiswritten at the University of Applied Sciences (HTW) and German Aerospace Center (Willich, 2016).

In project MATHECO, researchers of the Institute of Planetary Research of the German Aerospace Center, Berlin have employed the code Gaia to investigate the formation and stability of early built geochemical reservoirs in the interior of Mars, as suggested by the geochemical analysis of Martian meteorites. In addition, the focus of their research was to understand the consequences that such chemical reservoirs have on the long-term evolution of the Martian mantle, since chemical heterogeneities in planetary mantles can strongly influence the convection structure, hence the heat transport and with it the overall planetary evolution.

A variety of different scenarios involving a magma ocean have been proposed to explain chemical heterogeneities in planetary mantles. It is suggested that early in the evolution of a planet, the large amount of primordial heat due to accretion and core formation can give rise to a magma ocean as a consequence of significant or perhaps even complete melting of the mantle. Fractional crystallization of the mantle leads to a gravitationally unstable layering that may cause a mantle overturn. In this scenario, chemical convection plays a crucial role and can even suppress the onset of thermal convection. However, up to now most of the studies that investigated the consequences of a magma ocean crystallization and subsequent overturn have made strong assumptions on the domain geometry and mantle rheology.

Assuming scaling parameters appropriate for Mars, systems heated either solely from below or from within have been investigated by varying systematically the buoyancy ratio B, which measures the relative importance of chemical to thermal buoyancy, and the mantle rheology, by considering systems with constant, strongly temperature-dependent and plastic viscosity. A large set of simulations, spanning a wide parameter space, have been employed in order to understand the basic physics governing the magma ocean cumulate overturn and its consequence on mantle dynamics. Moreover, scaling laws have been derived that relate the time over which chemical heterogeneities can be preserved (mixing time) and the critical yield stress (maximal yield stress that allows the lithosphere to undergo brittle failure) to the buoyancy ratio. The results show that the mixing time increases exponentially with B, while the critical yield stress shows a linear dependence. This study has been published in the Journal of Geophysical Research, Planets, (Tosi et al., 2013).

In a further study, the researchers investigated Mars' early thermo-chemical evolution assuming a detailed magma ocean crystallization sequence available in the literature. An immobile lid (the stagnant lid) forms rapidly because of the strong temperature dependence of the viscosity and traps the densest cumulates and radioactive heat sources, which have been enriched during the freezing-phase of the magma ocean in the uppermost layers of the planet. The formation of the stagnant lid prevents the uppermost dense cumulates to sink, even when allowing for a plastic yielding mechanism. Below this dense stagnant lid, the mantle chemical gradient settles to a stable configuration. The convection pattern is dominated by small-scale structures which are difficult to reconcile with the large-scale volcanic features observed over Mars' surface. Assuming that the stagnant lid will break, a stable density gradient is obtained, with the densest material and the entire amount of heat sources lying above the core-mantle-boundary. This configuration leads to a strong overheating of the lowermost mantle, whose temperature increases and leads to the formation of a basal magma ocean, while the upper mantle, in the absence of heat sources, cools to a nearly conductive state. The results of this study have been published in the Earth and Planetary Science Letters (Plesa et al., 2014).

The above-mentioned studies suggest that chemical gradients obtained from the fractional crystallization of a global and deep magma ocean can entirely suppress thermal convection and lead to a conductive mantle. This is in particular difficult to reconcile with the volcanic history of Mars which indicates that the planet has been volcanically active up to the recent past, this being at best explained by an active interior.

An alternative scenario assumes that the magma ocean crystallized homogeneously without large chemical gradients established upon solidification. In fact, a recent study has shown that a sufficiently high Rayleigh number, indicating a vigorously convecting interior, guarantees the onset of solid-state convection prior to complete crystallization of the mantle (Maurice et al., 2017). Thus, convective mixing may reduce density gradients prior to complete solidification and in such case, the density contrast at the end of the crystallization sequence may be strongly diminished. Density variations may be produced shortly after the mantle solidified, by partial melting of the mantle being caused by large scale impacts and/or heating by the decay of radioactive isotopes of potassium, thorium and uranium. In fact, density variations of only 40 – 60 kg/m3 are sufficient to maintain chemical heterogeneities in an actively convecting mantle in particular if the chemical reservoirs are associated with viscosity variations (the reservoirs have a higher viscosity than the surrounding mantle material). In addition, if such reservoirs become trapped in the stagnant lid, they can be preserved over geological time scales. The results of such alternative scenarios have been discussed and published in a recent review study in Meteoritics and Planetary Science (Breuer et al., 2016).

In the MATHECO project, researchers of the Institute of Planetary Research of the German Aerospace Center in Berlin have shown that the mantle convection code Gaia can be efficiently employed to run large-scale numerical models on more than 50 thousand computational cores in parallel. Such large-scale numerical simulations can only be performed on dedicated high-performance computing centers like the HLRS supercomputing center in Stuttgart. Gaia has become one of the few codes world-wide able to efficiently run large-scale simulations of interior evolution and, thanks to the available computing resources of HLRS, high resolution mantle convection models can be employed to better understand the physical processes active in the interior of terrestrial bodies using realistic parameter values.

Project relevant publications:

Breuer D., Plesa A.-C., Tosi N., and Grott M. (2016): Water in the Martian interior — The geodynamical perspective. Meteoritics and Planetary Science 1–34; doi: 10.1111/maps.12727.

Hüttig C. and Stemmer K. (2008b): Finite volume discretization for dynamic viscosities on Voronoi grids. Physics of the Earth and Planetary Interiors; doi: 10.1016/j.pepi.2008.07.007.

Hüttig C. and Stemmer K. (2008a): The spiral grid: A new approach to discretize the sphere and its application to mantle convection. Geochemistry Geophysics Geosystems, Q02018, 9; doi: 10.1029/2007GC001581.

Hüttig C., Tosi N., and Moore W. B. (2013): An improved formulation of the incompressible Navier-Stokes equations with variable viscosity. Physics of the Earth and Planetary Interiors, 11–18, 220; doi: 10.1016/j.pepi.2013.04.002

Maurice, M., Tosi N., Samuel H., Plesa A.-C., Hüttig C., and Breuer D. (2017): Onset of solid-state mantle convection and mixing during magma ocean solidification. Journal of Geophysical Research: Planets, 122, 577–598; doi: 10.1002/2016JE005250.

Plesa A.-C., Tosi N., and Breuer D. (2014): Can a fractionally crystallized magma ocean explainn the thermo-chemical evolution of Mars? Earth and Planetary Science Letters403:225–235; doi: 10.1016/j.epsl.2014.06.034.

Plesa A.-C., Hüttig C., Maurice M., Breuer D., and Tosi N. (2016): Large Scale Numerical Simulations of Planetary Interiors. High Performance Computing in Science and Engineering ’15 ed. Wolfgang E. Nagel, Dietmar H. Kröner, Michael M. Resch; doi: 10.1007/978-3-319-24633-8_43

Tosi N., Plesa A.-C., and Breuer D. (2013): Overturn and evolution of a crystallized magma ocean: A numerical parameter study for Mars. Journal of Geophysical Research: Planets118:1512–1528; doi: 10.1002/jgre.20109.

Willich, F. (2016): MPI I/O compared to POSIX I/O for Numerical Simulations on Supercomputers: Differences in Performance and Operation Principle (Bachelor Thesis), University of Applied Sciences (HTW), Berlin; German Aerospace Center (DLR), Berlin.

Scientific Contact:

Ana-Catalina Plesa
German Aerospace Center
Institute of Planetary Research, Planetary Physics
Rutherfordstraße 2, D-12489 Berlin (Germany)
e-mail: ana.plesa[at]

Tags: DLR HLRS Astrophysics