GD10.1 | Advances in Numerical Modelling of Geological Processes: Methods, Applications and Tooling
Advances in Numerical Modelling of Geological Processes: Methods, Applications and Tooling
Convener: Ludovic Räss | Co-conveners: Boris Kaus, Ivan Utkin, Thibault Duretz, Dave May
| Wed, 17 Apr, 16:15–18:00 (CEST)
Room -2.21
Posters on site
| Attendance Wed, 17 Apr, 10:45–12:30 (CEST) | Display Wed, 17 Apr, 08:30–12:30
Hall X1
Orals |
Wed, 16:15
Wed, 10:45
Geological and geophysical data sets convey observations 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 in the subsurface.

The complexity in the physics of geological processes arises from their multi-physics nature, as they combine hydrological, thermal, chemical and mechanical processes. 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 numerical tools.

Integrating high-quality data into physics-based predictive numerical simulations may lead to a constructive workflow and allow to further constrain unknown key parameters within the models. Innovative inversion strategies, linking forward dynamic models with observables, and combining PDE solvers with machine-learning via differentiable programming is therefore an important research topic that will improve our knowledge of the governing physical parameters.

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
- Combining PDEs with AI / Machine learning-based approaches (physics-informed ML)
- Automatic differentiation (AD) and differentiable programming
- code and methodology comparisons (“benchmarks”)

#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

We plan on a special topical issue for the research output presented in this session in the EGU journal of Geoscientific Model Development (GMD)

Orals: Wed, 17 Apr | Room -2.21

Chairpersons: Ludovic Räss, Thibault Duretz
Multi-physics modelling
On-site presentation
Arne Spang, Marcel Thielmann, Daniel Kiss, and Casper Pranger

Thermal runaway is a ductile weakening mechanism that describes the feedback loop of shear heating, temperature-dependent viscosity and localization. It has been linked to deep-focus earthquakes which are unlikely to be caused by brittle failure due to the large lithostatic pressure.

We present one- and two-dimensional (1D and 2D) numerical, thermomechanical models that investigate the occurrence, nucleation and temporal evolution of thermal runaway in a simple shear setting. The models are characterized by a visco-elastic rheology where viscous creep is accommodated with a composite rheology of diffusion and dislocation creep as well as low-temperature plasticity. We implement the model in the Julia programming language and utilize the pseudo-transient iterative method to solve this nonlinear system of equations. Graphical processing unit (GPU) computing (by making use of the package ParallelStencil.jl) allows us to achieve high resolution models in 2D.

Like brittle plasticity, thermal runaway presents the challenge of modeling a very thin shear band in a continuum mechanics approach with finite resolution. To address this issue, we tested two different regularization techniques to provide stable solutions and introduce a grid-independent shear zone width. A viscosity regularization and a second-order gradient regularization achieve similar results. Loading and heating time scales of thousands of years in combination with relaxation time scales on the order of seconds also require an adaptive time stepping scheme that can span 12 orders of magnitude. We achieve this by adjusting time steps to the maximum gradients in temperature and stress, and by dynamically rescaling variables during computation to minimize rounding errors. Combining the aforementioned techniques allows us to cover loading and heating on geological time scales as well as near-instantaneous stress drop and local temperature surge in one model framework.

2D experiments show that thermal runaway allows highly localized ductile ruptures to nucleate at small heterogeneities and propagate like brittle fractures. The ruptures accelerate during propagation and reach the highest velocities when two tips link up. Rupture trajectories are usually parallel to the direction of background deformation but bend in the vicinity of other ruptures to allow for a link up.

How to cite: Spang, A., Thielmann, M., Kiss, D., and Pranger, C.: Thermal runaway and the challenges of rapid localization, EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-18807,, 2024.

On-site presentation
Sergio Zlotnik, Mariano Tomás Fernández, and Pedro Díez

Inverse problems in geophysics seeking to understand the current state of planet Earth use data from multiple observables and involve a variety of physical principles. One of the key fields to determine is the temperature, affecting almost all the other physical quantities involved (e.g. densities, viscosities, wave propagation velocities, among others).

The Lithosphere-Asthenosphere Boundary (LAB) is a boundary layer that affects most of the processes and properties in the Earth structure. Determining its location is one of the goals of inversions. It is usual within numerical studies to define the LAB as an isotherm. The need of determining the thermal field in accordance with a given LAB location leads to a mathematical problem with imposed interior Dirichlet conditions. In particular, the isotherm defining the LAB has to be located in the position tested by the inverse solver. Several approaches are sucessfully used to solve this kind of problem, but usually they lack of a sound physical model at least in some parts of the domain.

Here we analyze the mathematical structure of a thermal problem with known interior conditions and then propose several numerical procedures to solve it. The proposed methods are tailored to the geophysical case and are based on the certainty of the different boundary conditions that are imposed in the model.

Moreover, because this thermal solver is expected to be used many times within an inverse scheme, we want the numerical mesh to be fixed. The LAB, therefore, will not fit the mesh.

How to cite: Zlotnik, S., Fernández, M. T., and Díez, P.: Structure and numerical solution of a thermal problem with internal Dirichlet conditions, EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-5381,, 2024.

On-site presentation
Stefan Markus Schmalholz, Mattia Luca Mazzucchelli, Lisa Eberhard, and Oliver Plümper

Dehydration reactions play a pivotal role in the dynamics and seismicity at subduction zones and in the deep water cycle. These reactions often occur during rock deformation. The dehydration of antigorite serpentinite is particularly important at subduction zones. This dehydration has been investigated with laboratory experiments of serpentinite deformation. Yet, the reproduction of such laboratory deformation experiments of serpentinite dehydration with mathematical models is still a major challenge. Here, we test a two-dimensional (2D) hydro-mechanical-chemical (HMC) numerical model for serpentinite dehydration by comparing the numerical results with the results of laboratory experiments.

