Advances in Forward and Inverse Numerical Modelling of Geological Processes


Geological and geophysical data sets are in essence the result of physical processes governing the Earth’s evolution. Such data sets are widely varied and range from the internal structure of the Earth, plate kinematics, composition of geomaterials, estimation of physical conditions, dating of key geological events, thermal state of the Earth to more shallow processes such as natural and “engineered” reservoir dynamics and waste sequestration in the subsurface.

Combining such data with process-based numerical models is required for our understanding of the dynamical Earth. Process-based models are powerful tools to predict the evolution of complex natural systems resolving the feedback among various physical processes. Integrating high-quality data into numerical simulations leads to a constructive workflow to further constrain the key parameters within the models. Innovative inversion strategies, linking forward dynamic models with observables, is therefore an important research topic that will improve our knowledge of the governing physical parameters.

The complexity of geological systems arises from their multi-physics nature, as they combine hydrological, thermal, chemical and mechanical processes (e.g. thermo-mechanical convection). Multi-physics couplings are prone to nonlinear interactions ultimately leading to spontaneous localisation of flow and deformation. Understanding the couplings among those processes therefore requires the development of appropriate tools to capture spontaneous localisation and represents a challenging though essential research direction.

We invite contributions from the following two complementary themes:

#1 Computational advances associated with
- alternative spatial and/or temporal discretisation for existing forward/inverse models
- scalable HPC implementations of new and existing methodologies (GPUs / multi-core)
- solver and preconditioner developments
- AI / Machine learning-based approaches
- code and methodology comparisons (“benchmarks”)
- open source implementations for the community

#2 Physics advances associated with
- development of partial differential equations to describe geological processes
- inversion strategies and adjoint-based modelling
- numerical model validation through comparison with observables (data)
- scientific discovery enabled by 2D and 3D modelling
- utilisation of coupled models to explore nonlinear interactions

Co-organized by EMRP1
Convener: Ludovic Räss | Co-conveners: Thibault Duretz, Boris Kaus, Dave May
vPICO presentations
| Fri, 30 Apr, 13:30–15:00 (CEST)

vPICO presentations: Fri, 30 Apr

Chairpersons: Ludovic Räss, Boris Kaus
Carsten Uphoff, Dave May, and Alice-Agnes Gabriel

Earthquakes and aseismic slip are typically modelled as a displacement discontinuity on a prescribed infinitesimally thin fault surface embedded in linear elastic or viscoelastic media. The fault slip behaviour can be described by laboratory-derived rate and state friction laws, which are suitable to model frictional sliding throughout the complete seismic cycle, i.e. interseismic, coseismic, and post-seismic phase. The governing time scales vary from years in the interseismic phase to seconds in the coseismic phase and the respective spatial scales vary from hundreds of kilometres of tectonic structures  to metres (or less) on-fault. Therefore, simulating the entire seismic cycle is computational challenging and as such mandates utilization of high performance computing (HPC).

We present the open-source code tandem which is designed to model quasi-dynamic sequences of earthquakes and aseismic slip (SEAS). In tandem we explore the usefulness of the symmetric interior penalty Galerkin (SIPG) method using unstructured simplicial meshes for the computation of the elastostatic response to a displacement discontinuity. The potential of the SIPG method for SEAS models lies in (i) its geometric flexibility, (ii) its high-order approximation spaces, (iii) and its natural ability to deal with discontinuities.

Using a number of 2D and 3D SCEC community benchmarks (Erickson et al., 2020) we verify the tandem SIPG implementation. Based on the same reference models, we demonstrate benefits of using highly refined unstructured meshes and a high-order geometric representation of the fault. We also explore whether using a high-order discretisation in space is advantageous. Lastly, we outline how tandem may leverage modern supercomputing resources.

How to cite: Uphoff, C., May, D., and Gabriel, A.-A.: The discontinuous Galerkin method for sequences of earthquakes and aseismic slip, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-5408,, 2021.

Nicolas Riel, Boris Kaus, Eleanor Green, Nicolas Berlie, and Lisa Rummel

During the last decade, the development of numerical geodynamic tools helped the geosciences community to unravel complex thermo-mechanical processes at play during plate tectonics. Yet, the high computational cost of thermodynamic calculations, which simulates phase change, hampers our ability to integrate complex chemistry in such problems. This is particularly important for simulating magmatic processes, where the chemistry of differentiating melts can vary significantly from the mantle to the upper crust. The typical approach, currently used, is to precompute one or many phase diagrams and use them as look-up tables. For many geodynamic processes this is adequate but when the melt chemistry varies drastically it would be better to be able to do thermodynamic calculations on the fly, along with the geodynamic models.

For that, the thermodynamic computational approach must be sufficiently fast, should work fully automatically and be tuned for melting models of magmatic systems, for example by utilizing the recently developed thermodynamic melting model of Holland et al. (2018). Existing approaches do not fulfill all criteria, which is why we have developed a new computational library for this purpose. Our code is written in C, runs on massively parallel machines (MPI) and uses an adaptive mesh refinement strategy to compute phase diagrams. At the moment we have focused on the 'igneous set' of the Holland & Powell dataset (as defined in the thermocalc software) to calculate stable phase equilibria in the system K2O–Na2O–CaO–FeO–MgO–Al2O3–SiO2–H2O–TiO2–Fe2O3–Cr2O3 (KNCFMASHTOCr). The code uses pressure, temperature and bulk-rock composition as input and returns relevant petrological and geodynamic information such as (but not restricted to) stable assemblage, phase fractions and phase densities. Different than many of the existing approaches, our method can efficiently utilize initial guesses which naturally occur in geodynamic simulations where the changes in chemistry between timesteps are usually minor.

