Enter Zoom Meeting


Advances in Numerical Modelling of Geological Processes: Methods and Applications

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/TS9
Convener: Ludovic RässECSECS | Co-conveners: Boris Kaus, Thibault Duretz, Dave May
| Thu, 26 May, 15:55–18:30 (CEST)
Room -2.47/48

Thu, 26 May, 15:10–16:40

Chairpersons: Boris Kaus, Thibault Duretz



On-site presentation
Magali Billen et al.

Subduction is driven by difference in mass between the sinking plate and the surrounding mantle. The deformation calculated in numerical models of subduction is strongly dependent on the magnitude of the mass difference. The mass difference depends on the temperature of the slab. As the tectonic plate sinks it heats up, but it also cools down the surrounding mantle. The amount of heating and cooling is determined by conservation of thermal energy. Because the temperature also determines the thermal mass, conversing thermal energy also leads to conserving mass. For some studies, models of subduction are made to match the present day structure of a sinking plate. In this case, the temperature is defined to follow the observed geometry. In some previous studies, the temperature structure did not explicitly enforce conservation of energy or mass, and thus the density of the slab was not physically consistent, which is added source of uncertainty when analyzing the resulting flow and sensitivity of model results to mantle and slab rheology. Here we present a mass-conserving thermal structure for slabs that also creates a smoothly varying temperature structure. The thermal structure is based on a 1-D half-space cooling model (bottom) and an infinite space cooling model (top). It uses the age of the plate at the trench to determine the initial mass anomaly of the slab. The sinking velocity modifies the rate of heating and migration of the minimum temperature into the slab interior. The thermal model is calibrated against simple 2D subduction models in which the age and subduction velocity are held fixed. The new thermal structure has been implemented in the Geodynamic WorldBuilder (1), which can be used with different mantle convection software and is distributed as a plugin for ASPECT (2). Comparison of model results with the mass conserving slab thermal structure to the "plate" model from McKenzie (1970) is used to illustrate the differences in modeled results. References: 1. Fraters, M. R. T. et al., Solid earth, 2019. 2. Bangerth, W. et al., https://doi.org/10.5281/ZENODO.5131909, 2021.

How to cite: Billen, M., Fraters, M., and Babin, M.: Mass-Conserving Thermal Structure for Slabs in Instantaneous Models of Subduction, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-8849, https://doi.org/10.5194/egusphere-egu22-8849, 2022.

Virtual presentation
Emilie Macherel et al.

Power-law viscous flow describes 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. Calculating these stresses in a three-dimensional (3-D), geometrically and mechanically heterogeneous lithosphere requires high-resolution and high-performance computing.


