skip to main content


Title: 101 geodynamic modelling: how to design, interpret, and communicate numerical studies of the solid Earth
Abstract. Geodynamic modelling provides a powerful tool to investigate processes in the Earth's crust, mantle, and core that are not directly observable. However, numerical models are inherently subject to the assumptions and simplifications on which they are based. In order to use and review numerical modelling studies appropriately, one needs to be aware of the limitations of geodynamic modelling as well as its advantages. Here, we present a comprehensive yet concise overview of the geodynamic modelling process applied to the solid Earth from the choice of governing equations to numerical methods, model setup, model interpretation, and the eventual communication of the model results. We highlight best practices and discuss their implementations including code verification, model validation, internal consistency checks, and software and data management. Thus, with this perspective, we encourage high-quality modelling studies, fair external interpretation, and sensible use of published work. We provide ample examples, from lithosphere and mantle dynamics specifically, and point out synergies with related fields such as seismology, tectonophysics, geology, mineral physics, planetary science, and geodesy. We clarify and consolidate terminology across geodynamics and numerical modelling to set a standard for clear communication of modelling studies. All in all, this paper presents the basics of geodynamic modelling for first-time and experienced modellers, collaborators, and reviewers from diverse backgrounds to (re)gain a solid understanding of geodynamic modelling as a whole.  more » « less
Award ID(s):
1925677
NSF-PAR ID:
10347047
Author(s) / Creator(s):
; ; ; ; ;
Date Published:
Journal Name:
Solid Earth
Volume:
13
Issue:
3
ISSN:
1869-9529
Page Range / eLocation ID:
583 to 637
Format(s):
Medium: X
Sponsoring Org:
National Science Foundation
More Like this
  1. SUMMARY

    Phase transitions play an important role for the style of mantle convection. While observations and theory agree that a substantial fraction of subducted slabs and rising plumes can move through the whole mantle at present day conditions, this behaviour may have been different throughout Earth’s history. Higher temperatures, such as in the early Earth, cause different phase transitions to be dominant, and also reduce mantle viscosity, favouring a more layered style of convection induced by phase transitions. A period of layered mantle convection in Earth’s past would have significant implications for the secular evolution of the mantle temperature and the mixing of mantle heterogeneities. The transition from layered to whole mantle convection could lead to a period of mantle avalanches associated with a dramatic increase in magmatic activity. Consequently, it is important to accurately model the influence of phase transitions on mantle convection. However, existing numerical methods generally preclude modelling phase transitions that are only present in a particular range of pressures, temperatures or compositions, and they impose an artificial lower limit on the thickness of phase transitions. To overcome these limitations, we have developed a new numerical method that solves the energy equation for entropy instead of temperature. This technique allows for robust coupling between thermodynamic and geodynamic models and makes it possible to model realistically sharp phase transitions with a wide range of properties and dynamic effects on mantle processes. We demonstrate the utility of our method by applying it in regional and global convection models, investigating the effect of individual phase transitions in the Earth’s mantle with regard to their potential for layering flow. We find that the thickness of the phase transition has a bigger influence on the style of convection than previously thought: with all other parameters being the same, a thin phase transition can induce fully layered convection where a broad phase transition would lead to whole-mantle convection. Our application of the method to convection in the early Earth illustrates that endothermic phase transitions may have induced layering for higher mantle temperatures in the Earth’s past.

     
    more » « less
  2. null (Ed.)
    SUMMARY Evidence from seismology, geology and geodynamic studies suggests that regional-scale lower crustal flow occurs in many tectonic settings. Pressure gradients caused by mantle processes and crustal density heterogeneity can provide driving force for lower crustal flow. Numerically modelling such flow can be computationally expensive. However, by exploiting symmetry in the physical system, it is possible to represent the vertical component of flow in terms of its lateral components, thereby reducing the problem’s spatial dimension by one. Here, we present a mathematical formulation for flow in a viscous channel below an elastic upper plate, which is optimized for solution by common numerical methods. Our formulation drastically reduces the computational load required to simulate lower crustal flow over large areas and long timescales. We apply this model to two example problems, with and without an elastic upper plate, identifying combinations of parameters that are capable of generating measurable geologic uplift. 
    more » « less
  3. null (Ed.)
    Seismic observations indicate that the lowermost mantle above the core-mantle boundary is strongly heterogeneous. Body waves reveal a variety of ultra-low velocity zones (ULVZs), which extend not more than 100 km above the core-mantle boundary and have shear velocity reductions of up to 30 per cent. While the nature and origin of these ULVZs remain uncertain, some have suggested they are evidence of partial melting at the base of mantle plumes. Here we use coupled geodynamic/thermodynamic modelling to explore the hypothesis that present-day deep mantle melting creates ULVZs and introduces compositional heterogeneity in the mantle. Our models explore the generation and migration of melt in a deforming and compacting host rock at the base of a plume in the lowermost mantle. We test whether the balance of gravitational and viscous forces can generate partially molten zones that are consistent with the seismic observations. We find that for a wide range of plausible melt densities, permeabilities and viscosities, lower mantle melt is too dense to be stirred into convective flow and instead sinks down to form a completely molten layer, which is inconsistent with observations of ULVZs. Only if melt is less dense or at most ca. 1 per cent more dense than the solid, or if melt pockets are trapped within the solid, can melt remain suspended in the partial melt zone. In these cases, seismic velocities would be reduced in a cone at the base of the plume. Generally, we find partial melt alone does not explain the observed ULVZ morphologies and solid-state compositional variation is required to explain the anomalies. Our findings provide a framework for testing whether seismically observed ULVZ shapes are consistent with a partial melt origin, which is an important step towards constraining the nature of the heterogeneities in the lowermost mantle and their influence on the thermal, compositional, and dynamical evolution of the Earth. 
    more » « less
  4. SUMMARY

    Combining finite element methods for the incompressible Stokes equations with particle-in-cell methods is an important technique in computational geodynamics that has been widely applied in mantle convection, lithosphere dynamics and crustal-scale modelling. In these applications, particles are used to transport along properties of the medium such as the temperature, chemical compositions or other material properties; the particle methods are therefore used to reduce the advection equation to an ordinary differential equation for each particle, resulting in a problem that is simpler to solve than the original equation for which stabilization techniques are necessary to avoid oscillations.

    On the other hand, replacing field-based descriptions by quantities only defined at the locations of particles introduces numerical errors. These errors have previously been investigated, but a complete understanding from both the theoretical and practical sides was so far lacking. In addition, we are not aware of systematic guidance regarding the question of how many particles one needs to choose per mesh cell to achieve a certain accuracy.

    In this paper we modify two existing instantaneous benchmarks and present two new analytic benchmarks for time-dependent incompressible Stokes flow in order to compare the convergence rate and accuracy of various combinations of finite elements, particle advection and particle interpolation methods. Using these benchmarks, we find that in order to retain the optimal accuracy of the finite element formulation, one needs to use a sufficiently accurate particle interpolation algorithm. Additionally, we observe and explain that for our higher-order finite-element methods it is necessary to increase the number of particles per cell as the mesh resolution increases (i.e. as the grid cell size decreases) to avoid a reduction in convergence order.

    Our methods and results allow designing new particle-in-cell methods with specific convergence rates, and also provide guidance for the choice of common building blocks and parameters such as the number of particles per cell. In addition, our new time-dependent benchmark provides a simple test that can be used to compare different implementations, algorithms and for the assessment of new numerical methods for particle interpolation and advection. We provide a reference implementation of this benchmark in aspect (the ‘Advanced Solver for Problems in Earth’s ConvecTion’), an open source code for geodynamic modelling.

     
    more » « less
  5. SUMMARY EarthScope's USArray seismic component provided unprecedented coverage of the contiguous United States and has therefore spurred significant advances in tomographic imaging and geodynamic modelling. Here, we present a new global, radially anisotropic shear wave velocity tomography model to investigate upper mantle structure and North American Plate dynamics, with a focus on the contiguous United States. The model uses a data-adaptive mesh and traveltimes of both surface waves and body waves to constrain structure in the crust and mantle in order to arrive at a more consistent representation of the subsurface compared to what is provided by existing models. The resulting model is broadly consistent with previous global models at the largest scales, but there are substantial differences under the contiguous United States where we can achieve higher resolution. On these regional scales, the new model contains short wavelength anomalies consistent with regional models derived from USArray data alone. We use the model to explore the geometry of the subducting Farallon Slab, the presence of upper mantle high velocity anomalies, low velocity zones in the central and eastern United States and evaluate models of dynamic topography in the Cordillera. Our models indicate a single, shallowly dipping, discontinuous slab associated with the Farallon Plate, but there are remaining imaging challenges. Inferring dynamic topography from the new model captures both the long-wavelength anomalies common in global models and the short-wavelength anomalies apparent in regional models. Our model thus bridges the gap between high-resolution regional models within the proper uppermost mantle context provided by global models, which is crucial for understanding many of the fundamental questions in continental dynamics. 
    more » « less