The laboratory experiments are performed with a Griggs apparatus. Natural antigorite serpentinites with and without preferred orientation are deformed by vertical compression for a confining pressure of 1.5 GPa and maximum differential stresses between 350 and 700 MPa. For comparison, also experiments with hydrostatic stress are performed. The experimental temperature is between 620 and 650 °C. The applied confining pressure and temperature are in the olivine stability field according to the measured chemical composition of the serpentinite and thermodynamic Perple_X calculations. However, olivine only forms locally in the serpentinite if the serpentinite is deformed under differential stress. Olivine does not form in serpentinite under hydrostatic stress. Hence, we hypothesize that olivine formation is controlled by reaction kinetics and that the kinetics are locally faster in serpentinite that deforms under differential stress.

We elaborate a 2D HMC numerical algorithm that can simulate dehydration and olivine generation in a deforming serpentinite [Schmalholz et al., 2023]. The algorithm is based on a staggered finite difference discretization and employs a matrix-free, pseudo-transient iterative solver. Furthermore, the algorithm is programmed in the Julia language, employs the ParallelStencil package, and runs on GPUs. We discuss three major numerical challenges: First, the treatment of large changes in solid density during the generation of olivine by serpentinite dehydration. Second, the treatment of large temporal and spatial gradients in the unknowns, such as porosity and fluid pressure. Third, the treatment of strongly nonlinear relations between unknowns and parameters, such as the relations between density and fluid pressure, porosity and permeability, and porosity and rock viscosity. We implement several mathematical formulations for the reaction kinetics and discuss which formulation can explain the laboratory results best. We further discuss potential numerical benchmarks of such HMC algorithms for modelling the coupling of chemical reactions, fluid flow and rock deformation.       


Schmalholz, S. M., E. Moulas, L. Räss, and O. Müntener (2023), Serpentinite Dehydration and Olivine Vein Formation During Ductile Shearing: Insights From 2D Numerical Modeling on Porosity Generation, Density Variations, and Transient Weakening, Journal of Geophysical Research: Solid Earth, 128(11), e2023JB026985, doi:

How to cite: Schmalholz, S. M., Mazzucchelli, M. L., Eberhard, L., and Plümper, O.: Hydro-mechanical-chemical modelling of dehydration during serpentinite deformation: Comparing laboratory and numerical experiments, EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-6715,, 2024.

On-site presentation
Dániel Kiss, Viktoriya Yarushina, Boris Kaus, and Alexander Minakov

One of the continuing trends in geodynamics is to develop codes that are suitable to model thermo-hydromechanical processes with an increasing level of self-consistency. Developing such models is particularly challenging as some coupling terms can introduce strong non-linearities. We demonstrate that a finite difference discretization combined with pseudo transient solvers are well suited for such problems. Moreover, due to the inherent parallelism of the algorithm and thanks to the novel Julia language and the ParallelStencil.jl package, GPU implementation is not only feasible but also straightforward.

Here, we consider fluid flow in a deformable porous medium coupled to thermo-mechanical processes. We present a thermodynamically self-consistent set of governing equations, describing such processes. The governing equations consist of the conservation of mass, momentum, and energy in two phases. One phase represents the solid skeleton, which deforms in a poro-(visco-)elasto-plastic manner. The second phase represents a low viscosity fluid (water, CH, melt), percolating through the solid skeleton, that is described by Darcy’s law. A special process we will investigate is brittle failure of the matrix due to high fluid pressure (hydro-fracturing, fault reactivation, diking).

The system of equations is solved numerically, using the pseudo transient method, that is well suited to solve highly non-linear problems, as solving the global equations and iterating the non-linearities can be done at the same time. Moreover, the algorithm requires large number of local and cheap operations, which is ideal for GPU implementation. We will describe the governing equations, their numerical implementation, and show examples of numerical simulations that include mode-1 and mode-2 fractures. We apply the numerical model to study the effects of pore pressure in a siliciclastic reservoir on fault reactivation and associated integrity of cap rocks, and to provide a continuum analogue to magmatic dike emplacement.

How to cite: Kiss, D., Yarushina, V., Kaus, B., and Minakov, A.: Thermo-Hydromechanical modelling of poro-(visco-)elasto-plastic reservoir processes using GPUs in Julia    , EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-16070,, 2024.

On-site presentation
Yuan Li, Timothy Davis, Adina Pusok, and Richard Katz

Dykes are fluid-filled fractures that transport melt within the lithosphere. Their presence impacts lithospheric stresses and hence the forces of plate tectonics. Accurate representation of dykes in geodynamic models is crucial for modelling the dynamics of rifting. Li et al. (2023) approximated dyking as plastic tensile failure in a two-phase continuum with a poro-viscoelastic—viscoplastic (poro-VEVP) rheology. Li et al. only partially validated this model approach, and hence the extent to which it approximates a natural fracture remains unclear. In this study, we extend the comparison between our continuum formulation and the widely accepted theory of Linear Elastic Fracture Mechanics (LEFM) to the case of a buoyancy-driven dyke (Roper and Lister, 2007). We achieve this through detailed consideration of the dynamics and energetics.


Comparing the dynamics of the continuum and LEFM models is challenging due to their differing assumptions and the limitations of finite-difference numerical solutions. The continuum model treats the liquid phase as porous flow, while the LEFM model assumes a lubricated channel flow inside of an open fracture. A drawback of numerical computations is that the grid size can exceed the dyke width. In fact, one shortcoming of the previous computational framework is that the dyke width can grow to several grid cells, which further complicates the comparison. To overcome these challenges, we adjust both models to have the same geometry and mechanics: (1) we incorporate in the continuum model an anisotropic permeability treatment based on plastic strain components; (2) we introduce an intermediate 'poro-fracture' model that is a modified LEFM model for porous channel flow.