The methodology performs a Gibbs free energy minimization and involves two main steps. First, we use a combination of levelling methods (iterative change of base) to reduce the number of potential (pure and solution) phases and to bring the G-hyperplane close to solution. Second, we use a partitioning of Gibbs energy approach coupled with local minimization to satisfy the Gibbs-Duhem rule and to retrieve the final set of stable solution phases. To illustrate the efficiency of the library up to supra-solidus conditions we present a set of dry phase diagrams and compare results of our computations with thermocalc calculations.

Ongoing development includes the treatment of solvus to extend its applicability to complex wet systems involving solution phase such as amphibole.

How to cite: Riel, N., Kaus, B., Green, E., Berlie, N., and Rummel, L.: A new and efficient computational thermodynamics approach for magmatic systems, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-3910,, 2021.

Boris Kaus and Nicolas Riel

Thermodynamics plays an increasingly important role in computational geodynamics, as the community moves towards including chemical reactions in mechanical models of deformation. An example is mineral reactions that induce volume changes which affect the local state of stress that may trigger nonlinear feedbacks. Another example is numerical simulations of magmatic systems where the chemistry of melts changes dramatically as the melt differentiates on its way up through the lithosphere. In order to compare the results of numerical models with geochemical and petrological field data, it is crucial to predict the rock types and major element chemistry from numerical simulations of magmatic systems. Precomputing phase diagrams and including them as a lookup table, which is the current standard approach, works when the chemistry is more or less constant, but is nontrivial for magmatic systems where the chemistry changes drastically and involves a 13-dimensional system (11 oxides, plus pressure and temperature). Parameterizing the main reactions is a possibility, but this renders the comparison of simulation results with available data difficult.

In the context of magmatic systems, significant progress has been made in recent years and we now have thermodynamic melting models that are consistent with experiments and can simulate the full compositional range from mafic to felsic melt compositions. Yet, in order to be useful for geodynamic applications, we also need sufficiently fast Gibbs energy minimization software tools that can automatically determine the most stable assemblage for a given pressure, temperature and chemistry. This is a nonlinear constrained optimization problem, for which we have developed a new, parallel, solution approach. One novelty of our approach is that it can efficiently make use of good initial guesses, for example obtained from the previous timestep of a geodynamic simulation or from nearby points in pressure and temperature space. Yet, as with any gradient-based method, a risk remains to be trapped in a local minimum in the solution space that is not the overall most stable assemblage. It is thus important to explore a sufficiently broad range of starting parameters to ensure convergence toward the global minimum. Whereas this problem is well-known among users of existing thermodynamic software (such as THERMOCALC, where the starting values have to be adjusted manually), it is much less clear how nonlinear the parameter space actually is for real applications. Do we have thousands of local minima, or are there only a few (and if yes, can we precompute some of these)?

Here, we will discuss several examples of melting models and map out the nonlinearity of the parameter space for these cases to get better insights in how to further speed up such calculations. We also discuss how shallow and deep neural networks can be trained and implemented as part of the workflow.

How to cite: Kaus, B. and Riel, N.: Neural networks, local minima and computational thermodynamics, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-5173,, 2021.

Emilie Macherel, Yuri Podladchikov, Ludovic Räss, and Stefan M. Schmalholz

Power-law viscous flow describes well the first-order features of long-term lithosphere deformation. Due to the ellipticity of the Earth, the lithosphere is mechanically analogous to a shell, characterized by a double curvature. The mechanical characteristics of a shell are fundamentally different to the characteristics of plates, having no curvature in their undeformed state. The systematic quantification of the magnitude and the spatiotemporal distribution of strain, strain-rate and stress inside a deforming lithospheric shell is thus of major importance: stress is for example a key physical quantity that controls geodynamic processes such as metamorphic reactions, decompression melting, lithospheric flexure, subduction initiation or earthquakes.

Stress calculations in a geometrically and mechanically heterogeneous 3-D lithospheric shell require high-resolution and high-performance computing. The pseudo-transient finite difference (PTFD) method recently enabled efficient simulations of high-resolution 3-D deformation processes, implementing an iterative implicit solution strategy of the governing equations for power-law viscous flow. Main challenges for the PTFD method is to guarantee convergence, minimize the required iteration count and speed-up the iterations.

Here, we present PTFD simulations for simple mechanically heterogeneous (weak circular inclusion) incompressible 2-D power-law viscous flow in cartesian and cylindrical coordinates. The flow laws employ a pseudo-viscoelastic behavior to optimize the iterative solution by exploiting the fundamental characteristics of viscoelastic wave propagation.

