Following the long-standing paradigm in HPC, computational geodynamic codes have been typically written in high-level statically typed and compiled languages, namely C/C++ and Fortran. The low productivity rates of these languages led to the so-called two-language problem, where dynamic languages such as Python or MATLAB are used for prototyping purposes, before porting the algorithms to high-performance languages. The Julia programming language aims at bridging the productivity and advantages of such dynamic languages without sacrificing the performance provided by static languages. The high performance of Julia, combined with high-productivity rates and other powerful tools, such as advanced meta-programming (i.e. code generation), make Julia a suitable candidate for the next generation of HPC-ready scientific software.

We introduce the open-source and Julia-written package JustRelax.jl (https://github.com/PTsolvers/JustRelax.jl) as a way forward for the next generation of geodynamic codes. JustRelax.jl is a production-ready API for a collection of highly-efficient numerical solvers (Stokes equations, diffusion, etc.) based on the embarrassingly parallel pseudo-transient method. We rely on ParallelStencil.jl (https://github.com/omlins/ParallelStencil.jl), which leverages the advanced meta-programming capabilities of Julia to generate efficient computational kernels agnostic to the back-end system (i.e. Central Processing Unit (CPU) or Graphics Processing Unit (GPU)). Using ImplicitGlobalGrid.jl (https://github.com/eth-cscs/ImplicitGlobalGrid.jl) to handle the MPI and CUDA-aware MPI communication, these computational kernels run seamlessly in local shared-memory workstations and distributed memory and multi-GPU HPC systems with little effort for the front-end user.

Efficient computation of the (local) physical properties of different materials is another key feature required in geodynamic codes, for which we employ GeoParams.jl (https://github.com/JuliaGeodynamics/GeoParams.jl). This package provides lightweight, optimised, and reproducible computation of different material properties (e.g. advanced rheological laws, density, seismic velocity, etc.), amongst other available features. GeoParams.jl is also carefully designed to support CPU and GPU devices, and be fully compatible with other external packages, such as ParallelStencil.jl and existing auto-differentiation packages.

We finally show high-resolution examples of geodynamic models run on multi-GPU HPC systems employing the presented open-source Julia tools.

How to cite: de Montserrat Navarro, A., Kaus, B., Räss, L., Utkin, I., and Tackley, P.: Using Julia for the next generation of HPC software for geodynamic modelling, EGU General Assembly 2023, Vienna, Austria, 23–28 Apr 2023, EGU23-9242, https://doi.org/10.5194/egusphere-egu23-9242, 2023.

Mt. Etna, located in the north-eastern area of Sicily (Italy), is one of the most active and hazardous strato-volcano in the world, both in terms of paroxysmal events and continuous effusive activity from the summit area and hazardous flank eruptions. Long-term processes of deep magma recharge and storage within the upper crust, passive magma ascent along pre-existing weaknesses, and forceful dyke intrusions allow magma to rise to the surface. Past studies provided evidence supporting the view that the interplay between magma dynamics and storage and the thermomechanical response of the host medium control magma rise and the brittle seismic response of the volcano basement and edifice.

To further investigate this interplay, we have performed 3D thermomechanical simulations of the present-day state of the volcano using the Lithosphere and Mantle Evolution Model (LaMEM) code. The model is built between the volcanic surface and 30km depth and includes realistic topography. Magma storage zones within the model are inferred from seismic tomography and seismic source studies at ~30km (deep storage) and between 3 and 6 km (shallow storage). The characteristics of the molten zones are calibrated by physical and mechanical properties determined for the main representative lithologies (carbonates, basalts, clays) and the corresponding rheological laws. As we are interested in the present-day dynamics of the volcano, we ran our models for just a few timesteps to gain surface velocity and displacement data.

The LaMEM framework allows retrieving both deformation and gravity responses to the final model. These responses will be fit to real GPS, InSAR, and gravity data to define the most realistic properties of the Etna feeding systems. Future steps will include tectonic forces contributing to the sliding of the eastern volcanic flank in the simulation and propagation of seismic waves in the final model suitable to fit existing seismic data.

How to cite: Bensing, M., Vinciguerra, S., and De Siena, L.: Seismic response to volcanic processes at Mount Etna: coupling thermomechanical simulations with seismic wave-equation modeling, EGU General Assembly 2023, Vienna, Austria, 23–28 Apr 2023, EGU23-9515, https://doi.org/10.5194/egusphere-egu23-9515, 2023.

Strain localization is a crucial process for lithosphere and mantle deformation as it allows for the formation of faults and shear zones that enable plate tectonics. In the crust, strain localization usually occurs via brittle failure (i.e., breaking the rock). The deeper and/or hotter the setting, the less likely brittle failure becomes as the critical stress increases with the increasing overburden pressure while the temperature-dependent rheology of rocks limits the stresses that can be accumulated before being relaxed by slow, viscous flow.

Yet, we do observe fast and localized deformation (i.e., earthquakes) at depths of several hundred kilometers. These deep earthquakes either require local differential stresses of several Gigapascal (GPa) to trigger brittle failure or a different, ductile failure mechanism that significantly reduces rock strength while at the same time creating highly localized shear zones. Here, we investigate the feedback loop of visco-elastic deformation and shear heating to determine whether their combination can lead to a localized viscosity reduction and allow for fast slip.

Modeling this feedback loop and the accompanying strong localization of deformation poses a challenge for continuum modeling approaches, in particular when highly nonlinear rheologies such as dislocation creep and low-temperature plasticity are employed. Here, we present a collection of 1D and 2D numerical codes written in the Julia language which use the pseudo-transient approach and graphical processing unit (GPU) computing to model the process of ductile localization and thermal runaway in a simple-shear setting. Our models employ a nonlinear, visco-elastic rheology, including grain-size-dependent diffusion creep, stress-dependent dislocation creep and low-temperature plasticity. We find that the combination of the aforementioned mechanisms is sufficient for deformation to localize on a small perturbation and then propagate through the model similar to a brittle rupture. Our models show that low-temperature plasticity acts as a stress-limiting mechanism that facilitates numerical stability during thermal runaway. In a systematic series of models, we investigate under which conditions thermal runaway occurs and which role each of the rheological components plays in the localization process.

How to cite: Spang, A., Thielmann, M., and Kiss, D.: GPU-based numerical models of rapid ductile strain localization due to thermal runaway, EGU General Assembly 2023, Vienna, Austria, 23–28 Apr 2023, EGU23-11594, https://doi.org/10.5194/egusphere-egu23-11594, 2023.

On Earth, different geodynamic features form in response to a tectonic event. Continental plateaus, such as the Tibetan Plateau, are formed in a collisional environment and they are characterized by an unusually large crustal thickness, which generates lateral variations of gravitational potential energy per unit area (GPE). These GPE variations cause the thickened crust to flow apart and thin by gravitational collapse. Although plateau and lowland are in isostatic equilibrium, the lateral GPE variations must be balanced by horizontal differential stresses, which prevent the plateau from flowing-apart instantaneously. However, the magnitude and distribution of differential stress around plateau corners for three-dimensional (3-D) spherical geometries relevant on Earth remain disputed. Due to the ellipticity of the Earth, the lithosphere is mechanically analogous to a shell, characterized by a double curvature. Shells exhibit fundamentally different mechanical characteristics compared to plates, having no curvature in their undeformed state. Understanding the magnitude and the spatial distribution of strain, strain-rate and stress inside a deforming lithospheric shell is thus crucial but technically challenging. Resolving the stress distribution in a 3-D geometrically and mechanically heterogeneous lithosphere requires high-resolution calcuations and high-performance computing.

Here, we present numerical simulations solving the Stokes equations under gravity. We employ the accelerated pseudo-transient finite-difference (PTFD) method, which enables efficient simulations of high-resolution 3-D mechanical processes relying on a fast iterative and implicit solution strategy of the governing equations. The main challenges are to guarantee convergence, minimize the iteration count and ensure optimal execution time per iteration. We implemented the PTFD algorithm using the Julia language. The Julia packages ParallelStencil.jl and ImplicitGlobalGrid.jl enable optimal parallel execution on mulitple CPUs and GPUs and ideal scalability up to thousands of GPUs.

The aim of this study is to quantify the impact of different lithosphere curvatures on the resulting stress field. To achieve this, we use a simplified plateau geometry and density structure implemented in a spherical coordinate system. The curvature is modified by varying the radius of the coordinates system, without altering the initial geometry. We particularly focus on stress magnitudes and distributions in the corner regions of the plateau.

How to cite: Macherel, E., Podladchikov, Y., Räss, L., and Schmalholz, S. M.: GPU-based finite-difference solution of 3-D stress distribution around continental plateaus in spherical coordinates, EGU General Assembly 2023, Vienna, Austria, 23–28 Apr 2023, EGU23-12450, https://doi.org/10.5194/egusphere-egu23-12450, 2023.

While the last decade has seen significant progress in thermo-mechanical modelling of complex multiphase systems, the coupling with petrological thermodynamic modelling approaches, when addressed at all, remains a difficult task. First, most phase equilibria modeling tools have been developed with the primary focus to produce phase diagrams (e.g., Perple_X, Theriak_Domino, geoPS, MELTS) and do not offer useful interfaces for (parallel) geodynamic codes. Second, phase equilibrium modelling is generally achieved by solving a Gibbs energy minimization problem. This problem is computationally challenging as it involves solving a nested optimization problem subject to both equality and inequality constraints. As a result, the single point calculation of stable phase equilibrium is slow, and to our knowledge, >150 milliseconds for a compositional system involving a large number of chemical components. This limitation effectively precludes direct coupling of phase equilibria calculation with geodynamic models, which requires performing 1000s to 100'000s of such calculations every timestep.

We have recently developed a new open-source code, MAGEMin, that improves on part of this. In MAGEMin 75 to 90% of the computation time is dedicated to local minimization of solution models. Therefore, it becomes critically important to improve the minimization time of individual solution phase models to further speed-up phase-equilibria computational time.

Here, we present a reformulation of the solution phase model from Holland et al., (2018) that eliminates both equality and inequality constraints. Eliminating these constraints allows the utilization of faster unconstrained optimization methods, thus yielding much higher performance and stability. We compare the accuracy and performance of several unconstrained gradient-based optimization methods namely the conjugated gradient (CG), the Broyden-Fletcher-Goldfarb-Shanno (BFGS) methods and a hybrid combination (CG-BFGS).

How to cite: Riel, N., Kaus, B., and Green, E.: An unconstrained formulation for thermodynamic complex solution phase minimization, EGU General Assembly 2023, Vienna, Austria, 23–28 Apr 2023, EGU23-13828, https://doi.org/10.5194/egusphere-egu23-13828, 2023.

Regional geodynamic models require to impose boundary conditions that best represent the physical information exchanged between the modelled and a larger, non-modelled domain. Depending on the nature of the physical information exchange, the internal evolution of the regional system may differ. Nevertheless, the first and foremost observation is that the deformation in tectonic plates boundaries is three-dimensional, i.e., non-cylindrical, oblique.

To model 3D non-cylindrical deformation, regional geodynamic models mostly use initial conditions through oblique or offset weak zones together with cylindrical boundary conditions implying free slip. However, the problem with the free-slip boundary condition is that it enforces cylindrical behaviours in the vicinity of the boundary, limiting the obliquity of the whole system or forcing to consider very large domains to avoid a too strong influence of the boundary condition.

A way to work around this problem is to impose obliquity through boundary conditions. Until now, the main approach to impose oblique boundary conditions involves strong Dirichlet constraints, i.e., directly providing the solution for the velocity (or displacement) along the boundary.

However, the choice of velocity values can lead to arbitrarily imposed velocity gradients particularly in the tangential direction of the boundary when the velocity vectors point in different directions. Such boundary effects can then influence the strain localization and produce non-physical results.

In this work, we propose a formulation to impose oblique boundary conditions by enforcing the velocity direction but without constraining the magnitude of the velocity vectors. We seek to impose a slip-type boundary condition. The formulation is a generalszation of Nitsches’ method (Nitsche, 1971) thereby allowing Navier-slip constraints to be enforced independently of the orientation domain boundary. We refer to this new formulation as the generalised Navier-slip boundary condition..

In order to demonstrate that the method works as well as to illustrate the differences it produces on the evolution of a geodynamic system compared to the use of more classical boundary conditions, we show two 3D oblique rift models. The first uses Dirichlet boundary conditions and the second uses the generalised Navier-slip method to enforce an oblique extension at 45°.

The models show differences not only along and near their boundaries but also in the centre of the modelled domain during its tectonic evolution in terms of strian localization, basin architecture and topography. Moreover, the model using the generalised Navier-slip method to impose oblique extension shows a more natural evolution of the strain localization and tectonic features as the velocity along and near boundaries can vary in time and space to adapt to the internal evolution of the model.

Finally, we show that the generalised Navier-slip method provides a better approach to impose oblique boundary conditions than the classical methods as it does not require to impose an arbitrary velocity function directly into the solution.

Nitsche, J., 1971. Über ein variationsprinzip zur lösung von dirichlet-problemen bei verwendung von teilräumen, die keinen randbedingungen unterworfen sind. Abhandlungen aus dem Mathematischen Seminar der Universität Hamburg 36, 9–15. doi:10.1007/BF02995904.

How to cite: Jourdon, A., May, D. A., and Gabriel, A. A.: Generalization of the Nitsche method to apply oblique boundary conditions in regional geodynamic models, EGU General Assembly 2023, Vienna, Austria, 23–28 Apr 2023, EGU23-14864, https://doi.org/10.5194/egusphere-egu23-14864, 2023.

The processes that govern rock (trans)formation (deposition, deformation, segregation, metamorphism) can result in the development of layering and rock fabrics. Rocks can thus exhibit extrinsic or intrinsic anisotropy at various spatial scales. Anisotropy has important mechanical consequences, in particular, for strain localisation in the lithosphere. This effect is typically not included in geodynamic models. Mechanical anisotropy can be modelled by explicitly modelled by numerically resolving layers of different strengths. Due to the expensive computational cost, this approach is not suitable for large scale geodynamic models. The latter may rather benefit from an upscaling approach that involves anisotropic constitutive laws.  To model the evolution of such material Mühhlaus, (2002) proposed the use of the director vector which corresponds to a single orientation that is changing throughout the process of deformation. We have implemented visco-elasto-plastic anisotropic constitutive laws and the director vector approach in the geodynamic simulation tool MDoodz7.0. Here we present  the rheological implementation, we show some simple simulations involving anisotropic flow and discuss the potential role of anisotropy for large-scale geodynamic processes.

How to cite: Kulakov, R., Halter, W., Schmalholz, S., and Duretz, T.: Modelling lithosphere deformation with non-linear anisotropic constitutive models, EGU General Assembly 2023, Vienna, Austria, 23–28 Apr 2023, EGU23-17307, https://doi.org/10.5194/egusphere-egu23-17307, 2023.

Understanding the dynamics of magmatic systems requires numerical models that take the physics of the involved processes into account and allows interpreting geophysical and geological data in a consistent manner.  In climate science, a similar venture started over 5 decades ago with the generation of the first quantitative climate models, which has been indispensable in our understanding of the ongoing climate change. A similar effort for magmatic systems does not yet exist, even when many processes can already be described quantitatively.

Here, we will discuss recent progress towards creating a modelling framework to simulate magmatic systems, developed as part of the ERC MAGMA project. We initiated several open-source packages in the Julia programming language that significantly simplifies creating new codes that simulate different processes and run on both workstations and high-performance GPU systems.

This makes it straightforward to create a 3D model of a particular system taking available data into account (using GeophysicalModelGenerator.jl), use that as input for 3D models that link uplift/gravity data with dynamic models (using LaMEM), or simulate the thermal evolution and zircon age distribution following the intrusion of dikes & sills (using MagmaThermoKinematics.jl). One can easily switch the employed rheologies/parameterisations in the FEM or finite difference simulations, create synthetic seismic velocity models from the output (using GeoParams.jl) or account for the evolving chemistry of the magmatic system (using MAGEMin_C.jl).

In this presentation, we will discuss implementation details and show that the use of GeoParams, for example, slows down pseudo-transient codes (as expected) but not substantially, whereas it results in much shorter codes.

How to cite: Kaus, B., Kiss, D., de Montserrat, A., Riel, N., Berlie, N., and Spang, A.: Towards integrated numerical models of lithospheric-scale magmatic systems, EGU General Assembly 2023, Vienna, Austria, 23–28 Apr 2023, EGU23-14034, https://doi.org/10.5194/egusphere-egu23-14034, 2023.