Our results show that the poro-VEVP model converges to the LEFM poro-fracture model at a toughness value that is large relative to natural rock and at a propagation speed that is slow relative to the standard LEFM formulation. Firstly, we attribute the slow propagation speed in the continuum model to the high Darcy’s drag force, which can be improved by augmenting the permeability. Secondly, we show that plasticity in the continuum model relates quantitatively to the toughness in the fracture model. It is evidenced by the good agreement between the two models in stress field prediction and also the dyke porosity profiles when a specific toughness value is applied to the fracture model. Intriguingly, this toughness is derived from equating the total energy dissipation rate in the continuum model to the energy required to open a fracture, thereby establishing an energy-based connection between the two models.


Li, Y., Pusok, A., Davis, T., May, D., and Katz, R., (2023). Continuum approximation of dyking with a theory for poro-viscoelastic–viscoplastic deformation, Geophysical Journal International.

Roper, S.M. and Lister, J.R., (2007). Buoyancy-driven crack propagation:the limit of large fracture toughness. Journal of Fluid Mechanics.

How to cite: Li, Y., Davis, T., Pusok, A., and Katz, R.: Continuum models of dykes: a comparison to discrete fracture models, EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-16931,, 2024.

Methodological advances
On-site presentation
Thomas Poulet and Pouria Behnoudfar

Traditional geomechanical modelling involves a lot of a priori knowledge regarding, among others, the distribution of material properties and boundary conditions to apply to a model. Setting those parameters is usually time-consuming, starting with extensive literature reviews to get initial estimates, which are then refined by trial and error by calibrating the numerical simulations against observations. In particular, the boundary conditions typically remain poorly constrained since precise data is rarely available.

In this contribution, we present a physics-based machine learning approach to infer the current displacements and full stress tensor distribution of an effective 2D linear elastic model, based on stress orientation and Global Navigation Satellite System (GNSS) data. This allows for automatic retrieval of a reasonable approximation of the model's elastic material properties, consistent displacement, and stress values over the whole physical domain, including at the boundaries, which could be used for instance in following forward simulations.

We show an application to the Australian continent, for which a rich dataset of stress orientation is available from the World Stress Map project and the GNSS measurements are particularly steady. This allows us to compare various options to account for stress orientation and displacement information as input data. Interestingly, we recover the smoothest stress field compatible with the (very accurate) GNSS observations and consequently identify areas where the resulting stress orientation differs from current estimates.

How to cite: Poulet, T. and Behnoudfar, P.: Retrieving stress and elastic properties distributions from stress orientation and satellite data using physics-informed machine learning, EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-1400,, 2024.

On-site presentation
Lijun Liu

Quantifying surface manifestations of deep mantle dynamics represents an ultimate goal of Earth science. This task has remained challenging due to uncertain initial and boundary conditions of the 4D nature of Earth evolution. A promising solution is through data assimilation, which substantiates a large number of model parameters with available data, thus greatly reducing model uncertainties. During the past decade, our group has developed both forward-in-time and backward-in-time data assimilation techniques that have greatly improved the quantitative expressions on the linkage between deep mantle dynamics and surface tectonics. Here, I will demonstrate how these data-assimilation models work in practice, and how oceanic subduction and the resulting mantle flow influence the Earth surface through generating 3D lithospheric deformation, intraplate volcanism, and earthquakes. Specific examples include the circum-Pacific plates subduction below continents, and evolution of the cratonic lithosphere.

How to cite: Liu, L.: Linking surface tectonics with mantle dynamics using numerical models with data assimilation, EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-16051,, 2024.

On-site presentation
Beatriz Martinez Montesinos, Yujiro J Suzuki, Leonardo Mingari, and Antonio Costa

Explosive volcanic eruptions can inject high quantities of magmatic materials into the atmosphere representing a risk to life and society. To quantify the potential impacts of a future eruption or to forecast what will happen in the next few hours when a volcano is erupting, atmospheric dispersion models are commonly used providing important information for civil protection and other stakeholders. Spatiotemporal distributions of volcanic ash in volcanic plumes are used as an input by a numerical simulation of tephra dispersal and are called eruption source parameters (ESPs). They have been poorly constrained and therefore their variation between models affects volcanic tephra hazard assessment. Since the goodness of ESPs increases with knowledge of the dynamics of eruptive columns, reconstruction of volcanic columns from past eruptions will improve the assessment of volcanic hazards for future eruptions.

In this work we take advantage of recent advances in computational capabilities and modeling in order to robustly reconstruct the dynamics of eruption columns from past eruptions. We do that by applying a synergistic approach between atmospheric dispersion models capable of reproducing the transport of volcanic ash due to atmospheric wind, eruption cloud dynamics models that resolve the ascending and the horizontal spreading of umbrella cloud, and inversion methods able to estimate ESPs using geological data information of tephra deposits.

Specifically, we use the latest version of the ash dispersal model FALL3D that allows us to determine EPSs by inverting field data using the novel GNC (Gaussian with non-negative constants) ensemble-based inversion method, and the eruption cloud dynamics model SK-3D that accurately resolves the turbulence of the volcanic plumes.

As an application, we focus on Campi Flegrei (CF) caldera, in Italy. CF is currently a densely populated area under busy air traffic routes where the monitoring system of the Vesuvius Observatory highlights some variations in the state of the volcanic activity. CF has generated several explosive eruptions in recent geological times, including the ~39 ka Campanian Ignimbrite (CI) super-eruption that is the largest explosive eruption in Europe in the last 200 ka. To reconstruct the CI super-eruption and assess such a huge eruption in CF, we apply our methodology to the geological data associated with this eruption.

