PaperlayerView exact revision on arXiv

Source-linked AI summary

High Accuracy Mantle Convection Simulation through Modern Numerical Methods. II: Realistic Models and Problems

Timo Heister, Juliane Dannberg, Rene Gassmöller, Wolfgang Bangerth

arXiv:1702.05075v2cs.CE

TL;DR

Mantle-convection methods were often tested on simplified problems, while realistic simulations require compressibility, discontinuities, phase changes, and complex geometries. This paper revisits numerical algorithms using experience with realistic Aspect applications and reports techniques that jointly support accurate large-scale simulations. The authors also identify limitations in automatic mesh adaptation and pressure approximation.

  • Problem

    Existing numerical methods were commonly developed and tested on simplified problems, whereas realistic mantle models require compressibility, discontinuous coefficients, complex geometries, latent heat, and chemical transport.

  • Method

    The paper revisits time stepping, mesh refinement, compressibility, discontinuous coefficients, latent heat, and additional-quantity advection using experience with complex Aspect applications.

  • Results

    Excellent agreement was obtained for Nusselt number, root mean square velocity, average temperature, total shear heating, and adiabatic heating against previously reported compressible-convection results.

  • Takeaways & Limitations

    Taken together with earlier methods, the discussed techniques help solve complex mantle-convection models on large-scale compute clusters.

  • Takeaways & Limitations

    Automatic mesh adaptation that accounts for the quantity of interest remains unclear, leading to a flexible two-tier system driven primarily by user input.

Abstract

from arXiv · show

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.

1 Introduction

The paper addresses the gap between numerical methods tested on simplified mantle-convection problems and the realistic, physically complex models researchers now seek to solve. It revisits algorithmic choices using experience with complex applications and identifies methods suited to accurate, efficient computation.

  • 1 Introduction: Existing mantle-convection methods have progressed from simple uniform 2d discretizations to adaptive 3d meshes, higher-order elements, and complex solvers.The paper builds on a previous article describing state-of-the-art methods.
  • 1 Introduction: Realistic applications require compressibility, discontinuous coefficients, complex and time-varying geometries, latent heat, and transport of chemical inhomogeneities.These requirements differ substantially from the simplified model problems commonly used to test numerical methods.
  • 1 Introduction: The paper revisits time stepping, mesh refinement, compressible formulations, and other computational-code choices for complex mantle-convection problems.The authors consider when earlier methods need adaptation or replacement.
  • 1 Introduction: The discussion is based on five years of experience applying the Aspect code to increasingly complex and realistic mantle-convection problems.The authors emphasize that the lessons concern general numerical methods, not only Aspect.
  • 1 Introduction: The paper also identifies open methodological and mathematical problems relevant to numerical-methods and numerical-analysis researchers.The intended audience includes both code developers and readers seeking to understand modern mantle-convection codes.
  • 1 Introduction: The paper covers time stepping, compressibility, discontinuous coefficients, latent heat, mesh refinement, and advection of additional quantities.These topics organize the treatment of realistic mantle-convection computations.

2 Formulation of the problem

The paper formulates mantle convection as compressible, anelastic flow with thermal and compositional buoyancy and several energy-balance effects. Material properties may vary with state and location, including discontinuous changes associated with phase transitions.

  • 2 Formulation of the problem: The model treats mantle convection as compressible, anelastic flow driven by thermal or compositional buoyancy.It includes friction, adiabatic heating, radiogenic heat production, and latent heat in the energy balance.
  • 2 Formulation of the problem: Inertial and elastic effects are neglected because the modeled velocities are very low and the time scales are long.This is an explicit modeling assumption rather than a consequence of the numerical method.
  • 2 Formulation of the problem: Aspect is an open-source project providing a modern, parallel, extensible code for mantle-convection simulation.The code is supported by the Computational Infrastructure for Geodynamics initiative and the National Science Foundation.
  • 2 Formulation of the problem: The formulation uses velocity, pressure, and temperature together with stress and thermodynamic material properties.The displayed equations introduce the primary fields and constitutive quantities used by the model.
  • 2 Formulation of the problem: Viscosity, density, diffusivity, heat production, expansion, and gravity may depend on temperature, pressure, strain rate, and location.Viscosity depends additionally on strain rate, while gravity depends on location.
  • 2 Formulation of the problem: The model permits coefficients to jump discontinuously, as in thermodynamically consistent models with phase changes.This discontinuity is a central distinction from many simplified benchmark problems.
  • 2 Formulation of the problem: ALA, TALA, and Boussinesq approximations can be derived by assuming density variations are small relative to hydrostatic density increase.The paper discusses how these approximations relate to the broader formulation.

3 Numerical methods

