Development of a Stokes flow simulation code for core formation
Mikito Furuichi 1 and Takashi Nakawaga 111

Institute for Research on Earth Evolution, Japan Agency for Marine-Earth Science and Technology (JAMSTEC),
3173-25 Showa-machi, Kanazawa, Yokohama, Japan

Introduction

Development of the simulation code for the Stokes flow
problems is an interesting challenge for the computational
geoscience, because the Stokes system is relevant to thenumerical study of long-time scale geodynamics of the planetaryinterior such as mantle convection [1-4]. Over thirty-five years,Stokes flow simulation code has been progressed to investigate
the global scale geodynamical evolutions with realistic settings,
for example, in the treatment of extremely viscosity variations,nonlinear rheology, transport of material, spherical shellgeometry and adaptive grid resolutions.
We have developed a new simulation code for solving theStokes flow motion named after “Nplat” to investigate the coreformation process of the planetary interior [2]. The formation ofa metallic core is widely accepted as the biggest differentiation
event during the final stage of the planetary formation [5]. In ourscenario, the dense metal-rich materials are formed by theimpact-induced events and sink to the center in the long timescales.
In order to deal the realistic mechanical setting of the coreformation process, our simulation employs free surface
condition under a self-gravitating field. We use the finite
difference discretization on a staggered grid. Expressing the freesurface motion, a stick air layer which is the low viscosity andzero density layer surrounding the planetary surface, is assumed.
An ill conditioned Stokes problem arsing in a sharp viscosity
contrast is solved by iterative Stokes flow solver, robust to large
viscosity jumps, using a strong Schur complement
preconditioner and mixed precision arithmetic utilizing thedouble-double method [3].
The core formation process could be triggered by impacts ofplanetesimals, but the impact induced process in the planetaryinterior had not been fully understood yet [5-7]. Here in order toinvestigate the various possible scenarios, our simulation cantreat the several simplified models of impact event.
In this presentation, we introduce the overview of oursimulation algorithms [1-4], and demonstrate the preliminary
simulation results on massively parallel supercomputer Earth
Simulator 2.
2

Problem settings

In this section, we introduce basic simulation settings of ourcore formation simulation. Our scenario starts with a Mars-sized
protoplanet (PP) assumed as an undifferentiated planet forsimplicity [6]. In the accretionary stage, the planet grows withmany times impacts. As a consequence of each impact, molten
regions such as magma ponds or oceans are formed on theplanetary surface. During solidification of such regions, metal–
silicate separation locally occurs over short time-scales. Thetarget of our simulation is to reveal the essential global
dynamics of the planetary interior after the solidification overlong time-scales [7].

2.1

Stokes flow problem of Core formation simulation
We assume the Stokes flow problem in simple box domain  (0, 0, 1)  (0, 1, 0)  (1, 0, 0) of the three dimensional
Cartesian geometry under the self-gravitational field. The freeslip boundary condition is applied to every wall. Three materialcomponents which are form the protoplanet (PP) (s = 1), metalrich (s = 2) and silicate-rich materials (s = 3), are assumed in oursimulation. The surface of the material is tracked as the scalarcolor functions 0 < cs < 1, which is discretized at the center ofthe cell. The variables on the grid cells for example
density  are obtained by,(1)  c  .
sss

We define the color function of sticky air as cair  1  s cs thathas zero mass density and very small viscosity [4, 7].
The equation of the Stokes flow problem under anincompressible approximation is described as the balance ofviscous stress, pressure gradient and gravitational force given as,  u j ui  p  fi .
x j  xi x j  xiu  0(2)where xi is position (i=1,2,3 for three dimensions), and  ij , p,and fi are the deviatoric stress, pressure and gravity force
respectively. Here we assume Newtonian rheology. The
transport of color function is also solved.
The self-gravitating field is obtained by solving the gravity
potential  which satisfies 2   . The gravity force issimply derived by fi  i  .
2.2

Impact models
Our simulation code can deal with three types of impact
events. First model is the three-layered model to investigate theoverturn process which may occur when a dense metal-richlayer surround the central PP after the solidification of global
magma ocean [5].
The second model is the small iron diapers randomly
distributed in the planetary body, which represents thesolidification of the local magma ponds caused by the manyimpacts within a short time scale [7].
In the third model, we incorporate planetary growth as aresult of impacts. During the simulation, the impactor comesfrom the random directions with a constant time interval andcreates hemispheric magma pond on the surface. When amagma pond is solidified, the metal-rich material is segregatedfrom silicate-rich material, which sinks at the base of the pond.
33.1

MethodAdvection
In our simulation, the material transport including the color
function of sticky air can be expressed as the advection equation.