How to cite: Martinez Montesinos, B., Suzuki, Y. J., Mingari, L., and Costa, A.: Synergistic Approach to Robustly Reconstruct Eruption Plume Dynamics: Application to Campi Flegrei, Italy, EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-10708,, 2024.

On-site presentation
Albert de Montserrat Navarro, Ludovic Räss, Boris Kaus, Ivan Utkin, and Pascal Aellig

Traditionally, the Earth science community has relied on statically compiled and explicitly typed, long-lived programming languages, namely C/C++ and Fortran, for computationally intensive production runs. In contrast, dynamic languages such as Python and MATLAB have served as alternatives for less computationally demanding tasks, prototyping, visualisation and teaching. The lower entry-level of the latter is often used as a “glue language” where the performance-critical code is written in a static language. In the recent decades other modern languages (C#, Nim, Go, Rust…) have been developed, but so far none of them have had a significant impact on the general scientific community.  Among the new modern programming languages there is Julia -with the 1.0 version being released just in 2018- which was designed with the computational scientific community as the main target user base. 
Julia’s main appeal for computational science are (i) performance, it produces compiled machine code with performance potentially similar to other statically compiled languages; (ii) interactivity,  it is a dynamic language, and most of the development and prototyping can be done in interactive sessions using its built-in REPL (read-eval-print loop); and (iii) readability, code is often easier to read than e.g. C++ or Rust. Julia also natively supports linear algebra operations (e.g. it ships its own wrappers for   libraries such as OpenBlas, SuiteSparse or LAPACK), multi-threading and distributed parallelism, which are crucial for a vast number of scientific applications. 

Additionally, Julia offers a series of features that are not intrinsic to the language itself which can ease and improve the user experience compared to traditional C/C++/Fortran: (i) Julia has a well-built package manager where one can easily add and install any version of any registered third-party package by typing two words in the REPL;  (ii) Julia itself is open source, as well as all the registered packages (hosted in Git repositories), and discourages black-box software; (iii) reproducibility is easy thanks to .toml files containing what exact version packages are needed to be installed.

Here we will discuss how the combination of the previous points, along with other important features of the language (multiple dispatch, metaprogramming,…), and existing scientific libraries (GPU support, auto-differentiation,…) make Julia an excellent platform for the geoscientific community to build an ecosystem of open-source  and highly modular software, opposite to the classic lone-wolf large monolithic code. Finally, we will present an overview of the current status of the JuliaGeodynamics and related organisations, as well as the general Earth sciences ecosystem.

How to cite: de Montserrat Navarro, A., Räss, L., Kaus, B., Utkin, I., and Aellig, P.: Building an ecosystem of computational geosciences software: why Julia?, EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-16244,, 2024.


Posters on site: Wed, 17 Apr, 10:45–12:30 | Hall X1

Display time: Wed, 17 Apr 08:30–Wed, 17 Apr 12:30
Chairpersons: Boris Kaus, Ivan Utkin, Dave May
Samuel Omlin, Ludovic Räss, and Ivan Utkin

The ongoing digitalisation in geosciences, along with higher resolution data and more powerful computers, necessitates efficient workflows and tools for processing large amounts of data, modelling processes, and addressing new problems and challenges.

The Julia language offers an great basis for this by combining the benefits of a high-level language, such as ease of use and interactivity, with the features of a low-level language, including speed, efficiency, scalability, and native GPU support.

We will present recent efforts to accelerate geocomputing using Julia at scale. These efforts are supported by the PASC-funded GPU4GEO project, in collaboration with the Swiss National Supercomputing Centre. We will showcase the current HPC building blocks that we have developed, which enables geoscientists to write high-performance stencil codes that can scale from their laptops to the largest supercomputers. Additionally, we will demonstrate preliminary results on using automatic differentiation to perform inverse modelling. This allows us to efficiently constrain models with data.

Additionally, we discuss how these recent developments, along with future ones, will accelerate geocomputing and enable the education of the next generation of geoscientists.

How to cite: Omlin, S., Räss, L., and Utkin, I.: Accelerating Geocomputing with Julia at Scale, EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-19607,, 2024.

Mathis Bergogne, Laetitia Le Pourhiet, Ludovic Räss, Yury Podladchikov, Ivan Utkin, and Alexis Plunder

Pegmatites (and rare metal granites) are igneous rocks with a granitic composition, characterised by crystal growth dominated texture. They are often enriched in rare elements (such as Li, Cs, Be, Nb, Ta…) and offer valuable deposit of economic interest that belong to the list of critical raw material defined by the European commission. The aim is to use modeling to understand the formation of pegmatite, and especially the parameters that control the migration distance between their sources (granite, migmatites) and their level of emplacement.

We use a two phase flow finite difference code, in Julia, based on the porosity waves with compressible fluid in compressible medium, where the porosity is interpreted as melt [1,2], to model the magma migration inside migmatitics domes. To improve the yet existing codes, we implement temperature in our two phase flow formulation, it will be calculated from energy conservation to take into account the latent heat of the partial melting and the crystallisation. It will allow the to stop the experiments upon magma crystallisation. It will also be use to increase the viscosity of the melt while it cools down.

We here present the results of a preliminary 1D version of our model (without temperature). It shows that an increase in the ratio of matrix permeability over the fluid viscosity results in a greater distance travelled by the melt, for a constant number of time step. For a fluid viscosity of 104 Pa.s the increase of matrix (dynamic) permeability from 10-13 to 10-11 m-2 fasten the migration of the melt of a factor 2 and the travelled distance by a factor 1. When the energy will be implemented, we will compare the impact of matrix permeability, the fluid viscosity and the geothermal gradient on the height of the magmatic migration.



[1] L. Räss, T. Duretz & Y.Y. Podladchikov (2020).

[2] A. Plunder & al. (2022).

How to cite: Bergogne, M., Le Pourhiet, L., Räss, L., Podladchikov, Y., Utkin, I., and Plunder, A.: 3D Modelling of pegmatites migration at the onset of partial melting, EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-19212,, 2024.

Manuele Faccenda, Brandon VanderBeek, Albert de Montserrat, and Jianfeng Yang

In this contribution we introduce ECOMAN 2.0, an open-source software package for (1) modelling the strain-/stress-induced rock fabrics and related mechanical anisotropy, and (2) performing isotropic and anisotropic inversions using real/synthetic P- and S-wave travel-times and S-wave splitting parameters.  

The strain-induced intrinsic mantle fabrics are modelled inputting the velocity, pressure, temperature and dominant creep mechanism fields from large-scale mantle flow simulations into D-Rex (Kaminski et al., 2004). This open-source software has been parallelized using a hybrid MPI and OpenMP scheme and modified to account for combined diffusion-dislocation creep mechanisms, LPO of transition zone and lower mantle polycrystalline aggregates (Wadsleyite, Bridgmanite, post-Perovskite), P-T dependence of single crystal elastic tensors, advection and non-steady-state deformation of crystal aggregates in 2D/3D cartesian/spherical grids (Faccenda, 2014; Faccenda and Capitanio, 2013). The new version of D-Rex can solve for the LPO evolution of 100.000s polycrystalline aggregates of the whole mantle in a few hours, outputting the full elastic tensor of poly-crystalline aggregates as a function of each single crystal orientation, volume fraction and elastic moduli scaled by the local P-T conditions.

Extrinsic elastic anisotropy due to grain- or rock-scale fabrics or fluid-filled cracks can also be estimated with the Differential Effective Medium (DEM) (Faccenda et al., 2019). Similarly, extrinsic viscous anisotropy can be modelled yielding viscous tensors to be used in large-scale mantle flow simulations (de Montserrat et al., 2021). 

The elastic tensors can then be interpolated in a tomographic grid for (i) visual inspection of the mantle elastic properties (such as Vp and Vs isotropic anomalies; radial, azimuthal, Vp and Vs anisotropies), (ii) generating input files for large-scale synthetic waveform modelling (e.g., SPECFEM3D format), or (iii) P- and S-wave isotropic and anisotropic inversions (e.g., Faccenda and VanderBeek, 2023). The latter can be performed with the new PSI (Platform for Seismic Imaging) module, which includes recently developed techniques for seismic anisotropic inversions of body waves (VanderBeek and Faccenda, 2021; VanderBeek et al., 2023).

How to cite: Faccenda, M., VanderBeek, B., de Montserrat, A., and Yang, J.: ECOMAN 2.0: an open-source software package for coupling geodynamic and seismological modelling., EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-3547,, 2024.

Eugenio D'Ascoli, Andreas Burkhart, Nils Kohl, Hans-Peter Bunge, and Marcus Mohr

Simulating the Earth’s mantle convection at full convective vigor on planetary scales is a fundamental challenge in Geodynamics even for state of the art high-performance computing (HPC) systems. Realistic Earth mantle convection simulations can contribute a decisive link between uncertain input parameters, such as the mantle viscosity structure, and testable preconditions, such as dynamic topography. The vertical deflections predicted by such models may then be tested against history of dynamic topography from stratigraphic observations. Considering realistic Earth like Rayleigh numbers (∼ 108 ) a resolution of the thermal boundary layer of 10 − 50 km is necessary considering the volume of the Earth’s mantle.

Simulating Earth’s mantle convection at this level of resolution requires solving sparse indefinite systems with more than 1012 degrees of freedom, computationally feasible only on exascale HPC systems. This is achievable only by mantle convection codes providing high degrees of parallelism and scalability. Earlier approaches with prototype frameworks using hierarchical hybrid grids (HHG) as solvers for such systems demonstrated the scalability of the underlying concept for future generations of exascale computing systems.

Building up on the TerraNeo project here we report on the progress of utilizing the improved framework HyTeG (Hybrid Tetrahedral Grids) based on matrix-free multigrid solvers in combination with highly efficient parallelisation and scalability. This will allow to solve systems with more than a trillion degrees of freedom on present and future generations of exascale computing systems. We also report on the advances in developing the scalable mantle convection code TerraNeo using the HyTeG framework to realise extreme-scale mantle convection simulations with realistic, Earth like parametrisation and a resolution in the order of ∼ 1km.

How to cite: D'Ascoli, E., Burkhart, A., Kohl, N., Bunge, H.-P., and Mohr, M.: TerraNeo: Development of a scalable mantle convection code for exascale computing, EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-6508,, 2024.

Thibault Duretz and Laetitia Le Pourhiet
Frictional plasticity governs the state of stress and the deformation style of the Earth’s upper crust. Frictional models depend on both deviatoric stress and pressure and include coupling between shear and volumetric stress via dilatancy and compaction. Elastic-plastic strains can takes the form of shear bands that can further be compared to geological observations. The study of frictional shear banding is non trivial because (1) this process is highly non-linear and (2) the models may require enrichment/regularisation to yield successful integration. 2D or 3D models may be used to model shear banding, they are however computationally expensive, which hinders systematic exploration of shear banding processes. We here present a 1D dimensional model of frictional plastic strain deformation. These models allow for exploring the process of strain localisation and the coupled effects of fluid flow and material microstructure.  

How to cite: Duretz, T. and Le Pourhiet, L.: One dimensional models of frictional plastic strain localisation, EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-8080,, 2024.

Jon Bryan May, Peter Bird, and Michele Matteo Cosimo Carafa

Geodynamic forward models are typically computationally expensive, they generally draw on a large set of input and testing data and will normally have quite large model domains in order to ensure that the results are practical and applicable. To reduce simulation time many programs will use some form of parallel computing in their calculations, e.g., programs may use the message passing interface (MPI) communication protocol or OpenMP API to partition the domain into reasonable chunks divided between multiple processes, this can improve performance by having multiple processes solve smaller parts of a large problem in parallel.

Geodynamic inverse modelling involves using known phenomena to constrain a model, or set of models, to improve knowledge of unknowable or unmeasurable parameters involved within a dynamic system. For example, one could use measured GPS velocities and seafloor spreading rates to infer changes to mantle convection or combine GPS velocities with near-surface temperatures to monitor the growth of magma chambers.

Since geodynamic inversion involves running multiple constrained forward models it naturally suffers from the performance issues linked to the forward models themselves, briefly: a poorly performing forward model will mean a poorly performing inverse model. More than that, the inverse model must perform multiple (the more the better) forward models with updated parameters. These requirements, a good enough forward model, and an efficient method of performing multiple models in parallel while optimizing the desired parameters, leads to an obvious conclusion: to perform forward models in parallel while searching for an optimal model.

To this end we present a geodynamic inversion model, which we call ShellSetHPC, which uses a combination of existing, well known, and robust software to model the neotectonics of planetary lithosphere. This is further combined with an efficient, de-centralised random search algorithm able to generate testable models within a user defined N-dimensional space. This algorithm is also able to launch multiple models in parallel thanks to an MPI framework. The forward model makes use of Intel’s thread safe math kernel library (MKL) to solve the linear system in parallel. This hybrid approach lends itself to use on high performance computing (HPC) machines which would allow a more complete utilisation of these features.

In this work we will present scaling results on an HPC cluster and compare these with results obtained from more typical search algorithms. All tests are performed within a realistic geological setting, the results of which will be used to gather insight into the performance of the driving model generators and their settings.

How to cite: May, J. B., Bird, P., and Carafa, M. M. C.: ShellSetHPC – Parallel dynamic neotectonic modelling, EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-11777,, 2024.

Eric Hutton, Michael Steckler, and Gregory Tucker

Landlab is an open-source Python package that streamlines the creation, combination, and reuse of 2D numerical models and is a key element of the Community Surface Dynamics Modeling System (CSDMS) Workbench. The Landlab Toolkit provides building blocks for model development such as grid data structures, input/output functions, and a library of several dozen components that each model a separate physical process. Additionally, it provides a framework for assembling integrated models from component parts. We've found that Landlab significantly accelerates model development, encourages user-developers to adopt standard practices and contribute new components to the library. It serves as a platform that nurtures a community of model developers, assisting them in creating coupled models to investigate non-linear interactions between geologic processes.

Using the Landlab toolkit, we developed a new model, Sequence, which is a modular 2D (i.e., profile) sequence stratigraphic model that incorporates key geophysical processes influencing accommodation space in both terrestrial and marine environments. These factors include tectonics and faulting, eustatic sea level changes, flexural isostatic compensation of sediment and water, sediment compaction, and hypopycnal sediment plumes. Each process is encapsulated as an individual, standalone Landlab component, providing flexibility in the construction of new models. Sequence serves not only as a distinct model, but also as a scaffold for the development of new models.

Sequence simulates the evolution of stratigraphy on a continental margin over time scales ranging from thousands to millions of years. Sediment transport and deposition primarily occur during infrequent, high-energy events like storms and floods. For these extended time frames, Sequence employs a scale-integral approach. This method utilizes differential equations to summarize the cumulative effect of sediment transport and deposition across different depositional environments over longer periods (e.g., on the order of a hundred years). The model features a moving-boundary formulation to track shoreline changes and partitions the domain into distinct areas: coastal plain, continental shelf, and upper and lower slope/rise. Submarine sediment transport and deposition are modeled through nonlinear diffusion, with a diffusion coefficient that varies inversely with water depth. The model tracks evolving stratigraphic layers and sediment lithology that is a mixtuer of two grain sizes (sand and mud) each with separate transport functions.

How to cite: Hutton, E., Steckler, M., and Tucker, G.: Sequence: A coupled sequence-stratigraphic model built using Landlab, EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-14445,, 2024.

Zenan Huo, Michel Jaboyedoff, Yury Podladchikov, Ludovic Räss, and Emmanuel Wyser

We present a high-performance Material Point Method (MPM) numerical solver for simulating the landslide run-out process with high resolution. The solver has been adapted to run on multi-GPU platforms. The current version is backend-agnostic, operating efficiently across various CPU and GPU hardware from different vendors, utilizing the same codebase. We evaluate multiple performance metrics and ensure minimal data synchronization between different devices at each iteration. We validate the solver's accuracy by comparing the simulation results of an aluminum-bar collapse with the corresponding experimental outcomes. Consistency is observed between numerical and experimental results for the free and failure surfaces. The results also indicate favorable performance scalability in high-resolution models, significantly enhancing computational efficiency. Further, we simulate the slumping mechanic problem with a simplified landslide geometry. The results show that the shear band develops within the high plastic strain area. The failure surface is in good agreement with the solution reported by Huang[1] et al., demonstrating that MPM can accurately handle failure and large deformation problems as they occur in landslides. Our multi-GPU implementation using MPI makes it possible to perform large-scale simulations that enable to tackle research in the field of geotechnical engineering.

[1]. Huang, Peng, Shun-li Li, Hu Guo, and Zhi-ming Hao. “Large Deformation Failure Analysis of the Soil Slope Based on the Material Point Method.” Computational Geosciences 19, no. 4 (August 2015): 951–63.

How to cite: Huo, Z., Jaboyedoff, M., Podladchikov, Y., Räss, L., and Wyser, E.: Multi-GPU Material Point Method Solver for Landslide Simulation, EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-15096,, 2024.

Marcel Thielmann, Boris J.P. Kaus, Arne Spang, Christian Schuler, and Pascal Aellig

Geophysical datasets and their interpretations form the basis of geodynamic simulations of the Earth’s mantle and lithosphere. Yet, creating geodynamic models from these datasets is often non-trivial, particularly in complex regions such as active orogens. Creating self-consistent three-dimensional models from these datasets poses a challenge due to technical issues such as different data set formats and different spatial resolutions, but also due to discrepancies in the data itself. At the same time, the different datasets obtained through initiatives such as AlpArray contain a wealth of data that can help to constrain subsurface models to an unprecedented extent. Yet interpreting these different data still involves subjective steps and ideally different datasets are combined in the process.

To facilitate the joint interpretation of these datasets and the generation of geodynamic model setups, we therefore developed an open-source package - the Geophysical Model Generator (GMG) - to assist with unifying these datasets in a common data format that can then be further used to visualize, compare and interpret data. Within this package, we provide a set of routines to import different datasets, convert them to a common data format and to process them further (e.g., to create vote maps from different tomographies). These unified datasets can then be exported as vtk-files for further 3D visualization (e.g., Paraview). Moreover, with the Geophysical Model Generator it is also possible to create model setups for numerical models (such as the 3D geodynamic code LaMEM). This package thus covers the entire workflow from data import to numerical model generation. Key features of the Geophysical Model Generator include 1) the creation of 3D volumes from seismic tomography models, 2) the import of 2D data (e.g., surface or Moho topography or screenshots from published papers) and 3) the incorporation of point data such as earthquake locations or GPS measurements. Both scalar and vector data can be handled. With these tools, one can then create a consistent overview of the entire data available for a given region.