The numerical-methods section examines the modifications needed when moving from simplified mantle-convection benchmarks to realistic convective transport. It uses finite elements, stabilized advection, iterative solvers, and open-source high-performance software as the computational framework.

  • 3 Numerical methods: The comparison framework uses Taylor-Hood finite elements for Stokes flow and block-preconditioned GMRES for the resulting linear systems.The temperature equation is also discretized with finite elements and stabilized using nonlinear entropy viscosity.
  • 3 Numerical methods: The paper focuses on adapting time stepping, compressibility treatment, discontinuous-coefficient averaging, latent-heat incorporation, mesh adaptation, and additional-quantity advection.These modifications target realistic descriptions of convective transport in Earth’s mantle.
  • 3 Numerical methods: All computations use the open-source Aspect mantle-convection code built on deal.II, p4est, and Trilinos.The reported test computations use Aspect version 1.5.0, with setups made available by the authors.

3.1 Time stepping for the temperature equation revisited: explicit, semi-implicit, implicit

Complex, irregular meshes make the stability parameter and cell-size definition in semi-implicit temperature stepping difficult to choose safely. The authors therefore favor fully implicit stepping, which costs more per solve but reduces the number of time steps substantially.

  • 3.1 Time stepping for the temperature equation revisited: explicit, semi-implicit, implicit: Semi-implicit stepping treats thermal diffusion implicitly and advection explicitly, so the CFL condition limits the time step.The bound depends on cell size, temperature-element degree, maximum cell velocity, and method-dependent C.
  • 3.1 Time stepping for the temperature equation revisited: explicit, semi-implicit, implicit: Accuracy can still motivate time steps near the CFL scale because this limits transport to no more than the distance between adjacent nodes.This consideration becomes increasingly important with adaptive mesh refinement.
  • 3.1 Time stepping for the temperature equation revisited: explicit, semi-implicit, implicit: Uniform rectangular or box-shaped meshes permit several equivalent h_K definitions, with safe C values determined experimentally.Equivalent choices include cell diameter, shortest edge, minimum vertex distance, and volume-based measures.
  • 3.1 Time stepping for the temperature equation revisited: explicit, semi-implicit, implicit: On complex meshes, definitions of h_K are no longer equivalent, making the appropriate cell-size measure and stability constant C unclear.The issue arises for varying cell shapes, shell segments, free boundaries, and realistic topologies.
  • 3.1 Time stepping for the temperature equation revisited: explicit, semi-implicit, implicit: Poorly shaped cells can trigger instability only when steep solution features pass through them, requiring extensive testing of C.This makes stability difficult to assess from a single test configuration.
  • 3.1 Time stepping for the temperature equation revisited: explicit, semi-implicit, implicit: C = 1/5.9 in 2d and C = 1/43.6 in 3d were required for stability in the cited earlier work.Larger values produced geometries and setups in which the temperature eventually became unstable.
  • 3.1 Time stepping for the temperature equation revisited: explicit, semi-implicit, implicit: Fully implicit stepping is preferred because semi-implicit stepping needs roughly one sixth as many time steps in 2d and fewer than one fortieth in 3d only by taking much smaller steps.Each semi-implicit step still requires solving both the Stokes and temperature equations.
  • 3.1 Time stepping for the temperature equation revisited: explicit, semi-implicit, implicit: The recommended fully implicit scheme uses C = 1 and h_K equal to the minimum distance between cell vertices.It requires more costly nonsymmetric-system solvers such as GMRES, but reduces time steps by more than 5 in 2d or 40 in 3d.

3.2 Compressibility

Compressibility complicates mantle-convection solvers through nonlinear mass conservation, nonsymmetric systems, discontinuous coefficients, and solver trade-offs. The paper compares reformulations and formulations, showing accurate benchmark results and different cost profiles for stationary and time-dependent problems.

  • 3.2 Compressibility: Compressibility makes mass conservation nonlinear and can produce nonsymmetric Stokes matrices, complicating efficient solver and preconditioner design.These issues arise because density may depend on solution variables and the mass-conservation operator is no longer adjoint to the pressure gradient.
  • 3.2 Compressibility: Reformulating mass conservation with a prescribed or lagged density removes nonlinearity but generally sacrifices symmetry and can require nonlinear iteration when density depends on pressure.The reformulated operator is not adjoint to the pressure-gradient operator, while pressure-dependent density still requires iteration toward the desired solution.
  • 3.2 Compressibility: Discontinuous density makes expanded divergence formulations invalid because derivatives of density and velocity cannot be handled separately with the available pressure regularity.The issue is especially relevant when phase changes create coefficient jumps, and pressure test functions in L2 are insufficiently smooth for integration by parts.
  • 3.2 Compressibility: The proposed correction adjusts compressibility before Stokes assembly to enforce global mass conservation, but it is unnecessary for incompressible materials and undesirable with suitable normal-stress boundaries.For incompressible materials the correction is zero; with an open normal-stress boundary, velocity and pressure can adjust independently instead.
  • 3.2 Compressibility: For stationary computations, the implicit formulation is computationally cheaper and likely more stable, whereas the explicit formulation may be cheaper for time-dependent problems.A single implicit solve has more GMRES iterations, but explicit solves accumulate more iterations across Picard steps; the explicit advantage can emerge when prior-step velocity provides a good approximation.
  • 3.2 Compressibility: The compressible benchmark shows excellent agreement with prior results, while the Stokes-flow benchmark achieves optimal third-order velocity convergence.The reported comparisons include Nusselt number, velocity, temperature, shear heating, and adiabatic heating; heating terms converge at least as fast as velocity in the separate benchmark.

