Acknowledgements Introduction 1. The continuity equation 2. Density and gravity 3. Numerical solutions of partial differential equations 4. Stress and strain 5. The momentum equation 6. Viscous rheology of rocks 7. Numerical solutions of the momentum and continuity equations 8. The advection equation and marker-in-cell method 9. The heat conservation equation 10. Numerical solution of the heat conservation equation 11. 2-D thermomechanical code structure 12. Elasticity and plasticity 13. 2-D implementation of visco-elastic-plastic rheology 14. The multi-grid method 15. Programming of 3-D problems 16. Numerical benchmarks 17. Design of 2-D numerical geodynamic models Epilogue: outlook Appendix: MATLAB program examples References Index.

/pdf/introduction-to-numerical-geodynamic-modelling-27sy9p5rm3.pdf

Introduction to Numerical Geodynamic Modelling

The area of a spherical region can be easily measured by considering which sampling points of a lattice are located inside or outside the region. This point-counting technique is frequently used for measuring the Earth coverage of satellite constellations, employing a latitude-longitude lattice. This paper analyzes the numerical errors of such measurements, and shows that they could be greatly reduced if the Fibonacci lattice were used instead. The latter is a mathematical idealization of natural patterns with optimal packing, where the area represented by each point is almost identical. Using the Fibonacci lattice would reduce the root mean squared error by at least 40%. If, as is commonly the case, around a million lattice points are used, the maximum error would be an order of magnitude smaller.

Measurement of areas on a sphere using Fibonacci and latitude-longitude lattices

[1] The Moon possesses a clear dichotomy in geological processes between the nearside and farside hemispheres. The most pronounced expressions of this dichotomy are the strong concentration of radioactive heat sources on the nearside in a region known as the Procellarum KREEP Terrane (PKT) and the mare basaltic lava flows that erupted in or adjacent to this terrane. We model the thermochemical evolution of the Moon using a 3-D spherical thermochemical convection code in order to assess the consequences of a layer enriched in heat sources below the PKT on the Moon's global evolution. We find that in addition to localizing most of the melt production on the nearside, such an enriched concentration of heat sources in the PKT crust has an influence down to the core-mantle boundary and leaves a present-day temperature anomaly within the nearside mantle. Moderate gravitational and topographic anomalies that are predicted in the PKT, but not observed, may be masked either by crustal thinning or gravitational anomalies from dense material in the underlying mantle. Our models also predict crystallization of an inner core for sulfur concentrations less than 6 wt %.

https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1002/jgre.20103

Asymmetric thermal evolution of the Moon

Computations have helped elucidate the dynamics of Earth's mantle for several decades already. The numerical methods that underlie these simulations have greatly evolved within this time span, and today include dynamically changing and adaptively refined meshes, sophisticated and efficient solvers, and parallelization to large clusters of computers. At the same time, many of these methods -- discussed in detail in a previous paper in this series -- were developed and tested primarily using model problems that lack many of the complexities that are common to the realistic models our community wants to solve today. 
With several years of experience solving complex and realistic models, we here revisit some of the algorithm designs of the earlier paper and discuss the incorporation of more complex physics. In particular, we re-consider time stepping and mesh refinement algorithms, evaluate approaches to incorporate compressibility, and discuss dealing with strongly varying material coefficients, latent heat, and how to track chemical compositions and heterogeneities. Taken together and implemented in a high-performance, massively parallel code, the techniques discussed in this paper then allow for high resolution, 3d, compressible, global mantle convection simulations with phase transitions, strongly temperature dependent viscosity and realistic material properties based on mineral physics data.

High accuracy mantle convection simulation through modern numerical methods – II: realistic models and problems

https://www.brown.edu/research/projects/scientific-computing/sites/brown.edu.research.projects.scientific-computing/files/uploads/High%20order%20WENO%20and%20DG%20methods%20for%20time-dependent%20convection-dominated%20PDEs.pdf

High order WENO and DG methods for time-dependent convection-dominated PDEs