Here, we present numerical simulations of 3-D power-law viscous flow. We employ the pseudo-transient finite difference (PTFD) method, which enables efficient simulations of high-resolution 3-D deformation processes by implementing an iterative implicit solution strategy of the governing equations. The main challenges for the PTFD method are to guarantee convergence, minimize the required iteration count and speed-up the iterations. We implemented the PTFD algorithm using the Julia language (julialang.org) to enable optimal parallel execution on multiple CPUs and GPUs using the ParallelStencil.jl module (https://github.com/omlins/ParallelStencil.jl). ParallelStencil.jl enables execution on multi-threaded CPUs and Nvidia GPUs using a single switch.


We present PTFD simulations of mechanically heterogeneous (weak and less dense spherical inclusion), incompressible 3-D power-law viscous flow under gravity in cartesian, cylindrical and spherical coordinates systems. The viscous flow is described by a linear combination of a linear viscous and a power-law viscous flow law, representing diffusion and dislocation creep, respectively. The iterative solution strategy builds upon pseudo-viscoelastic behavior to minimize the iteration count by exploiting the fundamental characteristics of viscoelastic wave propagation. We performed systematic numerical simulations to investigate the impact of (i) buoyancy versus shear forces and (ii) linear versus power-law viscous flow on the vertical velocity of the spherical inclusion under bulk strike-slip shearing. We report the systematic results using the controlling dimensionless numbers and compare the numerical results with analytical predictions for buoyancy-driven flow of inclusions in a power-law matrix. We also aim to unveil preliminary results for a vertically and locally loaded power-law viscous lithosphere showing the impact of different lithosphere curvatures on the resulting stress field.

How to cite: Macherel, E., Podladchikov, Y., Räss, L., and Schmalholz, S. M.: GPU-based pseudo-transient finite difference solution for 3-D gravity- and shear-driven power-law viscous flow, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-8816, https://doi.org/10.5194/egusphere-egu22-8816, 2022.

On-site presentation
Yuan Li et al.

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 (Keller et al., 2013). Here we present our progress in consistent numerical modelling of poro-viscoelastic-viscoplastic (VEVP) flow. We show that a VEVP model with a new, hyperbolic yield surface can help to robustly simulate both shear and tensile modes of plastic failure in a two-phase system. 

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 with shear failure only. In two-phase systems, tensile failure of rocks can also occur due to an overpressured liquid phase. 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.

An appropriate rheological model is required to meet this challenge. The most straightforward choice is a Maxwell visco-elasto-plastic model, but this leads to grid-scale localisation and hence mesh-dependence. To obtain mesh-independent shear localisation, we employ the visco-elasto-viscoplastic model by introducing a viscous dashpot in parallel to the plasticity element. Whilst this formulation has shown promise in regularising shear failures in a single-phase flow model (de Borst and Duretz, 2020), its incorporation within two-phase systems has not been examined. We will show that the linear Griffith criteria for the tensile failure can lead to convergence issues whereas a new, hyperbolic yield surface is proposed to resolve these numerical issues. This yield surface provides a smooth transition between the two modes of failure.

The underlying PDEs are discretised using a conservative, finite-difference, staggered-grid framework implemented with PETSc (FD-PDE) that supports single-/two-phase flow magma dynamics. Here, we present simplified model problems using the FD-PDE framework for poro-viscoelastic-viscoplastic models designed to characterise the solution quality and assess both the discretisation and solver robustness. It has been observed that employing the hyperbolic yield surface improved the robustness in simulating plastic failures in both modes.



Keller, T., May, D. A., & Kaus, B. J. P., (2013). Numerical modelling of magma dynamics coupled to tectonic deformation of lithosphere and crust, Geophysical Journal International, v195, 1406-1442, https://doi.org/10.1093/gji/ggt306.

de Borst, R., Duretz, T., (2020). On viscoplastic regularisation of strain-softening rocks and soils. International Journal for Numerical and Analytical Methods in Geomechanics, v44, 890-903. https://doi.org/10.1002/nag.3046.

How to cite: Li, Y., Pusok, A., May, D., and Katz, R.: Simulation of partially molten rocks with visco-elasto-viscoplastic rheology and a hyperbolic yield surface for plasticity, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-5594, https://doi.org/10.5194/egusphere-egu22-5594, 2022.

On-site presentation
Lukas Fuchs et al.

The formation and maintenance of narrow, lithospheric shear zones and their importance in plate-tectonics remain one of the major problems in geodynamics. While the cause and consequence of strain localization and weakening within the lithosphere remain debated, it is clear that these processes play an essential role in lithospheric deformation across a wide range of spatio-temporal scales. Here, we analyze the efficiency of strain localization in a 2-D visco-elasto-plastic medium for a strain-dependent weakening and healing (SDWH) rheology using 2-D numerical, thermo-mechanical experiments with kinematic boundary conditions. Such a parameterized rheology successfully mimics more complex transient weakening and healing processes, akin to a grain-size sensitive composite (diffusion and dislocation creep) rheology. In addition, the SDWH rheology allows for memory of deformation. This enables self-consistent formation and reactivation of inherited weak zones within the lithosphere and sustains those weak zones over an extended period of time. We further analyze the resulting shear zone patterns and seek to answer the questions: What is the typical, effective intensity of strain localization? What are the dimensions of the resulting shear zones? Are such shear zones mesh-dependent in numerical models and, if so, can we exploit existing regularization approaches for the SDWH rheology?

How to cite: Fuchs, L., Duretz, T., and Becker, T. W.: Strain localization in a visco-elasto-plastic medium using strain-dependent weakening and healing rheology, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-11469, https://doi.org/10.5194/egusphere-egu22-11469, 2022.

On-site presentation
Nicolas Berlie et al.

One of the great challenges involved in modelling the lithosphere is its plastic behaviour, especially when dealing with compressible materials. Shear fractures are designated as mode 2 and 3 and can be described using a Linear Mohr Coulomb envelope  or a simplification of it like Drucker-Prager. Meanwhile, mode 1 fractures are created when the normal stresses become tensile  and require another yield function, such as the Griffith criterion or a tension cap function.

While the governing equations are well known and widely employed in engineering codes, they are usually expressed with a displacement formulation. Most geodynamic codes, on the other hand, use pressure and velocity as their primary variables. A numerically robust method that takes all plasticity modes into account in a staggered finite difference discretization remains an open task. Here we present a composite yield function implemented with pressure-velocity formulation, capable of producing produce shear and tensile failure.

We have implemented this in a new code that employed PETSc through the recently updated PETSc.jl Julia interface, while utilizing the automatic differentiation tools in julia. We found this workflow to significantly reduce the development time of complex nonlinear coupled codes.  

We will describe the implementation, propose regularization schemes and discuss benchmark cases and simple applications. We demonstrate Newton convergence for most cases and will discuss different methods to combine multiple plastic flow laws.

How to cite: Berlie, N., Kaus, B., Popov, A., Kiss, D., and Riel, N.: How to break the lithosphere: a compressible pressure-velocity formulation for elasto-visco-plastic rheologies that includes shear and tensile failure with dilation, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-6566, https://doi.org/10.5194/egusphere-egu22-6566, 2022.

Thu, 26 May, 17:00–18:30

Chairpersons: Ludovic Räss, Dave May


On-site presentation
Casper Pranger et al.

The rate- and state-dependent friction (RSF) laws (Dieterich, 1979; Ruina, 1983) have been widely successful in capturing the behavior of sliding surfaces in laboratory settings, as well as reproducing a range of natural fault slip phenomena in numerical models.

Studies of exhumed fault zones make it clear that faults are not two-dimensional features at the junction of two distinct bodies of rock, but instead evolve into complex damage zones that show clear signs of multi-scale fracturing, grain diminution, hydro-thermal effects and chemical and petrological changes. Many of these observed factors have been experimentally verified, and several studies have furthered our theoretical understanding of earthquakes and other seismic phenomena as volumetric, bulk-rock processes, including Sleep (1995, 1997), Lyakhovsky and Ben-Zion et al. (2011, 2014a,b, 2016), Niemeijer and Spiers et al. (2007, 2016, 2018), Roubicek (2014), and Barbot (2019).

While the established numerical modeling approach of simulating faults as planar features undergoing friction can be a useful and powerful homogenization of small-scale volumetric processes, there are also cases where this practice falls short -- most notably when studying faults that grow and evolve in response to a changing tectonic environment. This is mainly due to the computational challenges associated with automating the construction of a fault-resolving conformal mesh.

Motivated by this issue, we formulate a generalization of RSF as a plastic or viscous flow law with generation, diffusion, and healing of damage that gives rise to mathematically and numerically well-behaved finite shear bands that closely mimic the behavior of the original laboratory-derived formulation (Pranger et al., submitted). The proposed formulation includes the well-known RSF laws for an infinitely thin fault as a limit case as the damage diffusion length scale tends to zero. We will show the behavior of this new bulk RSF formulation with results of high-resolution 1D and 2D numerical simulations.

Dieterich, J.H. (1979), J. Geophys. Res., 84 (B5), 2161.
Ruina, A. (1983), JGR: Solid Earth, 88 (B12), 10359–10370.
Sleep, N.H. (1995), JGR, 100 (B7), 13065–13080.
Sleep, N.H. (1997), JGR: Solid Earth 102 (B2), 2875–2895.
Roubíček, T. (2014), GJI 199.1, 286–295.
Lyakhovsky, Hamiel and Ben-Zion (2011), J. Mech. Phys. Solids, 59, 1752-1776.
Lyakhovsky and Ben-Zion (2014a), PAGeoph 171.11, 3099–3123.
Lyakhovsky and Ben-Zion (2014b), J. Mech. Phys. Solids 64, 184–197.
Lyakhovsky, Ben-Zion et al. (2016), GJI 206.2, 1126–1143.
Barbot (2019), Tectonophysics 765, 129–145.
Niemeijer and Spiers (2007), JGR 112, B10405,
Chen and Spiers (2016), JGR: Solid Earth 121, 8642–8665.
van den Ende, Chen et al. (2018), Tectonophysics 733, 273-295.
Pranger et al. (202X), ESSOAr (https://www.essoar.org/doi/10.1002/essoar.10508569.1)

How to cite: Pranger, C., Sanan, P., May, D., Le Pourhiet, L., Räss, L., and Gabriel, A.: Rate and state friction on spontaneously evolving faults, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-10412, https://doi.org/10.5194/egusphere-egu22-10412, 2022.

On-site presentation
Nicolas Riel et al.

Modelling stable mineral assemblage is crucial to calculate mineral stability relations in the Earth’s lithosphere e.g., to estimate thermobarometric conditions of exposed rocks and to quantify the fraction and composition of magma during partial melting. Accurate prediction models of stable phase are also fundamental to model trace element partitioning and to extract essential physical properties such as, fluid/melt/rock densities, heat capacity and seismic velocities. This thus forms a crucial step in linking geophysical observations with petrological constraints.

Here, we present a new Mineral Assemblage Gibbs free Energy Minimizer (MAGEMin). The package has been developed with the objective to provide a minimization routine that is easily callable and fulfilling several objectives. Firstly, the package aims to consistently compute for single point calculations at given pressure, temperature and bulk-rock composition with no needed a priori knowledge of the system. Secondly, the package has been developed for stability, performance and scalability in complex chemical systems. Finally, the code is fully parallel and we directly translate THERMOCALC formulation of solution models which yields easier and faster updates, less prone to implementation mistakes.

As a proof of concept we apply our new approach to the thermodynamic dataset for igneous systems of Holland et al. (2018). The database works in the NCKFMASHTOCr chemical system and has been updated to account for the new plagioclase model Holland et al. (2021).

How to cite: Riel, N., Kaus, B., Green, E., and Berlie, N.: MAGEMin, a new and efficient Gibbs free energy minimizer: application to igneous systems, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-11494, https://doi.org/10.5194/egusphere-egu22-11494, 2022.

On-site presentation
Boris Kaus et al.

Julia(https://julialang.org) recently emerged as a very powerful high-level computer language for (parallel) scientific computing, which allows you to “write codes like in MATLAB”, while “achieving the speed of Fortran/C”. A particular strength of Julia is that it is easy to write composable software packages that talk to each other. Here we will discuss our efforts in making Julia a development platform for geodynamic applications that significantly simplifies the process of going from a working solver to a production code which runs on massively parallel (GPU) machines.  We are working on a number of open-source packages that simplify certain steps that many geodynamics codes have in common:

  • GeoParams.jl (https://github.com/JuliaGeodynamics/GeoParams.jl) is a package in which you can specify constitutive relationships (e.g., creeplaws). It automatically handles the (non-)dimensionalization of all input parameters, includes pre-defined creep laws (e.g., dislocation and diffusion creep laws), plotting routine and includes computational routines that can be directly integrated in your code.
  • PETSc.jl (https://github.com/JuliaParallel/PETSc.jl) is the main interface from Julia to PETSc, including MPI support and automatic installations of PETSc (one of the main hurdles that existing users faced). We have recently extended the package to include an interface to DMSTAG, such that you create a staggered finite difference grid and assemble the stiffness matrix in a straightforward manner. You can use automatic differentiation tools in Julia to create the Jacobians for nonlinear equations, which again minimizes the required lines of code (compared to their C counterparts). At the same time, the full range of (nonlinear multigrid) PETSc solvers is available. This is thus very well suited to write implicit solvers.
  • ParallelStencil.jl (https://github.com/omlins/ParallelStencil.jl) and ImplicitGlobalGrid.jl (https://github.com/eth-cscs/ImplicitGlobalGrid.jl) are packages that are devoted to solving stencils in a very efficient manner on (parallel) GPU or CPU machines, which scales to very large GPU-based computers. It is particularly efficient in combination with pseudo-transient iterative solvers and allow running codes on modern architectures.
  • GeophysicalModelGenerator.jl (https://github.com/JuliaGeodynamics/GeophysicalModelGenerator.jl) is a package that gives you a simple way to collect geophysical/geological data of a certain region and combine that to construct a 3D geodynamic input model setup.

Ongoing efforts include the development of a grid generation and a marker and cell advection package that work, seamlessly with both ParallelStencil and PETSc. This will allow developers to apply both direct-iterative and pseudo-transient implicit solvers to the same problem, while only having to make minimal changes to the model setup. Combined, these packages will make the step from developing a new (nonlinear) solver to having an efficient (3D) production code much easier.

How to cite: Kaus, B., Berlie, N., Churavy, V., Cosarinsky, M., Duretz, T., Kiss, D., Kozdon, J., de Montserrat, A., Moser, L., Medinger, N., Omlin, S., Räss, L., Sanan, P., Spang, A., Thielmann, M., and Utkin, I.: How composable software tools in Julia help developing multi-physics codes in geodynamics, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-5704, https://doi.org/10.5194/egusphere-egu22-5704, 2022.

Ivan Utkin et al.

The development of highly efficient, robust, and scalable numerical algorithms lags behind the rapid increase in massive parallelism of modern hardware. In this work, we address this challenge with the accelerated pseudo-transient iterative method. This method is motivated by the physical analogy between numerical iterations and transient processes converging to a steady state.

We analytically determine optimal iteration parameters for a variety of basic physical processes such as diffusion, diffusion-reaction and non-inertial viscous fluid flow featuring Maxwell viscoelastic rheology. We further confirm the validity of theoretical predictions with numerical experiments.

We provide an efficient numerical implementation of various pseudo-transient solvers on graphical processing units (GPUs) using the Julia language. We achieve a parallel efficiency over 96% on 2197 GPUs in distributed memory parallelisation weak scaling benchmarks. 2197 GPUs allow for unprecedented terascale solutions of 3D variable viscosity Stokes flow involving over 1.2 trillion degrees of freedom.

We verify the robustness of the method by handling contrasts up to 9 orders of magnitude in material parameters such as viscosity, and arbitrary distribution of viscous inclusions for different flow configurations. Moreover, we show that this method is well suited to tackle strongly nonlinear problems such as shear-banding in a visco-elasto-plastic medium.

We additionally motivate the accessibility of the method by its conciseness, flexibility, physically motivated derivation, and ease of implementation. This solution strategy has thus a great potential for future high-performance computing applications, and for paving the road to exascale in the geosciences and beyond.

How to cite: Utkin, I., Rass, L., Duretz, T., Omlin, S., and Podladchikov, Y.: Assessing the robustness and scalability of the accelerated pseudo-transient method towards exascale computing, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-9815, https://doi.org/10.5194/egusphere-egu22-9815, 2022.

Virtual presentation
Thibault Duretz et al.

The Face-Centered Finite Volume (FCFV) is a newly developed discretisation technique that has been applied to a variety of engineering problems. This approach is based on the hybridisable discontinuous Galerkin formulation with constant degree approximations. The FCFV is particularly attractive approach since it meets numerous essential criteria for successful geodynamic modelling. It offers full geometric flexibility, natural free surface boundary condition, second order accuracy velocity-field solutions, no oscillatory pressure modes, relatively low computational cost and adequate treatment of jump conditions at material interfaces. Here we present the implementation of Poisson and Stokes solvers in the Julia computing language. Here we present the implementation of Poisson and Stokes solvers using the performant Julia language. We discuss several solving strategies including direct-iterative and iterative pseudo-transient approaches, the latter executing efficiently on Graphical Processing Units. We extend the original FCFV Stokes formulation to account for discontinuous viscosity case and discuss the implementation of complex visco-elasto-plastic rheologies.

How to cite: Duretz, T., Räss, L., and Sevilla, R.: The Face-Centered Finite Volume method for Geodynamic Modelling, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-11131, https://doi.org/10.5194/egusphere-egu22-11131, 2022.

Virtual presentation
Vinicius Silva et al.

We present a machine learning strategy to accelerate the nonlinear solver convergence for multiphase porous media flow problems. The presented approach dynamically controls an acceleration method based on numerical relaxation. The methodology is implemented and demonstrated in a Picard iterative solver; however, it can also be used with other types of nonlinear solvers. The goal of the machine learning acceleration is to reduce the number of iterations required by the nonlinear solver by adjusting the value of the relaxation factor to the complexity/physics of the system. A set of dimensionless parameters is used to train and control the machine learning. In this way, a simple two-dimensional layered reservoir can be used for training while still exploring a large portion of the dimensionless parameter space. As a result, the training process is simplified, and the machine learning model can be applied to any type of reservoir models.

We demonstrate that the presented technique dramatically reduces the number of nonlinear iterations without sacrificing the quality of the results, even for models that are far more complex than the training case. The average reduction in the number of nonlinear iterations obtained due to the presented method is 24% and the reduction in runtime is 37%. It is worth noting that the optimum value of the relaxation factor is not known a-priori and it is problem specific. Hence, having an acceleration that adapts itself to the complexity/physics of the system throughout the numerical simulation is extremely valuable and has driven several publications in multiple fields.

The method presented here provides an easy way to deal with nonlinear system of equations that does not necessitate as much effort as a custom nonlinear solver while producing outstanding results. We believe that the machine learning acceleration is not limited to the multiphase porous media flow but extendable to any other system that can be studied based on dimensionless numbers, and that a relaxation technique can be used to stabilize the nonlinear solver.

How to cite: Silva, V., Salinas, P., Jackson, M., and Pain, C.: Nonlinear solver acceleration based on machine learning applied to multiphase porous media flow, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-12398, https://doi.org/10.5194/egusphere-egu22-12398, 2022.

On-site presentation
Peter Mora et al.

Transient superstructures in mantle convection whose life and morphology vary with Rayleigh and Prandtl number have recently been demonstrated. These superstructures appear to be a two-scale phenomenon where smaller scale rolls organize into larger scale convection cells. Simulation of such superstructures requires the ability to model 3D convection in box with very large width/height ratio of order greater than 10, and with resolution to resolve the thermal boundary layer at Rayleigh numbers of 108 to 1010, respectively at least 100 height levels and 200 height levels. We achieve this with an efficient parallel implementation of the Lattice Boltzmann Method using Python which operates with high efficiency and linear speedup on thousands of cores. We present simulations with Rayleigh numbers of up to 1010 and Prandtl numbers from 1 to 100 to illustrate covering regimes from a magma ocean to solid mantle convection. We further present simulations using the LBM to model variable viscosity – specifically, temperature dependent– and illustrate the existence of pulsating plumes. We further demonstrate power law scaling between Nusselt number and Rayleigh number Nu  ~ Rag, which to first order is consistent with the Grossmann and Lohse theory.

How to cite: Mora, P., Morra, G., and Yuen, D.: Simulation of 3D transient superstructures in mantle convection and variable viscosity via the Lattice Boltzmann Method, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-9069, https://doi.org/10.5194/egusphere-egu22-9069, 2022.

Virtual presentation
Tobias Rolf et al.

The core-mantle boundary (CMB) is the most prominent compositional boundary inside the Earth. Its topography provides insight on lower mantle flow and the thermochemical structure above the CMB. Yet, CMB topography remains challenging to observe and estimates from seismology vary substantially. Numerical models of mantle convection provide complementary means to estimate CMB topography. Classically, topography is determined from the normal stresses acting on the CMB. However, this is known to face severe complications when applied to the surface boundary of the mantle, leading to non-Earth-like topographic scales and a different style of subduction. A (quasi-)free surface yields more Earth-like predictions, but for the CMB this comparison has never been made.

Here, we compare CMB topography predicted from mantle convection modelling using different treatments of the CMB. Specifically, we test the role of a ‘sticky core’, a quasi-fluid approximation the core. We compare results predicted by different codes (with either sticky core or true free base) and compare to a simple analytical case. Also, we simulate the evolution of subduction and deep thermochemical provinces to compare the topography of the (quasi-)free CMB and the free-slip approach. Initial results indicate that the sticky core approach can reproduce CMB topography reasonably well, but has rather high computational cost (grid resolution, number of particles). In analogy to the sticky air at the surface, the viscosity contrast of the sticky core layer determines the quality of predicted topography, with larger contrasts (≥103) leading to acceptable levels of artificial CMB topography. In dynamic flow cases with vigorous mantle convection, entrainment by plumes further complicates application of the sticky core, but can be tackled with an unmixing procedure. A true free base tends to better accuracy than the sticky core approach and avoids the problem with entrainment, but it also comes with additional computational costs as various forces at the CMB have to be taken into account.

How to cite: Rolf, T., Crameri, F., Heyn, B. H., and Thielmann, M.: Testing a (quasi-)free base for modelling core-mantle boundary topography, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-9133, https://doi.org/10.5194/egusphere-egu22-9133, 2022.

Nicolas Coltice et al.

Motions within the Earth mantle and tectonics constitute a single self-organized system which is cooling the planet over its geological history. Since the end of the XXth century, models of mantle convection self-generating plate tectonic behavior have progressed to a state that makes them applicable to global tectonic problems. The possibility of combining geological and geophysical data with dynamic models to retrieve the recent history of mantle flow and tectonics becomes realistic. Therefore, it is a challenge to build inverse methods to study inverse and sensitivity problems in the Earth's mantle convection. We have automatically generated the tangent-linear and the adjoint source code from the StaggYY code (Tackley, Phys. Earth Planet. Int. 171, 7-18, 2008). The Fortran code of the model was translated to the corresponding derivative codes using TAF (Transformation of Algorithms in Fortran), source-to-source translator. All codes run in parallel mode, using MPI (Message Passing Interface). The economic taping strategy of TAF, including re-computations, and checkpointing, helps to keep the memory footprint of the adjoint code low and the performance high. We highlight some key features of the automatic differentiation, evaluate the performance of the adjoint code, and show first results from 2D and 3D sensitivity fields, focusing on the relationships between temperature in the mantle and tectonics. Ultimately the addjoint code shall be applied to inversion and assimilation problems using a bayesian framework.

How to cite: Coltice, N., Blessing, S., Giering, R., and Tackley, P.: Automatic generation of the adjoint of the StagYY mantle convection model, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-9232, https://doi.org/10.5194/egusphere-egu22-9232, 2022.

On-site presentation
Zhouji Liang et al.

Gravity is one of the most widely used geophysical data types in subsurface exploration. In the recent developments of stochastic geological modeling, gravity data serves as an additional constraint to the modeling construction and can be included in the modeling process as the likelihood function in a Bayesian workflow. A fast but also precise forward gravity simulation is key to the success of the geological modeling inverse problem.

In this study, we present a gravity kernel method, which is based on the widely adopted analytical solution on a discretized grid. As opposed to a globally refined regular mesh, we construct local tensor grids for each sensor, respecting the gravimeter locations and the local sensitivities. The kernel method is efficient in terms of both computing and memory use for meshless implicit geological modeling approaches. This design makes the method well suited for many-query applications like Bayesian machine learning using gradient information calculated from Automatic Differentiation (AD). Optimal grid design without knowing the underlying geometry is not straightforward before evaluating the model. Therefore, we further provide a novel perspective on a refinement strategy for the kernel method based on the sensitivity of the cell to the corresponding receiver. Synthetic results are presented and show superior performance compared to the traditional spatial convolution method.

How to cite: Liang, Z., De La Varga, M., and Wellmann, F.: Gravity kernel method for implicit geological modeling, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-3730, https://doi.org/10.5194/egusphere-egu22-3730, 2022.

Meissam L. Bahlali et al.

Density-driven flows in porous media are frequently encountered in natural systems and arise from the gravitational instabilities introduced by fluid density gradients. They have significant economic and environmental impacts, and numerical modelling is often used to predict the behaviour of these flows for risk assessment, reservoir characterisation or management. However, modelling density-driven flow in porous media is very challenging due to the nonlinear coupling between flow and transport equations, the large domains of interest and the wide range of time and space scales involved. Solving this type of problem numerically using a fixed mesh can be prohibitively expensive.  Here, we apply a dynamic mesh optimisation (DMO) technique along with a control-volume-finite element method to simulate density-driven flows. DMO allows the mesh resolution and geometry to vary during a simulation to minimize an error metric for one or more solution fields of interest, refining where needed and coarsening elsewhere. We apply DMO to the Elder problem for several Rayleigh numbers. We demonstrate that DMO accurately reproduces the unique two-dimensional (2D) solutions for low Rayleigh number cases at significantly lower computational cost compared to an equivalent fixed mesh, with speedup of order x16. For unstable high Rayleigh number cases, multiple steady-state solutions exist, and we show that they are all captured by our approach with high accuracy and significantly reduced computational cost, with speedup of order x6. The lower computational cost of simulations using DMO allows extension of the high Rayleigh number case to a three-dimensional (3D) configuration and we demonstrate new steady-state solutions that have not been observed previously. Early-time, transient 3D patterns represent combinations of the previously observed, steady-state 2D solutions, but all evolve to a single, steady-state finger in the late time limit.

How to cite: Bahlali, M. L., Salinas, P., and Jackson, M. D.: Dynamic mesh optimisation for efficient numerical simulation of density-driven flows: Application to the 2- and 3-D Elder problem, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-12156, https://doi.org/10.5194/egusphere-egu22-12156, 2022.