3.3 Averaging of material properties

Abrupt viscosity and density contrasts create pressure oscillations and solver difficulties, motivating cellwise averaging of material properties. Across benchmark tests, averaging improves pressure behavior and computational cost, with harmonic averaging presented as a reasonable accuracy compromise.

  • Abrupt viscosity contrasts cause long solver times, large pressure gradients, and over- and undershoots in discrete pressure approximations.
  • Cellwise averaging replaces quadrature-point material values with averaged or projected representations, but changes the problem and can affect accuracy.
  • Pressure spikes near the discontinuous sinker interface become far less pronounced with averaging, although pressure extrema can continue growing slowly under mesh refinement.
  • All averaging schemes improve pressure approximation, while some deteriorate velocity approximation; harmonic averaging therefore appears to be a reasonable compromise.
  • 5870s to 240s: harmonic averaging reduces Taylor-Hood sinker runtime on a 256 × 256 mesh by roughly a factor of 25.
  • 145 hours to 17 hours: averaging reduces runtime in a roughly 17-million-unknown, 64-processor model without significantly changing the computed solutions.
  • 169 to 77: averaging reduces outer GMRES iterations in the complex model, while inner solves require 2 to 5 times fewer Conjugate Gradient iterations per preconditioner application.

3.4 Latent heat

The paper examines how to represent latent heat released or consumed across sharp phase transitions and how mesh resolution and transition smoothing affect accuracy. A benchmark shows that resolving the transition is essential for predictable convergence.

  • Implementation: Latent heat is represented as a heating term associated with entropy changes as material crosses a phase boundary.The paper discusses analytic phase functions, thermodynamic p-T tables, and effective thermal expansivity approaches.
  • Implementation: Narrow latent-heat zones require smoothing because their widths are often below geodynamic grid resolution and create strong temperature gradients.The smoothing must remain numerically tractable while preserving accuracy.
  • Benchmark setup: The benchmark prescribes inflow at the top, outflow at the bottom, and a central phase transition where latent heat raises the steady-state temperature below the transition.The analytical reference follows from the one-dimensional downward-flow energy equation.
  • Numerical results: Approximately 4 mesh cells across the phase transition are required for small error with bi-quadratic finite elements.When the transition is resolved, the remaining error is dominated by the chosen smoothed transition width.
  • Numerical results: Choosing a phase-transition width of approximately four times the smallest cell size yields the highest accuracy and predictable convergence.If the transition is under-resolved, refinement and reducing the transition width have little significant effect.

3.5 Mesh refinement

The paper presents adaptive mesh refinement as a flexible, goal-dependent process for complex mantle convection models. Standard solution-based indicators may miss discontinuities, while derived-field indicators can resolve interfaces but lack a universally optimal choice.

  • Refinement indicators: The Kelly criterion estimates interpolation error from finite-element derivative jumps and is widely used to mark cells for refinement or coarsening.It can be applied to temperature or velocity fields in mantle convection.
  • Refinement goals: Energy-norm-based refinement may fail to resolve quantities of interest, such as phase-related density, viscosity, or seismic structures.The relevant target may be only tangentially related to the energy-norm error.
  • A practical approach: A two-tiered system combines selectable indicators for solution and derived fields with user-defined choices about the information that matters.Indicators include Kelly estimates, weighted gradients, and gradients of density, viscosity, or thermal energy density.
  • A practical approach: The flexible refinement approach provides high accuracy where needed but has little theoretical underpinning that it is optimal.The method is therefore practical and adaptable rather than supported by a universal optimality guarantee.
  • 2-D spherical convection: In a spherical-shell model with phase transitions, temperature-only refinement largely misses phase boundaries while thermal-energy-density or viscosity gradients resolve them.The model includes discontinuous density and viscosity at 410 km and 660 km depths.
  • 2-D spherical convection: Which mesh is best cannot be quantified without specifying the simulation goal, because interface accuracy and energy-norm accuracy may differ.Meshes resolving tomographically visible interfaces may have larger velocity or temperature energy-norm errors.