Numerical simulations of thermal convection in the Earth’s mantle often employ a pseudoplastic rheology in order to mimic the plate-like behavior of the lithosphere. Yet the benchmark tests available in the literature are largely based on simple linear rheologies in which the viscosity is either assumed to be constant or weakly dependent on temperature. Here we present a suite of simple tests based on nonlinear rheologies featuring temperature, pressure, and strain rate-dependent viscosity. Eleven different codes
based on the finite volume, finite element, or spectral methods have been used to run five benchmark cases leading to stagnant lid, mobile lid, and periodic convection in a 2-D square box. For two of these cases, we also show resolution tests from all contributing codes. In addition, we present a bifurcation analysis, describing the transition from a mobile lid regime to a periodic regime, and from a periodic regime to a stagnant lid regime, as a function of the yield stress. At a resolution of around 100 cells or elements in both vertical
and horizontal directions, all codes reproduce the required diagnostic quantities with a discrepancy of at
most $3% in the presence of both linear and nonlinear rheologies. Furthermore, they consistently predict
the critical value of the yield stress at which the transition between different regimes occurs. As the most
recent mantle convection codes can handle a number of different geometries within a single solution
framework, this benchmark will also prove useful when validating viscoplastic thermal convection simula-
tions in such geometries.

/pdf/a-community-benchmark-for-viscoplastic-thermal-convection-in-3qqps505rb.pdf

A community benchmark for viscoplastic thermal convection in a 2-D square box

[1] This paper presents a new method to generate a three-dimensional spherical grid using natural neighbor Voronoi cells distributed by spiral functions. A unique property of this grid is the complete removal of symmetries with arbitrary selectable lateral and radial resolution, which are not restricted to discrete radial levels or geometrical constraints as compared to the commonly used grids based on projected triangulated platonic solids such as a cube, a rhomboid, or an icosahedron. The spiral grid can be refined in certain areas of interest and makes it possible to have a very small inner radius to outer radius ratio. Cell volumes can be made almost constant throughout the computational domain. Analysis, statistics, and computation methods are described in detail, as well as a possible domain decomposition suitable for parallel computing. Conductive temperature profiles were numerically calculated in the spherical shell and directly compared with the analytic solution as verification. The grid is applied to numerical simulations of mantle convection using a finite volume scheme. The model is validated by a comparison of steady state cubic and tetrahedral convection patterns with other published models.

https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2007GC001581

The spiral grid: A new approach to discretize the sphere and its application to mantle convection

https://hal.archives-ouvertes.fr/hal-00532161/document

Finite volume discretization for dynamic viscosities on Voronoi grids

Mantle convection in terrestrial planets is strongly influenced by the temperature dependence of its viscosity. Considering this property, the mantle convection is basically divided into three different regimes depending on the viscosity contrast: stagnant lid, sluggish and mobile lid regime [1]. In the first regime the mantle is divided into two parts an active part which is convecting and another one on top which is immobile. In the last two regimes the surface material can move and is incorporated into the mantle convection. The temperature dependence of the mantle viscosity can be modeled using the so called Arrhenius Law. In the Arrhenius formulation the temperature dependence of the viscosity for a silicate mantle is given by the activation energy. Using realistic values for the activation energy in the Arrhenius formulation will result in large viscosity contrasts (∼ 10) which cannot be handled very well by the numeric. Therefore, an approximation of the viscosity is commonly used to model the mantle convection, i.e., the Frank-Kamenetskii approximation. This approximation linearises the Arrhenius law, suggesting a viscosity which is many orders of magnitude smaller at the surface. Although the approximation has been shown to represent only the stagnant lid regime correctly, it is widely used in the literature also for the other regimes. We present a comparison of the mantle convection for stagnant lid cases using either the Arrhenius law or the Frank-Kamenetskii approximation with a 3D spherical code, GAIA [2,3].The results confirm earlier studies that the Frank-Kamenetskii approximation can be used in the stagnant lid regime for a fixed Rayleigh number and a sufficiently thick stagnant lid. However, several problems arise when using this approximation in numerical simulations for a cooling mantle or for convection with a thin stagnant lid. A systematic study is presented highlighting the differences in various control parameters of mantle convection, e.g. degree of convection, Nusselt number, and stagnant lid thickness.

/pdf/mantle-convection-in-a-spherical-shell-the-problems-of-using-1hqov37qkb.pdf

C. Hüttig

Papers

A community benchmark for viscoplastic thermal convection in a 2-D square box

The spiral grid: A new approach to discretize the sphere and its application to mantle convection

Finite volume discretization for dynamic viscosities on Voronoi grids

Mantle Convection in a Spherical Shell: The problems of using Frank-Kamenetskii Approximation for the viscosity law

Thermal evolution and heat-pipe melt transport: implications for one plate planets