Our code can use two types of advection method, particle-in-cell(PIC) and CIP-CSLR method [2, 8]. Although both methods
provide low diffusive solutions, the CIP-CLSR is generally lessexpensive in a computational cost, but has more difficulty in thetreatment of sub-grid scale profile and multiple components,
than the PIC method.
Stokes flow solver
Since the viscosity contrast at the material boundary is sharpand large, an algebraic equation for (3) has strongly illconditioned structure. In order to solve such an ill-conditioned
problem in 3-D, we have developed the robust iterative Stokes
flow solver [3, 4]. We use fully coupled approach for thevelocity-pressure saddle point system of (3) and (4) utilizing
Schur complement preconditioner dependent on the local
viscosity. Here the robust preconditioned Krylov subspacemethod can be a key for solving strongly ill-conditioned saddle
point problems [9]. We therefore choose Arnoldi type GCR
method as the Krylov subspace method combined owing to its
robustness compared to the alternative Lanzos method. We also
improve the robustness of the solution method by a mixed
precision technique for the inner problem of velocity. In themixed precision approach, the quad precision arithmetic isapplied throughout the GCR algorithm except for during theapplication of the preconditioner. This is consistent with the ideathat the preconditioner is used for a rough estimation of the realsolution.
The quad precision arithmetic implements into the doubledouble precision algorithms (DD) [3], which operates a pair ofdouble precision terms to express the quad precision. The loop
calculation in DD precision can be vectorized, but not in theofficially supported for the quad precision arithmetic. This is anadvantage on the DD for running the code on the ES2. Details of
this solver and numerical examples on the problems with locallyand highly varying viscosity are presented in the references [3].

Fig. 1. Snapshots of the steady state calculation. Outer blue semitransparent iso-surface and half cropped white isovolume represent thesilicate-rich and metal-rich materials respectively. Color on the slice
shows the distribution of color functions.

3.24Results

In this section, we briefly introduce three examples of coreformation simulation results. Fig. 1 shows the test calculation of
three layer model in which the initial cubic layer (left figure)deforms to the steady state of self-gravitation (i.e. sphere ofright figure). In this calculation we assume uniform viscosityand density for three layers. The global overturn scenario in thecore formation can be simulated by switching the viscosity anddensity to the target material parameters.
The result of iron diapir model is shown in Fig. 2. The
diapirs which are initially placed inside the planet at random
sink to the center. This simulation setting is useful to investigate
the behaviour of small diaper controlled by the rheologicalproperties of planetary interior.
The simulation of the scenario of the growing planet withmany impact events over long time scale is also presented in Fig.
3. This simulation can the change the mass balance of the planetwith increasing the total volume by the impact events. This is
the first 3-D core formation simulation code to incorporate
planetary growth as a result of impacts.
Our simulation is found to successfully capture the large
deformation of the materials following the three possible
scenario of core formation with maintaining a sharp interface
between them.

Fig. 2. Example of small metal diaper model calculation. Outer semitransparent glass, brown, and metric objects represent silicate-rich, PP,and metal-rich material layer surfaces. The material propertiesare  =1020[Pa.s],   3,900 [Kg/m3],  =1020[Pa.s],   3,500
11

2

2[Kg/m3] and  =1019[Pa.s],   7,000 [Kg/m3]
33Fig. 3 Snapshot of the core formation caused by a number of impacts.
The semi-transparent blue surface represents the silicate-rich material.
The white and orange object shows metal-rich and PP material
respectively.

5Acknowledgements

All calculations presented here are performed on Earth
Simulator 2 at the Japan Agency for Marine-Earth Science andTechnology. The present research was made possible by Grantin-Aid for Young Scientist (24740310) by Japan Society for thePromotion of Science (JSPS).
References
Tackley, P. J., Xie, S., Nakagawa, T., Hernlund, J. W., “Numericaland laboratory studies of mantle convection: Philosophy,accomplishments and thermo-chemical structure and evolution”,
AGU monography, Vol. 160, pp.83-99, (2005).
[2] Furuichi, M., Kameyama, M., Kageyama, A., " Three-Dimensional
Eulerian method for large deformation of viscoelastic fluid:
Toward plate-mantle simulation ", J. Compt. Phys, Vol 227,
pp.4977-4997, (2008).
[3] Furuichi, M., May, D., Tackley. P.J., " Development of a Stokes
Flow Solver Robust to Large Viscosity Jump Using a Schur
Complement Approach with Mixed Precision Arithmetic ", J.
Compt. Phys, 230, 8835-8851, (2011)[4] Furuichi, M., " Numerical modeling of three dimensional selfgravitating Stokes flow problem with free surface", Procedia.
Computer. Science. Vol 4, pp.1506-1515, (2011).
[5] Stevenson, D. J., (1990), “Fluid dynamics of core formation”, In“Origin of the Earth”, Cambridge Univ. Press, pp.231-249.
[6] Senshu, H., Kuramoto, K., Matsui, T., “Thermal evolution of agrowing Mars”, J. Geophys. Res., Vol. 107, (2002).
[7] Golabek, G. J., Gerya, T. V., Kaus, B. J. P., Ziethe, R., Tackley, P.
J., " Rheological controls on the terrestrial core formation
mechanism", Geochem. Geophys. Geosyst., Vol 10, Q11007,(2009).
[8] Xiao, F., Yabe, T., Peng, X., Kobayashi, H., (2002). “Conservativeand oscillation-less atmospheric transport schemes based onrational functions”, J. Geophys. Res. (D22), pp.4609.
[9] May, D.A., Moresi, L., " Preconditioned iterative methods forStokes flow problems arising in computational geodynamics",
Phys. Earth Planet Interiors. Vol 171, pp.33–47, (2008).
[10] Bailey, D. H., " High-Precision Software Directory",
http://crd.lbl.gov/~dhbailey/mpdist/[1]