The package is written in Julia and hosted as a public open-source repository on GitHub ( To assist the joint interpretation of different geophysical datasets, we furthermore provide a graphical user interface that allows to view and compare them ( The GUI provides an interactive interface, allows loading different datasets and facilitates the manual interpretation of different structures (such as subducting slabs) along profiles and visualize them in 3D while taking different data into account. 

How to cite: Thielmann, M., Kaus, B. J. P., Spang, A., Schuler, C., and Aellig, P.: The Geophysical Model Generator: a tool to unify and interpret geophysical datasets, EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-15998,, 2024.

Boris Kaus, Albert de Montserrat, Nicolas Riel, and Pascal Aellig

Julia is a relatively new scientific computing language with which you can write code that is nearly as fast as Fortran or C, but being a dynamic language, it remains very compact. Since the language is newer than say Python, many of the shortcomings of earlier languages have been fixed. Installing and updating packages with external dependencies is straightforward, for example, and exact reproducibility of computational results can be done by uploading a single additional file.

Here we will discuss our recent experiences with using Julia for teaching and research:

  • The binarybuilder package simplifies the precompilation of non-Julia open-source packages (C, C++, Fortran), allowing cross-compilation for Apple, Linux, and Windows. Uploaded centrally, these precompiled binaries can be effortlessly installed via the Julia package manager. This expedites the setup of complex research code in teaching, reducing the time spent on package installation during the initial lecture to just a matter of minutes. Examples where this was recently implemented was GMT.jl (the Julia interface to the Generic Mapping Tools) which now comes with precompiled version of GMT. Likewise, we provide precompiled packages for LaMEM (a parallel 3D geodynamic code) and MAGEMin (a new thermodynamic code).
  • It is quite easy to call such binary packages, either by calling the executable itself or by automatically creating wrapper (using Clang.jl) and calling the functions within dynamic libraries directly from julia.
  • There is a large and ever-growing ecosystem of existing packages that are useful for computational geosciences. Whereas it is still not as mature as in Python or MATLAB, much of what is required in our daily life is already there, including plotting, machine learning, data I/O.
  • It is easy to develop additional packages, and testing is build-in the language (and CI/CD is free on GitHub for open-source packages). That encourages users to implement tests and helps maintaining packages.
  • Packages can talk to each other in a straightforward manner, which facilitates building composable software stacks.
  • There is great support for GPU computing.
  • It is possible to create GUI’s that run in a webbrowser (for example by using Dash.jl). Recent examples where we used that is in InteractiveGeodynamics.jl (which provides user interfaces to simulate typical geodynamic problems such as convection), GeoDataPicker.jl (a 2D/3D tool to analyze and visualize geodynamic data) and MAGEMin_app (a web-based tool to compute phase diagrams).
  • It can directly be used in Jupyter notebooks (as the “Ju” in its name stands for Julia).
  • It is easy for students to start writing simple codes without having to think about definition of types or having to load packages to create vectors or matrixes, which makes it well-suited for teaching quantitative classes. In many cases, these simple codes already run very fast, as the Julia compiler does an excellent job in optimisations. With a few additional tricks (ensuring type stability and reducing allocations) the simple examples can be turned into a high-performance code that runs at nearly the speed of the more classical languages, while having the advantage to still be very compact and readable.

We will highlight these aspects including with interactive demos.

How to cite: Kaus, B., de Montserrat, A., Riel, N., and Aellig, P.: Using Julia for teaching and research in the Geosciences, EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-16725,, 2024.

Timothy Gray, Taras Gerya, and Paul Tackley

The study of coupled Earth systems, and in particular the coupled interactions between the lithosphere, atmosphere, and biosphere, have received greater attention in recent years (Gerya et al. 2020). Interactions between these systems occur primarily at the surface, and are driven on the large scale by topographic and bathymetric evolution controlled by deep mantle processes. However, due to the large difference in length scales between the mantle and the surface, it is difficult to capture topographic evolution to a high degree of accuracy in existing global mantle convection models including a free surface boundary condition.

Global mantle convection models incorporating a free surface often employ a marker-in-cell technique with a layer of “sticky air” (i.e. material with the density of the air and sufficiently low viscosity, which is still much higher than that of real air) to characterise the surface. However, accurate topographic evolution using this method requires a high density of markers near the surface. This need for additional computational resources motivates alternative methods of tracking the interface between the air and rock layers, as is done frequently in existing multiphase fluid flow codes. We demonstrate the implementation of two such methods of modelling the surface. The first is a Lagrangian marker chain (or mesh in 3D models) which, when combined with an appropriate remeshing procedure, directly tracks the rock-air interface. The second is a volume of fluid approach adapted from the open source code gVOF using the unsplit volume of fluid library gVOF (López & Hernández, 2022).

We demonstrate toy models and benchmarks (based on those in Crameri, 2012) comparing the Lagrangian marker method and the volume of fluid methods as implemented in the global scale mantle convection code StagYY (Tackley, 2008). Models of global scale topography and evolution produced using StagYY may then be used as a tool for further studies on the coupling of mantle dynamics with modelling of the landscape, and the evolution of the atmosphere and biosphere. Initial applications include modelling hypsometric curves from global scale numerical models, and the tracking of sea level changes over time.

How to cite: Gray, T., Gerya, T., and Tackley, P.: Free surface methods applied to global scale numerical geodynamic models, EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-18423,, 2024.

Peter Mora, Gabriele Morra, Leila Honarbakhsh, Christian Huttig, and Nicola Tosi

The Thermal Lattice Boltzmann Method (TLBM) for geodynamical simulation research offers an alternative to classical PDE based methods for 2D and 3D geodynamics simulation research. It is based on modelling the Boltzmann equation on a discrete lattice which involves the movement of number densities carrying mass and energy density on a discrete lattice and their relaxation to equilibrium which model collisions. We present examples in 2D and 3D to illustrate the capabilities, performance, and accuracy of this method for geodynamics research, namely: (1) ability to handle highly nonlinear rheology, ultra-high Rayleigh numbers, a wide range of Prandtl numbers, and multiphase flow, (2) linear scaling up to 300K cores on HPC CPU clusters, and (3) ability to closely match the Blankenbach benchmarks demonstrating the LBMs accuracy. Examples in 2D include high Rayleigh number simulations to Ra = 1015, highly nonlinear rheology leading to the emergence of plate-tectonic like behaviour, and planetary accretion. Examples in 3D include modelling of a mantle with an aspect ratio of 25x25x1 representing a case from a recent nature paper, and modelling a case of an aspect ratio of 14.4x14.4x1 which is similar that of the Earth for Ra = 106 and Pr = 100. Potential benefits of the TLBM include an ability for higher resolution simulations than can be achieved using classical methods, and faster simulations which may allow phase space studies to determine which parameter combinations lead to which class of behaviour. As the TLBM is a new method for geodynamical simulation, it will take some time to determine the limits of this method. For example, a simulation can be made to run faster by increasing the physical time step, but eventually, if the time step is too large, the Mach numbers on the lattice become too high leading to lower accuracy and eventually instability due to non-convergence of the collision step which involves a relaxation of the number densities to equilibrium. We believe that over time, these limitations will become well understood and that the outstanding parallel scaling performance on HPC CPU clusters of the TLBM - which makes possible 3D models up to 50003 - will  open up exascale computing to geodynamics research and will lead to fundamental advances in geodynamics research. As such, the TLBM may become a valuable tool to advance geodynamics research into the future through large to exascale simulations that may lead to new insights into the dynamics and evolution of the earth and exoplanets from the early lava world stage through to plate tectonics or other regimes.

How to cite: Mora, P., Morra, G., Honarbakhsh, L., Huttig, C., and Tosi, N.: Progress and perspectives on using the Lattice Boltzmann Method for geodynamics simulation research, EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-14391,, 2024.

Anton A. Popov, Nicolas Berlie, and Boris J.P. Kaus

Developing a deep understanding of magmatic processes, such as the ascent of magmatic melts through the lithosphere, is a notoriously complex interdisciplinary task that involves contributions from various branches of geosciences. Here, we focus on the implementation of mode-1 plasticity, which is highly relevant for the modeling of dyke propagation through the brittle crust under an excess of fluid (melt) pressure, as part of a nonlinear elasto-visco-plastic rheology that is also appropriate for modeling the ductile-brittle transition in the upper mantle.

For this, we adopt a coupled poro-elasto-visco-plastic description of the strain localization process to capture the onset and advance of vertical dykes originating from magma reservoirs located at various depth. We mostly aim to investigate the hydro-mechanical interactions of such a system (e.g. permeability increase and elastic modulii degradation with increasing plastic strain, and influence of fluid pressure on both total stresses and on the yield surface). Thermal effects are essentially ignored, apart from the background geothermal gradient. The reservoir is represented as a cavity subjected to an increased fluid pressure.

We describe the brittle deformation using multi-surface plasticity models in the framework of the flow plasticity theory. A number of challenging problems are usually encountered in this case from an algorithmic viewpoint, including a lack of convexity and continuity of both the yield surface and flow potential, spurious elastic domains, singularity points and loss of convergence. We address these issues by using a relatively simple Perzyna-type visco-plastic model that consists of a linear Drucker-Prager shear failure envelop, and a circular tensile cap. Both surfaces are combined with each other in a way that enforces dimensional consistency, convexity, and continuity throughout the entire stress space.

Finally, we discuss the algorithmic details of the model which is incorporated in an implicit finite element code and demonstrate the  application results.

How to cite: Popov, A. A., Berlie, N., and Kaus, B. J. P.: Modeling fluid-induced dyke propagation in elasto-visco-plastic rocks, EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-12283,, 2024.