The developed PTFD algorithm executes in parallel on CPUs and GPUs. The development was done in Matlab (, then translated into the Julia language (, and finally made compatible for parallel GPU architectures using the ParallelStencil.jl package ( We may unveil preliminary results for 3-D spherical configurations including gravity-controlled lithospheric stress distributions around continental plateaus.

How to cite: Macherel, E., Podladchikov, Y., Räss, L., and Schmalholz, S. M.: GPU-based pseudo-transient finite difference solution for power-law viscous flow in cartesian, polar and spherical coordinates, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-9108,, 2021.

Nicolas Berlie, Boris Kaus, Anton Popov, Mara Arts, Nicolas Riel, and Daniel Kiss

The dynamics of magmatic systems remain poorly understood, due to the lack of resolving power of geophysical methods to study active systems and the difficulty of interpreting exposed crystallized magma bodies. Numerical models are therefore helpful to connect the dots between classical geological studies, using rheological information and geometries derived from field or geophysical investigations to shed new lights on the mechanisms involved in such systems.

Taking advantage of the big CPU clusters currently available and the development of the DMStag framework as part of the PETSc infrastructure, the ERC-funded MAGMA project aims to build tools to analyse magmatic processes in the lithosphere. We developed a finite-difference staggered grid code solving the Stokes equations for visco-elasto-plastic rheologies and using analytical jacobians for linear and non-linear solvers, combined with regularized plasticity. The code is combined with both a marker and cell and semi-lagrangian advection schemes, is fully parallel and includes automated testing.

Here, we provide application examples ranging from simple benchmark validations against analytical solutions to more complex settings taking advantage of the broad rheologies and local heterogeneities permitted by high resolution settings and the finite difference method. Ongoing technical developments include adding two-phase flow and coupling to it with thermodynamic calculations to track the evolving chemistry of magmatic systems.

How to cite: Berlie, N., Kaus, B., Popov, A., Arts, M., Riel, N., and Kiss, D.: Towards modelling magmatic systems using a staggered grid finite difference method, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-4271,, 2021.

Yuan Li, Adina Pusok, Dave May, and Richard Katz

It is broadly accepted that magmatism plays a key dynamic role in continental and oceanic rifting. However, these dynamics remain poorly studied, largely due to the difficulty of consistently modelling liquid/solid interaction across the lithosphere. The RIFT-O-MAT project seeks to quantify the role of magma in rifting by using models that build upon the two-phase flow theory of magma/rock interaction. A key challenge is to extend the theory to account for the non-linear rheological behaviour of the host rocks, and investigate processes such as diking, faulting and their interaction. Here we present our progress in consistent numerical modelling of poro-viscoelastic–plastic modelling of deformation with a free surface.

Failure of rocks (plasticity) is an essential ingredient in geodynamics models because Earth materials cannot sustain unbounded stresses. However, plasticity represents a non-trivial problem even for single-phase flow formulations (Spiegelman et al. 2016). The elastic deformation of rocks can also affect the propagation of internal failure. Furthermore, deformation and plastic failure drives topographic change, which imposes a significant static stress field. Robustly solving a discretised model that includes this physics presents severe challenges, and many questions remain as to effective solvers for these strongly nonlinear systems. 

We present a new finite difference staggered grid framework for solving partial differential equations (FD-PDE) for single-/two-phase flow magma dynamics (Pusok et al., 2020). Staggered grid finite-difference methods are mimetic, conservative, inf-sup stable and with small stencil — thus they are well suited to address these problems. The FD-PDE framework uses PETSc (Balay et al., 2020) and aims to separate the user input from the discretization of governing equations. The core goals for the FD-PDE framework is to allow for extensible development and implement a framework for rigorous code validation. Here, we present simplified model problems using the FD-PDE framework for two-phase flow visco-elasto-plastic models designed to characterise the solution quality and assess both the discretisation and solver robustness. We also present results obtained using the phase-field method (Sun and Beckermann, 2007) for representing the free surface. Verification of the phase-field approach will be shown via simplified problems previously examined in the geodynamics community (Crameri et al, 2012).

Balay et al. (2020), PETSc Users Manual, ANL-95/11 - Revision 3.13.

Pusok et al. (2020) 

Spiegelman et al. (2016)

Sun and Beckermann (2007)

Crameri et al. (2012)

How to cite: Li, Y., Pusok, A., May, D., and Katz, R.: Free surface models of partially molten rock with visco-elasto–plastic rheology: numerical solutions using a staggered-grid finite-difference discretisation, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-7631,, 2021.

Alice-Agnes Gabriel, Duo Li, Simone Chiocchetti, Maurizio Tavelli, Ilya Peshkov, Evgeniy Romenski, and Michael Dumbser

Earthquake fault zones are more complex, both geometrically and rheologically, than an idealised infinitely thin plane embedded in linear elastic material.  Field and laboratory measurements reveal complex fault zone structure involving tensile and shear fractures spanning a wide spectrum of length scales (e.g., Mitchell & Faulkner, 2009), dense seismic and geodetic recording of small and large earthquakes show hierarchical volumetric faulting patterns (e.g., Cheng et al., 2018, Ross et al., 2019) and 2D numerical models explicitly accounting for off-fault fractures demonstrate important feedback with rupture dynamics and ground motions (e.g., Thomas & Bhat 2018, Okubo et al., 2019).

Here (Gabriel et al., 2021) we adopt a diffuse crack representation to incorporate finite strain nonlinear material behaviour, natural complexities and multi-physics coupling within and outside of fault zones into dynamic earthquake rupture modeling. We use a first-order hyperbolic and thermodynamically compatible mathematical model, namely the GPR model (Godunov & Romenski, 1972; Romenski, 1988),  to describe a continuum in a gravitational field which provides a unified description of nonlinear elasto-plasticity, material damage and of viscous Newtonian flows with phase transition between solid and liquid phases.

The model shares common features with phase-field approaches but substantially extends them. Pre-damaged faults as well as dynamically induced secondary cracks are therein described via a scalar function indicating the local level of material damage (Tavelli et al., 2020); arbitrarily complex geometries are represented via a diffuse interface approach based on a solid volume fraction function (Tavelli et al., 2019). Neither of the two scalar fields needs to be mesh-aligned, allowing thus faults and cracks with complex topology and the use of adaptive Cartesian meshes (AMR). High-order accuracy and adaptive Cartesian meshes are enabled in 2D and 3D by using the extreme scale hyperbolic PDE solver ExaHyPE (Reinarz et al., 2019).

We show a wide range of numerical applications that are relevant for dynamic earthquake rupture in fault zones, including the co-seismic generation of secondary off-fault shear cracks, tensile rock fracture in the Brazilian disc test, as well as a natural convection problem in molten rock-like material. We compare diffuse interface fault models of kinematic cracks, spontaneous dynamic rupture and dynamically generated off-fault shear cracks to sharp interface reference models. To this end, we calibrate the GPR model to resemble empirical tensile and shear crack formation and friction laws. We find that the continuum model can resemble and extend classical solutions, while introducing dynamic differences (i) on the scale of pre-damaged/low-rigidity fault zone, such as out-of- plane rupture rotation; and (ii) on the scale of the intact host rock, such as conjugate shear cracking in tensile lobes. 

Our approach is part of the TEAR ERC project ( and will potentially allow to fully model volumetric fault zone shearing during earthquake rupture, which includes spontaneous partition of fault slip into intensely localized shear deformation within weaker (possibly cohesionless/ultracataclastic) fault-core gouge and more distributed damage within fault rocks and foliated gouges.

How to cite: Gabriel, A.-A., Li, D., Chiocchetti, S., Tavelli, M., Peshkov, I., Romenski, E., and Dumbser, M.: A unified first order hyperbolic model for nonlinear dynamic rupture processes in diffuse fracture zones, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-15237,, 2021.

Mohsen Goudarzi, René de Borst, Taras Gerya, Meng Li, and van Dinther Ylona

Accurate representation of fault zones is important in many applications in Earth sciences, including natural and induced seismicity. The framework developed here can efficiently model fault zone localization, evolution, and spontaneous fully dynamic earthquake sequences in a continuum plasticity framework. The geometrical features of the faults are incorporated into a regularized continuum framework, while the response of the fault zone is governed by a rate and state-dependent friction. Although a continuum plasticity model is advantageous to discrete approaches in representing evolving, unknown, or arbitrarily positioned faults, it is known that either non-associated plasticity or strain-softening can lead to mesh sensitivity of the numerical results in absence of an internal length scale. A common way to regularize the numerical model and introduce an internal length scale is by the adoption of a Kelvin-type visco-plasticity element. The visco-plastic rheological behavior for the bulk material is implemented along with a return-mapping algorithm for accurate stress and strain evolution. High slip rates (in the order of 1 m/s) are captured through numerical examples of a predefined strike-slip fault zone, where a detailed comparison with a reference discrete fault model is presented. Additionally, the regularization effect of the Kelvin viscosity parameter is studied on the fault slip velocity for a growing fault zone due to an initial material imperfection.  The model is consistently linearized leading to quadratic convergence of the Newton solver. Although the proposed framework is a step towards the modeling of earthquake sequences for induced seismicity applications, the numerical model is general and can be applied to all tectonic settings including subduction zones.




How to cite: Goudarzi, M., de Borst, R., Gerya, T., Li, M., and Ylona, V. D.: A regularized continuum model for fault zones with rate and state friction, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-8109,, 2021.

Iskander Muldashev, Marta Pérez-Gussinyé, Mário Neto Cavalcanti de Araújo, and Zhonglan Liu

Rifts and rifted margins result from interaction of several physical processes, which produce a range of crustal structures, subsidence histories, and sedimentary architectures. Study of these processes in academia and industry includes kinematic modelling (i.e. cross-section restoration, backstripping) combined with simple thermomechanical models and dynamic modelling. In kinematic models, the thinning of the lower crust and mantle is kinematically imposed in the form of pure shear, which contradicts natural non-linear viscous behavior. Although, kinematic modelling can provide a crustal thinning profile, heatflow estimates, subsidence rates etc., imposed extension of the lower crust and mantle might strongly impact the result. On the other hand, a dynamic approach allows to model the whole range of possible physical processes, but it cannot be used to model particular extension histories.
Here, we show a new modelling technique, namely KineDyn, to combine the advantages of the above-mentioned approaches into a single modelling framework. Our method employs full non-linear visco-elasto-plastic rheology, surface process of erosion and sediment transport, decompression melting of the mantle, and serpentinization of mantle rocks. Faults are introduced as weak planes in the upper crust, in order to simulate faulting during the model run. In our approach, faults are initially controlled by prescribed initial locations, offsets and timings, while the rest of the model is resolved in a fully dynamic mode. Since fault planes are much weaker than the surrounding upper crust, extension of the model naturally leads to slip on the faults. We demonstrate that faults modelled this way reproduce a natural behavior, including rotation due to flexure and unloading of the fault plane.
In order to reconstruct the evolution of an existing rift or rifted margin we model extension of the lithosphere with controlled faulting. To do this we use the interpreted spatio-temporal evolution of the faulting from a seismic profile to guide the evolution of the dynamic model. After a trial-and-error process, where we correct the faults’ locations, the thicknesses of layers, surface process’s parameters, initial thermal gradient etc., we obtain the model that best fits the observations. Thus, KineDyn gives, in effect, the same results as existing section restoration techniques (i.e. the potential history of faulting) and forward modeling techniques (i.e. the likely history of sedimentation, thinning, heat flow and subsidence), while simultaneously taking into account non-linear interactions between processes occurring during rifting.
In this work we show the methodology, examples, tests and benchmarks of the technique. Finally, we present applications of KineDyn for the following rifts and rifted margins: Malawi Rift, East African Rift System, hyper-extended West Iberia Margin, and ultra-wide Santos-Benguela Rifted Margin.

How to cite: Muldashev, I., Pérez-Gussinyé, M., Neto Cavalcanti de Araújo, M., and Liu, Z.: KineDyn: Thermo-Mechanical Method for Validation of Seismic Interpretations with Dynamic Forward Models, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-8541,, 2021.

Taras Gerya, Thibault Duretz, and Rass Ludovic

In the marker-in-cell method combined with staggered finite differences, Lagrangian markers carrying information on material properties are advected with the velocity field interpolated from the staggered Eulerian velocity grid. With existing schemes, velocity interpolation from the grid points to markers violates (to some extent) mass conservation requirement that causes excess convergence/divergence of markers and opening marker gaps after significant amount of advection. This effect is especially well visible in case of diagonal simple shear deformation along planes that are oriented at 45 degrees to the grid and marker circulation through grid corners.

Here, we present a new second order velocity interpolation scheme that guarantees exact interpolation of normal strain rate components from pressure nodes (i.e. from the locations where these components are defined by solving of the mass conservation equation). This new interpolation scheme is thus applicable to both compressible and incompressible flow and is trivially expendable to 3D and to non-regular staggered grids.

Performed tests show that, compared to other velocity interpolation approaches, the new scheme has superior performance in preserving continuity of the marker field during the long-term advection including the diagonal simple shear deformation and marker circulation through grid corners. We showcase a performance-oriented implementation of the new scheme using Julia language's shared memory parallelisation features. The Julia implementation of the new advection schemes further augments the ParallelStencil.jl  related application collection with advection routines.

How to cite: Gerya, T., Duretz, T., and Ludovic, R.: New continuity-based velocity interpolation scheme for staggered grids, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-15308,, 2021.

Patrick Sanan

Many problems in computational geophysics can be answered by solving a discrete problem on a staggered grid - these define quantities on the elements, faces, edges, and vertices of a regular cell complex, embedded in space with an orthogonal coordinate system. A key example is the Stokes equations for creeping or highly-viscous fluid flow in Cartesian or spherical polar coordinates. These discretizations are attractive in many ways, both practically and computationally. They provide an intuitive way to pose the discrete problem, as they can be motivated as low-order finite volume methods or even as finite difference methods. They provide very regular data layout and relative easy of implementation, allowing for optimized and portable implementations in software. They couple particularly well with particle systems (as in the classical MAC scheme) which define material coefficients with arbitrary, non-grid-aligned interfaces; they provide first-order convergence with an extremely narrow stencil - in solving the discrete system, unknowns are coupled to fewer neighbors than in higher-order finite element or finite volume methods. However, especially in the context of geodynamical and seismic simulations, a severe limitation exists. Many open scientific questions relate to resolving localized features like shear bands, faults, subduction zone boundaries, and entrainment zones, and other boundary layers. With first order methods on regular grids, to resolve these fine spatial features, one must uniformly refine the entire domain (or resort to "Swiss cross" meshes, which are limited in the geometry they can resolve). These quickly becomes computationally intractable. One approach to resolving this tension is to define and use a non-uniform mesh. This is the basis for adaptive mesh refinement (AMR) methods, which are well-established in the context of finite element methods. Here, we present ongoing work in defining a stable, non-uniform, analogue to the MAC scheme for the 3D Stokes equations, attempting to retain the narrow stencil width characteristic of the uniform method while supporting arbitrary 2-to-1 refinement in a mesh described by an octree. Several researchers have made steps in this direction, so in this contribution we review those methods and discuss how they may be generalized or modified to provide a 3D scheme suitable for use in variable-viscosity Stokes flow.

How to cite: Sanan, P.: Towards locally refined 3D staggered grids for the Stokes equations, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-14305,, 2021.

Paul Tackley

In order to treat a free surface in models of lithosphere and mantle dynamics that use a fixed Eulerian grid it is typical to use "sticky air", a layer of low-viscosity, low-density material above the solid surface (e.g. Crameri et al., 2012). This can, however, cause numerical problems, including poor solver convergence due to the huge viscosity jump and small time-steps due to high velocities in the air. Additionally, it is not completely realistic because the assumed viscosity of the air layer is typically similar to that of rock in the asthenosphere so the surface is not stress free.  

In order to overcome these problems, Duretz et al. (2016) introduced and tested a method for treating the free surface that instead detects and applies special conditions at the free surface. This avoids the huge viscosity jump and having to solve for velocities in the air. They applied it to a two-dimensional staggered grid finite difference / finite volume scheme, a discretization that is in common use for modelling mantle and lithosphere dynamics. Here I document the application of this approach to a three-dimensional spherical staggered grid solver in the mantle simulation code StagYY. Some adjustments had to be made to the two-dimensional scheme documented in Duretz et al. (2016) in order to avoid problems due to undefined velocities for certain boundary topographies. The approach was applied not only to the Stokes solver but also to the temperature solver, including the implementation of a mixed radiative/conductive boundary condition applicable to surface magma oceans/lakes.


Crameri, F., H. Schmeling, G. J. Golabek, T. Duretz, R. Orendt, S. J. H. Buiter, D. A. May, B. J. P. Kaus, T. V. Gerya, and P. J. Tackley (2012), A comparison of numerical surface topography calculations in geodynamic modelling: an evaluation of the ‘sticky air’ method, Geophysical Journal International,189(1), 38-54, doi:10.1111/j.1365-246X.2012.05388.x.

Duretz, T., D. A. May, and P. Yamato (2016), A free surface capturing discretization for the staggered grid finite difference scheme, Geophysical Journal International, 204(3), 1518-1530, doi:10.1093/gji/ggv526.

How to cite: Tackley, P.: A Three-Dimensional Version of the Free Surface Capturing Discretization for Staggered Grid Finite Difference Schemes: Implementation into StagYY, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-16228,, 2021.

Sergio Zlotnik, Olga Ortega, Pedro Díez, and Juan Carlos Afonso
One of the main challenges in modern geophysics is the understanding and characterization of the present-day physical state of the thermal and compositional structure of the Earth’s lithospheric and sub-lithospheric mantle. In doing so, high resolution inverse problems need to be solved (with thousands of parameters to determine).
One of the most abundant and better constrained data used for the inversion is the Earth’s topography. Despite its quality, the topography models included in inversion schemes are usually very simplistic, based on density contrasts and neglecting dynamic components. The reason for this is simply computational efficiency; 3D dynamical models are too expensive to be embedded within inversion schemes.
In this context we propose the use of a greedy reduced basis strategy within a probabilistic Bayesian inversion scheme (MCMC) that makes feasible accounting for the fully dynamic topography model within the inversion.
We tested the proposed approach in a synthetic experiment aiming to recover the base of the African plate. It is well-agreed within the geophysical community that the dynamic component in the region is of first order importance. Our scheme is able to successfully recover the expected shape of the plate while reducing the computational time to less than 1% when compared to a full Finite Element approach.

How to cite: Zlotnik, S., Ortega, O., Díez, P., and Afonso, J. C.: A combined Reduced Order-Bayesian scheme to drastically accelerate stochastic inversions, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-8824,, 2021.

Dave May and Philip England

Subduction zones can give rise to severe natural hazards, e.g. earthquakes, tsunami & volcanism. Improved hazard assessment may be realised through physics based modelling. The thermal structure of a subducting plate has a first order control on many aspects of the subduction zone, including: dehydration reactions; intermediate depth seismicity; melt production; formation of arc volcanoes. Subduction zones exhibit a wide variability with respect to slab age, velocity, dip, rheology and mechanical behaviour of the overriding plate. For many subduction zones the assumption of a thermo-mechanical steady-state is reasonable, hence forward models often assume the form of a kinematically driven slab causing traction-driven mantle wedge flow. Even for this simplified forward model, our understanding of how the parameters and their uncertainties influence the thermal structure is incomplete. 

To address this uncertainty, here we use a data-driven model reduction technique, specifically the interpolated Proper Orthogonal Decomposition (iPOD), to define a fast-to-evaluate and surrogate model of a steady-state subduction zone that is valid over a high-dimensional parameter space. The accuracy of the iPOD surrogate model is controlled using a hyper-rectangle tree-based adaptive sampling strategy combined with a non-intrusive error estimator. To illustrate the applicability of the iPOD, we present examples in which reduced-order models are constructed for combinations of parameters related to the kinematics, rheology and geometry of the subduction zone. The examples will characterize the efficiency and accuracy of the iPOD reduced-order model when using parameter spaces that vary in dimension from 1 to 7.

How to cite: May, D. and England, P.: Quantifying thermal variability in subduction zones via data-driven reduced-order modelling, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-16138,, 2021.

Siddhant Agarwal, Nicola Tosi, Pan Kessel, Sebastiano Padovan, Doris Breuer, and Grégoire Montavon

The thermal evolution of terrestrial planets depends strongly on several parameters and initial conditions that are poorly constrained. Often, direct or indirect observables from planetary missions such as elastic lithospheric thickness, crustal thickness and duration of volcanism are inverted to infer the unknown parameter values and initial conditions. The non-uniqueness and non-linearity of this inversion necessitates a probabilistic inversion framework. However, due to the expensive nature of forward dynamic simulations of thermal convection , Markov Chain Monte Carlo methods are rarely used. To address this shortcoming, some studies have recently shown the effectiveness of Mixture Density Networks (MDN) (Bishop 1995) in being able to approximate the posterior probability using only the dataset of simulations run prior to the inversion (Meier et al. 2007, de Wit et al. 2013, Käufl et al. 2016, Atkins et al. 2016).

Using MDNs, we systematically isolate the degree to which a parameter can be constrained using different “present-day” synthetic observables from 6130 simulations for a Mars-like planet. The dataset – generated using the mantle convection code GAIA (Hüttig et al. 2013)- is the same as that used by Agarwal et al. (2020) for a surrogate modelling study.

The loss function used to optimize the MDN (log-likelihood) provides a single robust quantity that can be used to measure how well a parameter can be constrained. We test different numbers and combinations of observables (heat flux at the surface and core-mantle boundary, radial contraction, melt produced, elastic lithospheric thickness, and duration of volcanism) to constrain the following parameters: reference viscosity, activation energy and activation volume of the diffusion creep rheology, an enrichment factor for radiogenic elements in the crust, and initial mantle temperature. If all observables are available, reference viscosity can be constrained to within 32% of its entire range (1019−1022 Pa s), crustal enrichment factor (1−50) to within 15%, activation energy (105−5×105 J mol-1 ) to within 80%, and initial mantle temperature (1600−1800K) to within 39%. The additional availability of the full present-day temperature profile or parts of it as an observable tightens the constraints further. The activation volume (4×10-6 −10×10-6  m3 mol-1) cannot be constrained and requires research into new observables in space and time, as well as fields other than just temperature. Testing different levels of uncertainty (simulated using Gaussian noise) in the observables, we found that constraints on different parameters loosen at different rates, with initial temperature being the most sensitive. Finally, we present how the marginal MDN proposed by Bishop (1995) can be modified to model the joint probability for all parameters, so that  the inter-parameter correlations and the associated degeneracy can be capture, thereby providing a more comprehensive picture of all the evolution scenarios that fit given observational constraints.

How to cite: Agarwal, S., Tosi, N., Kessel, P., Padovan, S., Breuer, D., and Montavon, G.: Towards constraining Mars' thermal evolution using machine learning, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-4044,, 2021.

Arne Spang, Tobias S. Baumann, and Boris J. P. Kaus

Advanced numerical methods and increasing computational power have enabled us to incorporate numerical forward models into geodynamic inverse frameworks. We now have several strategies to constrain the rheological properties of the crust and lithosphere. Yet, the initial geometry of geological formations (e.g., salt bodies, magma bodies, subducting slabs) and associated uncertainties are, in most cases, excluded from the inverse problem and assumed to be part of the a priori knowledge. Usually, geometrical properties remain constant, or we employ simplified bodies like planes, spheres or ellipsoids for their parameterization.

Here, we present a simple method to parameterize complex three-dimensional bodies and incorporate them into geodynamic inverse problems. Our approach enables us to automatically create an entire ensemble of initial geometries within the uncertainty limits of geophysical imaging data. This not only allows us to account for geometric uncertainties, but also enables us to investigate the sensitivity of geophysical data to the geometrical properties of the geological structures.

We present 3 areas of application for our method, covering salt diapirs, magmatic systems and subduction zones, using both synthetic and real data.

How to cite: Spang, A., Baumann, T. S., and Kaus, B. J. P.: Geodynamic inversion with uncertain initial geometries, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-5680,, 2021.

Wolfgang Szwillus

Commonly, the physical properties of the Earth (e.g., velocity, density) are parameterized as continuous fields. The most popular representation are grids and basis functions like spherical harmonics or splines. In an inversion context it is quite common that not all the parameters are fully constrained by the available inputdata. This relates to the common issues of insufficient resolution, incomplete coverage, and trade-offs due tonon-uniqueness. By applying some form of regularization to the inverse problem, a well-behaved and unique solution can be obtained, but this solution depends on the details of the chosen regularization.

Transdimensional approaches address the regularization problem by using a model representation with a variable number of parameters. The number of parameters is adjusted according to the requirements of the input data using the reversible jump Monte Carlo Markov Chain (rj-MCMC) algorithm. The output is an ensemble of variable resolution models that provides insight into the required model complexity and trade-offbetween parameters.

Here, I present synthetic tests from a joint inversion of satellite gravity gradients and normal modes for the Earth's velocity and density structure. The mantle's seismic velocity and density inside a 2-D spherical annulus are described by a variable number of discrete anomalous volumes, each with a variable size, shape, location and strength of velocity and density anomaly. The discrete anomalies are adjusted using the transdimensional approach in order to fit the gravity and normal mode data. This synthetic example shows promising results, because the synthetic model can recovered reasonably well.

How to cite: Szwillus, W.: Possibilities of transdimensional inversion for estimating deep Earth velocity and mantle structure, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-14782,, 2021.

Hugo Bloem, Daniel Tetzlaff, and Andrew Curtis

Seismic tomography is used widely to image the Earth's interior structure and to infer subsurface properties. Tomography is a nonlinear inverse problem, so computationally expensive inversion methods are required to estimate uncertainties in tomographic results. Monte Carlo sampling explores the Bayesian posterior probability distribution function (pdf) which describes the solution to tomographic problems. However, this is expensive, and often intractable for high dimensional model spaces due to the curse of dimensionality – the exponential increase in computation required to explore parameter space with increasing number of degrees of freedom in the Earth model. Variational methods and some neural network inversion methods use optimisation to estimate the posterior pdf. These methods may be more efficient, but they still suffer from the ‘curse’ in high dimensional problems.

The Bayesian solution to tomographic problems combines information available prior to collecting the current data set with information from the geophysical data. To counteract the curse we wish to inject geological prior information that reduces the hypervolume of parameter space to be explored. We use geological process modelling software SedSimple to generate geological training images. This software is designed specifically to produce relatively simple three-dimensional geological structures at reduced computational cost compared to more sophisticated current process models. The output represents the general form of expected geological structures, but not specific detail that might be encountered in any particular volume of the Earth. These geological models are used to train a generative adversarial network (GAN): the GAN then performs high-dimensional regression between the models so that it can generate other, similar models extremely rapidly. The latent space of the trained GAN provides a reduced-dimensionality representation of prior information from SedSimple, which we wish to use as prior information to constrain geophysical imaging problems.

The GAN provides a mapping from latent parameters a to simplistic geological models. Real structure consists of a simplistic model, plus overlying geological complexity parametrised by vector m. We seek the posterior probability of m and a given geophysical data d, written . Assume d is divided into a part ds that is only sensitive to simplistic structures, the remaining data being mainly sensitive to the overlying complexity (e.g. wave travel times versus seismic waveforms). We evaluate the posterior pdf P’(ads) of latent parameters a given data ds using the GAN. This pdf is estimated without considering trade-offs between a and m since we limit the data to ds only.

We can expand the posterior pdf using identity P(m,ad) = P(ma,d) P(ad). Assuming the influence of m on a in the posterior is minimal, we can approximate this equation by P(m,ad) = P(ma,d) P’(ads). We can estimate P(ma,d) for each value of latent parameters a by using conventional linearised methods if m is very high dimensional, or using fully nonlinear methods if m can be decomposed into lower-dimensional dependencies on the data. This framework allows informative prior information from geological process models to constrain posterior pdf’s on the full complexity model (m,a) using geophysical data d, as will be illustrated in this talk.

How to cite: Bloem, H., Tetzlaff, D., and Curtis, A.: Geological Prior Information for Bayesian Tomography, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-7292,, 2021.

Andrea Tomassi, Fabio Trippetta, Roberto De Franco, and Roberta Ruggieri

Forward modeling is a fundamental prospecting method to understand reservoirs structure and architecture in the subsurface. Moreover, it also helps maximizing the hydrocarbons extraction of petroleum systems. Carbonate reservoirs result in non-univocal seismic response caused by the facies heterogeneity and due to the possible presence of infilling fluids. The carbonate ramp outcropping in the Majella Massif (Central Italy) is an excellent surface analogue of buried structures. It offers the opportunity to directly analyze a carbonate reservoir which clearly shows facies variations and natural hydrocarbon-impregnations allowing to quantify the induced petrophysical changes. In this study field and laboratory measurements are used to carry out 1D and 2D forward seismic models of the reservoir. Density and porosity of samples were measured through a helium pycnometer on both hydrocarbon-saturated and not-saturated natural samples. Data show density from 2,35 g/cm3 to 2,64 g/cm3 and porosity from 5,9% to 21%. Seismic velocity data from previous works and ongoing measurements show values from 3,24 km/s to 5,93 km/s and they are clearly related to porosity. After building a seismic velocity, density and acoustic impedance model, a low-frequency (40Hz) synthetic 1D seismogram was carried out simulating facies (porosity) and hydrocarbon-saturation variations. The presence of hydrocarbon causes an increase in acoustics impedance by 16,2% and it increases the amplitudes by 13,6% at the reservoir boundaries. A 12 km long synthetic profile from the platform top to the basin, oriented SSE-NNW, was then carried out simulating the outcropping architecture and spatial distribution of the facies. The synthetic 2D model demonstrates how Vp variations related to facies association and distribution slightly influence seismic facies. Sensitivity tests show that Ricker wavelet with 25 Hz of central frequency is the frequency that produces best seismic images. Ongoing laboratory measurements and simulations are expected to help in validating facies interpretations of real seismic profiles and detect the presence of hydrocarbons in carbonate-ramp petroleum systems.

How to cite: Tomassi, A., Trippetta, F., De Franco, R., and Ruggieri, R.: The impact of facies heterogeneity on the seismic properties of carbonates: forward modeling and reservoirs potential, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-12865,, 2021.

Jose Bastias, Gabriel Rau, Andy Wilkins, and Philipp Blum

Realistic modelling of tightly coupled hydro-geomechanical processes is relevant for the assessment of many hydrological and geotechnical applications. Such processes occur in geological formations and are influenced by natural heterogeneities. Current numerical libraries offer capabilities and physics coupling, that have proven to be valuable in simulating various applications in geotechnical fields such as underground gas storage, rock fracturing, land subsidence and Earth resources extraction. However, implementation and verification of full heterogeneity of subsurface properties using high resolution field data in coupled simulations has not been done yet. Hence, we develop, verify and document RHEA (Real HEterogeneity App), an open-source fully coupled finite element application capable of including node-resolution hydro-geomechanical properties in coupled simulations. We propose a simple, yet powerful workflow to allow the integration of fully distributed hydro-geomechanical properties. We then verify the code with analytical solutions in one and two dimensions and propose a benchmark semi-analytical problem to verify heterogeneous systems with sharp gradients. Finally, we exemplify RHEA's capabilities with a comprehensive example integrating realistic properties. With this we demonstrate that RHEA is a verified open-source application able to integrate complex geology to perform scalable fully coupled hydro-geomechanical simulations. Our work is a valuable tool to assess real world hydro-geomechanical challenging systems that may include different levels of complexity like heterogeneous geology with several time and spatial scales and sharp gradients produced by contrasting subsurface properties.

How to cite: Bastias, J., Rau, G., Wilkins, A., and Blum, P.: RHEA: A verified simulator for hydro-geomechanical heterogeneity, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-3440,, 2021.