3.6 Tracking chemical compositions and other quantities

The paper replaces or supplements tracer particles with advected compositional fields for tracking chemical composition and other transported quantities. The same framework can represent reactions and accumulated finite strain, while requiring stabilization for sharp gradients.

  • Motivation: Complex mantle simulations may need to track composition, trace elements, isotope ratios, water content, grain sizes, and other quantities alongside thermal energy.These quantities are transported with the velocity field.
  • Field-based tracking: Field-based descriptions offer established advection solvers and efficient data movement during adaptive mesh refinement and repartitioning.They also allow material properties to be evaluated at quadrature points.
  • Compositional fields: Compositional fields C_i are advected with source terms Q_i that can represent transport, reactions, partial melting, and freezing.The source terms may depend on velocity, pressure, temperature, and the compositional fields themselves.
  • Numerical challenges: Because the advection equations contain no diffusion, sharp chemical heterogeneities create numerical challenges and require stabilization.The nondiffusive formulation reflects the negligible diffusion of chemical species over mantle length scales.
  • Tracking finite strain: The compositional-field framework tracks accumulated finite strain by storing the components of the deformation gradient tensor F as fields.The velocity gradient updates F, whose polar decomposition yields the stretching tensor used to describe deformation.
  • Tracking finite strain: The finite-strain example shows shortening at plume heads, horizontal stretching around them, and concentrated shear at plume edges, tails, and the bottom boundary layer.The model visualizes temperature together with accumulated natural strain and stretching directions.

4 Application to a complex problem

The paper demonstrates a global mantle convection model combining compressibility, Earth-like material properties, strongly temperature-dependent viscosity, chemical heterogeneities, and prescribed surface velocities. The resulting realistic setup produces narrow plumes and complex material-property structures while retaining the computation’s basic results.

  • Model setup: The global model combines compressibility, Earth-like mineral-physics material properties, strongly temperature-dependent viscosity, compositional fields, and prescribed surface velocities.Its geometry resembles Earth’s mantle and uses an adiabat computed from mineral-property data.
  • Model setup: The simulation uses 23 km finest cells, about 100 million degrees of freedom per time step, and 3,700 time steps.Large portions of the model are adaptively coarsened.
  • Results: Strongly temperature-dependent viscosity produces narrow upwelling plumes about 100 km in diameter that rise from dense basal piles toward hotspot-like surface locations.The passage relates these plume structures to observed hotspot locations and prior studies.
  • Results: Despite greater physical complexity, the model preserves the computation’s basic results while focusing plumes into narrower structures and supporting regional and seismic-model applications.Velocities, temperatures, and compositions can constrain regional models, while compressible profiles facilitate comparisons with seismic models.

5 Conclusions

The paper argues that credible mantle-convection simulations require both more complex physics and larger computational scale. It summarizes numerical methods that jointly enable complex models on large compute clusters and support comparison with observational data.

  • 5 Conclusions: Improving mantle-convection model veracity requires handling highly variable coefficients, latent heat, and adaptively refined meshes at larger discretization scales.The authors frame both physical complexity and computational scale as challenges beyond conventional code operation.
  • 5 Conclusions: The paper summarizes methods for time stepping, compressibility, discontinuous coefficients, latent heat, adaptive finite-element refinement, and advecting additional quantities.These are presented as lessons from solving complex mantle-convection problems with state-of-the-art computational methods.
  • 5 Conclusions: Although no individual technique suffices, their joint use with earlier methods enables some of the most complex mantle-convection models on large-scale compute clusters.The authors expect the methods to remain useful for models incorporating additional material and deformation processes.
  • 5 Conclusions: The resulting models can produce outputs compared with seismic imaging, surface heat fluxes, plate velocities, and other measurements.The passage identifies these measurements as targets for evaluating more complicated models.

A King results

The paper reports compressible benchmark results for both ALA and TALA formulations using increasingly fine meshes. The results include convergence and extrapolated values for comparison with the established benchmark.

  • A King results: Tables 6 and 7 report compressible benchmark results for the ALA and TALA formulations defined by King et al. (2010).The benchmark is widely used and the tables provide a full account of the authors’ results.
  • A King results: The benchmark tables show convergence as mesh size decreases and provide extrapolated values for comparison with values reported by King et al. (2010).This comparison is stated for the compressible-equation strategy discussed in Section 3.2.
  • A King results: The ALA results are obtained on increasingly finer meshes and extrapolated from the 1/128 mesh using Richardson extrapolation.Table 6 identifies this procedure for the ASPECT results.
Loading 1702.05075v2…