The complexity 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 further constraining 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 processes.
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. Theoretical 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 new numerical modelling approaches
- Utilisation of coupled models to explore nonlinear interactions
Mantle convection simulations are central to understanding the thermal evolution of rocky planets. However, their high computational cost limits the feasibility of extensive parameter studies needed to constrain models with observations. While scaling laws provide low-cost alternatives, they are limited in the physical processes they can capture and typically predict only reduced quantities rather than the full spatio-temporal fields.
Machine learning (ML) promises to accelerate mantle convection simulations, yet purely data-driven approaches can fail to match the accuracy and stability of numerical solvers, even when trained on thousands of simulations. To address this, we propose a physics-based ML framework that combines neural networks with numerical time integration. The ML model predicts creeping-flow velocities as a function of temperature, thereby bypassing the numerical solution of the Stokes equations, which poses the primary computational bottleneck in mantle convection simulations. Mass conservation is enforced as a hard constraint through a stream-function formulation. The predicted velocity field is then used by a finite-volume solver to advect and diffuse temperature forward in time.
The model is trained on temperature–velocity snapshots from 94 two-dimensional simulations of statistically steady mantle convection, where three parameters are varied: internal heating as well as pressure- and temperature-dependence of viscosity. Compared to the direct numerical solver, our model is 89 times faster. For some parameter combinations, the model outperforms an under-relaxed iterative numerical solver in speed and accuracy, further underscoring the potential of ML in geodynamics. Ablation studies demonstrate the importance of mass conservation, learned boundary padding, and loss scaling in achieving stable and accurate predictions over long time-integration scales.
Despite being trained exclusively on snapshots from statistically steady simulations, the model successfully performs thermal evolution, demonstrating generalization in this unseen setting. Performance degrades, however, when additional compressibility effects are introduced at inference or when initial conditions deviate substantially from the training data, highlighting directions for future improvements.
How to cite: Agarwal, S., Bekar, A., Hüttig, C., Greenberg, D., and Tosi, N.: Physics-based machine learning for mantle convection , EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-20015, https://doi.org/10.5194/egusphere-egu26-20015, 2026.
Mantle convection ultimately drives Earth’s tectonics. Therefore, it is key to truthfully reconstruct the past mantle flow and associated dynamic evolution. Mathematically, mantle convection could be approximated as an advection-diffusion problem with a high effective viscosity. When supplemented with appropriate initial and boundary conditions, the past mantle flow can be numerically reconstructed. However, constructing physically and geologically valid initial conditions for Earth remains challenging due to the lack of direct observations. The adjoint inversion scheme demonstrates great potential for constraining such valid initial conditions through assimilating present-day mantle structures revealed by seismic tomography and other geophysical and geological observations. In practice, this method iteratively solves the forward and adjoint conservation equations until the mismatch between the model prediction and observation becomes small enough.
However, traditional grid-based energy equation solvers suffer from numerical oscillations due to instability of the advection term. Numerical errors will accumulate through the iterative process, resulting in a final mismatch locked to a local minimum rather than the global minimum. To address this problem, we introduce the particle-in-cell (PIC) method, which solves the advection problem using Lagrangian tracers, to improve the solutions to both the forward and adjoint energy equations. Here, we present a series of synthetic experiments to validate the new method while exploring its limitations. Our preliminary results show that incorporating the PIC method into the adjoint inversion scheme significantly reduces numerical oscillations, thus improving the reliability of the reconstructed initial conditions. Further experiments show that the updated adjoint inversion scheme converges much faster and produces a smaller mismatch between the prediction and observation, thereby reducing computational cost while improving the validity of numerical results. This updated adjoint scheme will be applied to real Earth problems and the result will improve our understanding of complex geodynamic problems.
How to cite: Liu, G., Liu, L., Zhao, H., Cao, Z., and Dong, H.: Improving the performance of the adjoint inversion scheme for mantle convection with the particle-in-cell method, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-6917, https://doi.org/10.5194/egusphere-egu26-6917, 2026.
How to cite: Burkhart, A., Wohlmuth, B., Schuberth, B., Robl, G., Papanagnou, I., Kohl, N., and Zawallich, J.: Matrix-free, robust linear solvers for realistic, large scale, high-contrast forward mantle convection simulations, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-7979, https://doi.org/10.5194/egusphere-egu26-7979, 2026.
The compositional record in mineral grains is used to characterize temperatures (geothermometry), pressures (geobarometry), dates (geochronology) and durations (diffusion chronometry). However, the accessible record is limited by the timing of formation of a grain, either by primary crystallization or by recrystallization due to mechanical or chemical influences. In this presentation we would like to focus on recrystallization. One effect which causes recrystallization is the effect of lattice strain caused by the diffusion of ions of a different size in a mineral. While the lattice strain itself has been studied [1] in detail, the effect on recrystallization has only been shown qualitatively in experimental settings [2, 3] and field observations [4]. A model for a single crystal has been developed recently [5], but the effect on a polycrystalline rock has not been studied yet.
In this presentation we build on the model introduced by Haddenhorst et al. [5] for the evolution of a single olivine crystal. We simulate multiple crystals simultaneously and introduce a maximum volume constraint. This enables us to predict the composition and crystal sizes in a rock bearing multiple olivine crystals in a small volume (e.g. a dunite). Results of simulations and comparisons with real world examples will be shown.
References:
[1]: Blundy, J., & Wood, B. (1994). Prediction of crystal–melt partition coefficients from elastic moduli. Nature, 372 (6505), 452–454.
[2]: Bestmann, M., Pennacchioni, G., Grasemann, B., Huet, B., Jones, M. W. M., & Kewish, C. M. (2021). Influence of deformation and fluids on ti exchange in natural quartz. Journal of Geophysical Research: Solid Earth, 126 (12), e2021JB022548. Retrieved from https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2021JB022548
[3]: Beyer, C., & Chakraborty, S. (2021). Internal stress-induced recrystallization and diffusive transport in catio3-pbtio3 solid solutions: A new transport mechanism in geomaterials and its implications for thermobarometry, geochronology, and geospeedometry. American Mineralogist: Journal of Earth and Planetary Materials, 106 (12), 1940–1949.
[4]: Nachlas, W., & Hirth, G. (2015). Experimental constraints on the role of dynamic recrystallization on resetting the ti-in-quartz thermobarometer. Journal of Geophysical Research: Solid Earth, 120 (12), 8120–8137.
[5]: Haddenhorst, H. H., Chakraborty, S., & Hackl, K. (2023). A model for the evolution size and composition of olivine crystals. Proceedings in Applied Mathematics and Mechanics, 00, e202300081. https://doi.org/10.1002/pamm.202300081
How to cite: Haddenhorst, H. H., Waimann, J., Chakraborty, S., and Hackl, K.: A material model for compositional resetting due to coupled mechanical and chemical effects, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-9491, https://doi.org/10.5194/egusphere-egu26-9491, 2026.
Accurately incorporating realistic phase transitions in geodynamic models is crucial but challenging. Phase transitions strongly influence mantle convection, as their effects on buoyancy can hinder or accelerate slabs and plumes. In a heterogeneous mantle, different mineral assemblages undergo phase transitions at different depths, leading to lateral buoyancy variations that can cause specific compositions to stagnate or accumulate within characteristic depth ranges. However, complex phase relations, abrupt changes in material properties, and the release or absorption of latent heat pose significant challenges for modeling phase transitions. Dannberg et al. (2022) addressed these challenges by solving the energy equation in terms of entropy rather than temperature, allowing it to capture realistic phase changes in geodynamic models. However, this method was limited to chemically homogeneous systems.
Building on our earlier work, we now present a new multi-component formulation that extends the entropy method to systems with compositional heterogeneities. Our formulation assumes thermal equilibrium below the scales resolved by our mesh, i.e., all components share a single temperature at each point represented in the model. Pressure changes during advection produce an isentropic temperature change. As different chemical components may have different isentropic temperature gradients, which would imply different temperatures for each component at the same location, our multicomponent formulation involves a thermal equilibration step. Using a series of benchmarks and test cases, we show that our implementation in the geodynamic modeling software ASPECT satisfies the conservation of total energy and captures phase transitions self-consistently, regardless of their sharpness, during advection of chemical heterogeneities.
We show the applicability of our new formulation in a series of global convection models. We compare: (1) a single-component pyrolite model, and (2) a two-component mechanical mixture of basalt and harzburgite with the same pyrolitic bulk composition. Our results reveal that small differences in the ringwoodite to bridgmanite + ferropericlase transition between these assemblages with the same composition can significantly affect slab and plume stagnation. Our results highlight the importance of accurately capturing the full effects of phase transitions in a chemically heterogeneous mantle, and our approach enables new investigations into how planetary mantles evolve.
References: Dannberg, J., Gassmöller, R., Li, R., Lithgow-Bertelloni, C., & Stixrude, L. (2022). An entropy method for geodynamic modelling of phase transitions: capturing sharp and broad transitions in a multiphase assemblage. Geophysical Journal International, 231(3), 1833-1849.
How to cite: Li, R., Dannberg, J., Gassmöller, R., and Myhill, R.: A Multi-Component Entropy Method for Modeling Phase Transitions in a Heterogeneous Mantle, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-12685, https://doi.org/10.5194/egusphere-egu26-12685, 2026.
Solifluction is an important geomorphological process, mainly driven by seasonal freeze–thaw cycles that mobilize the upper soil layers and contribute significantly to landscape evolution. Despite its relevance, there is currently no numerical package available for simulating solifluction in a physically based manner. In this study, we present a new three-dimensional numerical model, soli3d, developed to simulate this phenomenon by deriving the governing physical partial differential equations and explicitly solving them using numerical discretization.
The model integrates mass conservation, momentum conservation, and heat transfer. The thawed soil layer is treated as a viscous fluid, and the corresponding momentum equations are derived following principles from fluid mechanics. Topographic evolution is tracked using a Volume of Fluid (VOF) method. Soil temperature profiles are computed by solving the heat transfer equation in the vertical direction. The viscosity of the thawed soil is assumed to depend on temperature, soil and vegetation properties. The soli3d software package uses the finite difference method (FDM) to discretize and solve the governing equations. The vertical domain is represented by multiple layers, while the horizontal plane is discretized using uniform, structured rectangular grids. The soli3d model (see soli3d link [https://github.com/computationalgeography/soli3d]) is available as an open-source Python implementation utilizing the LUE environmental modelling framework (see link [https://zenodo.org/records/16792016]) for efficient parallel and distributed computing.
To assess the model performance, soli3d is validated using a set of benchmark test cases. The results demonstrate the capability of the model to reproduce key characteristics of solifluction processes and provide a foundation for future applications.
How to cite: Faraji Gargari, S., Draebing, D., van Beek, R., Schmitz, O., Karssenberg, D., and Eichel, J.: A Three-Dimensional Solifluction Model: From Governing Equations to a Numerical Model Package, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-8105, https://doi.org/10.5194/egusphere-egu26-8105, 2026.
Understanding the Earth’s internal structure remains a major challenge, as traditional geophysical methods face ambiguities in linking seismic observations to temperature, composition, or mass density variations. Atmospheric neutrinos offer a complementary probe: while traversing the Earth, they undergo flavor oscillations that depend on the local electron density, which reflects both mass density and composition. Here, we present EarthProbe, a forward-modeling framework for neutrino propagation and detection, providing a methodology to quantify neutrino sensitivities to the Earth’s interior. Using EarthProbe, we assess the detectability of localized electron-density perturbations, taking the Mantle Transition Zone (410–670 km depth) and the core as case studies. We consider idealized next-generation detectors representing fundamental sensitivity limits and the state-of-the-art instruments KM3NeT/ORCA, Hyper-Kamiokande, and DUNE. While studying the core is within reach of current detector capabilities, probing the MTZ would require improved detector performance. Our methodology lays the foundation for future joint inversion of neutrino and seismic data, providing a framework to advance Earth tomography with neutrinos.
How to cite: Goos, I. A., Fuji, N., Van Elewyck, V., Durand, S., Coelho, J. A. B., Mittelstaedt, E., and Deniz, Y. A.: Sensitivity of neutrino oscillations to the Earth’s interior properties, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-18386, https://doi.org/10.5194/egusphere-egu26-18386, 2026.
How to cite: de Montserrat Navarro, A., Kaus, B., and Duretz, T.: RheologyCalculator.jl: automatically solving arbitrary constitutive relationships, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-3772, https://doi.org/10.5194/egusphere-egu26-3772, 2026.
Understanding and quantifying the flow of multiple pore fluid phases is a key component to understanding many natural and engineered reservoir processes. A well-known example is the migration of hydrocarbon phases in sedimentary reservoirs, including natural migration, conventional production, and enhanced hydrocarbon recovery. As Carbon Capture and Storage (CCS) emerged as a potential strategy to combat climate change, understanding multiphase flow processes of supercritical CO2 and the various pore fluids is at the forefront of scientific interest.
Our governing equations are based on the conservation of mass and momentum in two immiscible fluid phases. The fluid phases may be compressible or incompressible. We assume the inertial terms to be negligible in both phases and momentum transfer to happen in the Darcy-flow regime. Multiple fluid phase effects are introduced in the mathematical model by saturation-dependent relative permeabilities and fluid viscosities, which uniquely determine the total mobility, relative mobility, and fractional fluid flow curves. Here we consider capillary effects negligible. This results in two independent equations: an elliptic equation for fluid pressure and a hyperbolic equation for saturation.
We choose a staggered-grid-based finite-difference discretization. The elliptic fluid pressure equation is solved using pseudo-transient iterations. The hyperbolic saturation equation is solved using a conservative first-order upwind scheme. Multiple coupling and time-stepping options (e.g., explicit/implicit, first order/higher order) are tested against an analytical solution. The numerical scheme can be implemented on GPUs in a straightforward manner using the ParallelStencil.jl package in Julia. We will provide some examples of various reservoir applications, which will also be used as performance benchmarks. Finally, we will discuss the potential of the presented numerical scheme to be implemented in a broader THMC framework.
How to cite: Kiss, D., Wang, L. H., and Yarushina, V.: A simple, conservative, staggered grid-based pseudo-transient scheme implemented on GPUs to solve two-fluid phase flow in porous reservoirs, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-20424, https://doi.org/10.5194/egusphere-egu26-20424, 2026.
The Lithosphere and Mantle Evolution Model (LaMEM; Kaus et al., 2016) is a parallel thermo-mechanical code designed for a wide range of geomechanical and geotechnical applications involving rocks and soils with highly nonlinear visco-elasto-plastic rheologies. The code primarily targets the simulation of geological processes related to crustal dynamics and lithosphere-mantle interactions, while also supporting several geotechnical applications. LaMEM employs a stable, lightweight staggered-grid finite-difference discretization combined with a marker-in-cell technique for material advection. The code is built on the highly efficient PETSc (Portable, Extensible Toolkit for Scientific Computation) library (Balay et al., 2025), enabling near-optimal scaling on massively parallel computing architectures. PETSc further provides a comprehensive suite of algorithms and data structures for the solution of linear and nonlinear systems. To accelerate large-scale simulations, custom Galerkin coarsening operators have been integrated into the PETSc multigrid framework. LaMEM is distributed under the MIT license and is openly available on GitHub at https://github.com/UniMainzGeo/LaMEM.
In this contribution, we announce an upcoming major release of LaMEM (v3.0.0), which introduces substantial enhancements in both functionality and computational performance. We provide a concise overview of the key developments in this release, including: periodic boundary conditions; a hybrid Galerkin/matrix-free multigrid framework; two-dimensional multigrid coarsening; block factorization and wBFBT preconditioners; improved default solver configurations; and automatic stopping tolerances.
References
Balay, S., Abhyankar, S., Adams, M. F., Benson, S., Brown, J., Brune, P., Buschelman, K., Constantinescu, E., Dalcin, L., Dener, A., Eijkhout, V., Faibussowitsch, J., Gropp, W. D., Hapla, V., Isaac, T., Jolivet, P., Karpeev, D., Kaushik, D., Knepley, M. G., Kong, F., Kruger, S., May, D. A., McInnes, L. C., Mills, R. T., Mitchell, L., Munson, T., Roman, J. E., Rupp, K., Sanan, P., Sarich, J., Smith, B. F., Suh, H., Zampini, S., Zhang, H., Zhang, H., Zhang, J., 2025. PETSc/TAO Users Manual, Tech. Rep. ANL-21/39 - Revision 3.23. Argonne National Laboratory.
Kaus, B. J. P., Popov, A. A., Baumann, T. S., Pusok, A. E., Bauville, A., Fernandez, N., Collignon, M., 2016. Forward and Inverse Modelling of Lithospheric Deformation on Geological Timescales. NIC Series, 48, 299-307.
How to cite: Popov, A. and Kaus, B.: Lithosphere and Mantle Evolution Model (LaMEM): v3.0.0, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-22075, https://doi.org/10.5194/egusphere-egu26-22075, 2026.
Strain localization in the crust surrounding magmatic systems plays a critical role in controlling magma storage, transport, and eruption pathways. In geoscientific modeling, plastic deformation remains a topic of active discussion, with studies focusing on reducing mesh-dependence of shear bands, implementing tensile (mode-1) plasticity, or analyzing fault orientations. While large-scale geodynamic models typically resolve shear zones on the order of kilometers, e.g., in subduction zones, the faults around magma chambers or sills occur at much smaller scales. These smaller-scale structures are crucial, as they form open pathways through which magma can propagate upward, potentially producing fissure eruptions (e.g., Iceland) or major explosive eruptions through a central conduit (e.g., Krakatoa).
Until recently, resolving these different scales in a single simulation was computationally expensive and nearly impossible. However, with the improvement of Graphical Processing Units (GPUs) and the focus on GPU-accelerated high-performance computing (HPC), these fine scales are now accessible at feasible computational cost. In this study, we use the thermo-mechanical visco-elasto-plastic code JustRelax.jl to model evolving shear bands around a circular inclusion in 2D. We simulate different grid resolutions ranging from 64x64 to ultra-high 81920x81920 cells with the latter resulting in a spatial resolution of 0.5 m if applied to a crustal scale model of 40 km. The resulting shear bands are then analyzed using the Fast Fourier Transform (FFT) to extract dominant wavelengths and power spectra, revealing features only captured at ultra-high resolution. We demonstrate the applicability of FFT analysis to a nonlinear visco-elasto-plastic setup representing a magma chamber experiencing thermal stresses during cooling and pressure deviations due to eruptions and recharge.
How to cite: Aellig, P., de Montserrat, A., Duretz, T., Kaus, B. J. P., and Räss, L.: Faults and fractures: Implications of ultra-high-resolution models for shear bands around magmatic chambers, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-3952, https://doi.org/10.5194/egusphere-egu26-3952, 2026.
The solid Earth is a complex system characterized by dynamic interactions among various tectonic components. Global mantle convection models equipped with data assimilation can effectively reproduce past subduction and associated mantle flow, providing a realistic framework for evaluating the intricate dynamic processes within the solid Earth. Despite recent advancements of data assimilation methods, their widespread application has been hindered by high computational costs due to the need for increasing model resolution and nonlinear rheology. Here, leveraging the mantle convection code ASPECT, we develop a multiscale global mantle convection model featuring adaptive data assimilation and employing nonlinear visco-plastic rheology. Our model successfully reproduces complex mantle evolution and structures, consistent with both observational constraints and previous model results. This represents the first published global mantle flow model built using ASPECT to achieve Earth-like subduction, with the aid of nonlinear rheology and adaptive data assimilation. Furthermore, the incorporation of adaptive mesh refinement and high-order finite element ensures high resolution and accuracy of model results. These advancements will contribute to a better understanding of plate tectonics and continental evolution.
How to cite: Li, X., Liu, L., and Cao, Z.: Four-dimensional multiscale global subduction models with data assimilation and realistic rheology using ASPECT, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-9034, https://doi.org/10.5194/egusphere-egu26-9034, 2026.
A major symptom of the dynamic evolution of planets is their geological map. Rock types express the imprint of the deep and surface processes. Structures often represent the tectonic styles of planets. Fossil content provides clues on past life evolution and environments. Although geologic maps are among the most used tool to infer planetary evolution, we still miss quantitative systemic links between planetary scale processes and the structure and nature of the maps. Also, the nature of rocks at the surface determines the chemistry of atmosphere-ocean composition.
The development of Earth system models that integrate geodynamic models open the possibility of generating synthetic geological mapping, a major step to link maps to dynamic evolution. Global geodynamic, climate and surface processes models provide analogous synthetic data. However, tracking these quantities through space and time and building progressively a geological map remains challenging, primarily due to the enormous quantity of information to proceed and interpret.
Here, we present an automated approach to produce geological maps from geodynamic models based on “offline” Lagrangian tracers. Tracers are introduced during post-processing and advected using a fifth-order Dormand-Prince scheme, while physical quantities are interpolated using radial basis functions (RBF). The mesh-free nature of RBF interpolation provides large flexibility in handling data of variable spatial resolution and discretization. This is particularly advantageous for coupling geodynamic models with surface processes, climate or biological evolution models, each relying on distinct spatial representation.
Like geologists combine observations to produce geological maps and interpret geodynamic systems, we integrate synthetic data to produce synthetic geological maps evolving over time. The developed software is parallel and HPC-ready, enabling efficient processing of large-scale models. Beyond geological mapping, the framework is fully generic and can interpolate and advect arbitrary fields defined on arbitrary discretizations, making it a versatile tool for multi-physics Earth system modelling.
How to cite: Jourdon, A., May, D. A., and Coltice, N.: Automatic geological mapping from geodynamic models using offline Lagrangian tracers and Radial Basis Function interpolation, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-9406, https://doi.org/10.5194/egusphere-egu26-9406, 2026.
The petrogenesis of rocks can be investigated by integrating multiple datasets and numerical methods. Petrology benefits from the abundance of field, petrographic, geochronological and geochemical data. Moreover, the growing number of numerical methods opens new opportunities for inversion. Inversion allows the quantification of useful parameters, offering deeper insights into natural processes. However, an increasing number of parameters to invert for is accompanied by significant computational costs. Therefore, relevant algorithms are needed to perform inversion and uncertainty quantification in high-dimensional spaces. Here, we demonstrate the use of Hamiltonian Monte Carlo for inverse diffusion modeling in a petrological framework. Hamiltonian Monte Carlo is a gradient-based method that efficiently explores high-dimensional parameter spaces. We implemented it using the Turing.jl package in Julia which makes use of Automatic Differentiation to efficiently explore the parameter space. Our analysis focused on the calculation of the initial cooling rate, equilibration temperature and effective mineral grain size. We fit compositional garnet data and Ar-muscovite geochronological data from the Pindos metamorphic sole (Greece) by using two different forward diffusion models. Our joint inversion shows an initial equilibration temperature of 632.3 ± 9.3 °C and a cooling rate of 202.6 ± 72.0 °C/Myr. These values reproduce not only thermobarometric and geochronological observations but also fit the garnet composition profiles and the 40Ar/39Ar age of muscovite. We finally aim to highlight the potential of Hamiltonian Monte Carlo as a robust method to perform high-dimensional inversion and constrain complex petrological processes.
How to cite: Moutzouris, D., Stroh, A., Schorn, S., and Moulas, E.: Hamiltonian Monte Carlo applied to inverse petrological problems, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-9493, https://doi.org/10.5194/egusphere-egu26-9493, 2026.
We present a modified version of optimally accurate operators for partial differential equations of arbitrary order and arbitrary dimension. Optimally accurate operators were originally proposed for seismic wave propagation in homogeneous media by Geller and Takeuchi (1995), who derived compact operator coefficients for specific wave equations in specific dimensions. Fuji and Duretz (2025) showed that these coefficients can be obtained by formulating a weak form of the PDE using basis functions defined as Taylor expansions about the centred grid point, yielding a reduction of the error by a factor of 100 for the 1D Poisson equation with three collocated grid points. However, the convergence rate varies from O4 to O2 depending on the degree of heterogeneity. Here, we generalise the theory by using Taylor expansions at all grid points involved as basis functions. The construction of symbolic expressions for the resulting coefficients requires nested loops over all grid points in space and time, which becomes intractable without GPU acceleration. In this contribution, we present the theoretical framework, benchmark results, and a Julia notebook implementing the proposed method.
How to cite: Fuji, N. and Duretz, T.: Optimally accurate operators for arbitrary PDEs using interpolated Taylor expansions, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-10044, https://doi.org/10.5194/egusphere-egu26-10044, 2026.
Geothermal reservoirs are governed by tightly coupled interactions between fluid flow, heat transport, and deformation of the host rock, posing significant challenges for numerical modeling and inversion. Capturing these processes is essential for understanding reservoir evolution, permeability development, and fluid circulation in high‐enthalpy systems. We present a forward and inverse modeling framework designed to simulate fluid migration in deforming, porous media under geothermal conditions.
The framework is implemented in the Julia programming language, exploiting its high performance and native support for automatic differentiation (AD). Forward simulations are solved on a staggered-grid, using an implicit finite-difference discretization where the jacobian matrix is assembled using AD and accelerated by automatic sparsity pattern detection. The model solves Darcy flow coupled to incompressible Stokes deformation and includes visco-elasto-viscoplastic rheology with both shear and tensile yielding, enabling the simulation of fracture-like permeability enhancement driven by thermo-mechanical stresses.
AD enables efficient adjoint-based sensitivity analysis, substantially reducing the computational cost of parameter estimation compared to traditional approaches. The resulting gradients provide the foundation for gradient-based optimization and adjoint inversions of geothermal reservoir properties. We demonstrate the capabilities of the framework through representative numerical experiments relevant to geothermal systems, illustrating its ability to capture the coupled thermal, hydraulic, and mechanical processes that control reservoir dynamics.
How to cite: Frasunkiewicz, J., Kaus, B., and Popov, A.: Forward and Adjoint Modeling of Coupled Fluid Flow, Heat Transport, and Deformation in Geothermal Reservoirs, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-10695, https://doi.org/10.5194/egusphere-egu26-10695, 2026.
PETSc is a widely used scientific computing library that allows you to write software that runs on massively parallel supercomputers. PETSc is written in C but provides interfaces for Fortran and Python. Here, I present PETSc.jl, which is a Julia interface to the (nearly) full PETSc library.
There are a number of advantages compared to other attempts:
- It can be very easily installed by typing "add PETSc" in the Julia package manager and is distributed with precompiled PETSc binaries and MPI.
- We provide both a high-level and a low-level interface. The low-level interface automatically wraps nearly the full PETSc library with over 3000 functions, whereas the high-level interface gives a more Julia-like experience but is currently limited to part of the library.
- Automatic testing and CI/CD is performed with currently >50’000 tests.
- Users can combine features from Julia, the Julia ecosystem, such as automatic differentiation and plotting with PETSc solvers.
- The resulting codes are much more compact than their counterparts in lower-level languages; yet, users still have access to all PETSc features, such as multigrid solvers for DMDA or DMStag grids.
- It allows running code on both a local workstation and on a large HPC system.
In the presentation, I will summarise some of the work done to achieve this and show scalability results of typical codes. I will also compare the timing with native compiled code.
How to cite: Kaus, B.: PETSc.jl - a julia interface to the Portable, Extensible Toolkit for Scientific Computation, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-13531, https://doi.org/10.5194/egusphere-egu26-13531, 2026.
Tailing dams constitute critical geotechnical infrastructure designed for the containment of mining waste. Its stability is governed by the complex coupling between pore fluid pressures and soil deformation [1]. Accurately modelling this behavior requires coupled porous hydro-mechanical simulations; however, the computational cost associated with high-fidelity finite element models is often unaffordable for extensive sensitivity analyses or real-time monitoring.
To overcome this bottleneck, dimensionality reduction techniques offer a robust and efficient alternative. The application of these techniques involves two phases: i) creation of a training set via a sampling of the parameter space, and ii) creation of a reduced space where to find new solutions within the family.
In this work, we explore the use of Singular Value Decomposition (SVD) and, specifically, High Order Singular Value Decomposition (HOSVD), to build multidimensional surrogates of the transient-hydromechanical models capable of accurately reproducing solutions within a fraction of the computational time.
REFERENCES
- Nasika C., P. Díez, P. Gerard, T.J. Massart and S. Zlotnik. Towards real time assessment of earthfill dams via Model Order Reduction. Finite Elements in Analysis & Design, Vol. 199, 103666, doi:10.1016/j.finel.2021.103666, 2022.
How to cite: García González, A., Zlotnik, S., and Muixí, A.: Surrogate models for transient phenomena in poromechanics, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-13625, https://doi.org/10.5194/egusphere-egu26-13625, 2026.
How to cite: Duvernay, T., Davies, R., Ghelichkhan, S., Garel, F., Scott, W., Kramer, S., Roberts, D., and Gibson, A.: Multi-material forward and inverse modelling in G-ADOPT via the conservative level-set approach, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-17242, https://doi.org/10.5194/egusphere-egu26-17242, 2026.
Earthquakes are often treated as isolated ruptures in simplified media. Growing geodetic, seismological, and geological datasets reveal deformation that extend far beyond individual events and reflect interactions between fault slip, plate motion, and mantle dynamics. Our objective is to adapt a High-Performance-Computing geodynamic code base to resolve earthquake rupture, plate tectonics forces, and mantle dynamics consistently across time scales withing a single modeling framework. Our new framework QuakeSystem.jl is being built upon an efficient and scalable matrix-free thermo-mechanical pseudo-transient framework JustRelax.jl. This requires updates to an advanced time-stepping algorithm, invariant rate-and-state dependent friction and the use of inertia. This approach allows us to explore system-scale interactions at high spatial resolution in 2D and 3D, including episodic variations in subduction rates, a flexural response of the overriding plate, and transient changes in mantle flow associated with large earthquakes. Matrix-free pseudo-transient code will significantly reduce computational cost and memory requirements while enabling massively parallel simulations on GPU architectures. This will allow us to produce many forward models efficiently, which will be incorporated in and adjoint-based inversion of megathrust slip and bulk rheological parameters from satellite-derived and ground-based deformation data, without requiring separate simulations for each parameter. By integrating observations across pre-, co-, and post-seismic periods, the framework enforces temporal and mechanical consistency and reduces non-uniqueness.
How to cite: Ibragimov, I., Koezler, A., van Amerongen, A., and van Dinther, Y.: Integrating Earthquake Cycles, Plate Tectonics and Mantle Dynamics Using Scalable GPU-accelerated Pseudo-transient Simulations, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-19510, https://doi.org/10.5194/egusphere-egu26-19510, 2026.
We are developing differentiable multi-physics solvers for extreme-scale geophysical simulations on GPUs (∂GPU4GEO project). These solvers exploit the massive parallelism of graphics processing units to significantly accelerate computations in geodynamic and related modelling applications.
As the hardware landscape evolves rapidly and GPUs are continuously enhanced with new capabilities, systematic performance benchmarking is required to better understand performance opportunities and limitations, particularly for complex workloads such as coupled Stokes and multi-physics solvers.
The Julia programming language, together with the JuliaGPU GitHub organisation, provides a unified framework for GPU programming in a high-level language. Owing to its expressive syntax and flexible compiler infrastructure, Julia enables productive development of GPU-accelerated codes without compromising performance.
We present performance benchmark results for GPUs from major vendors and assess their relative performance. We investigate benchmarks relevant for Geodynamics and ice flow codes, and consider forward and reverse passes on kernels as this is relevant for automatic gradient generation using automatic differentiation (AD). We further discuss performance per GPU price, which is relevant when evaluating the acquisition of local research computing infrastructure.
How to cite: Räss, L., Omlin, S., and Utkin, I.: Assessing high-performance GPU programming in a high-level language, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-19587, https://doi.org/10.5194/egusphere-egu26-19587, 2026.
Geodynamic numerical modelling is a common approach for investigating the mechanical evolution of large-scale tectonic features such as accretionary wedges, foreland fold-and-thrust belts, and continental rifting. Several numerical techniques have been applied in the field of geodynamics, such as the finite difference method, the finite volume method, the finite element method, and the spectral method. The long-term mechanical behaviours of large-scale tectonic processes are usually represented by the implementation of visco-elasto-plastic constitutive models to mimic the rheological mechanisms of geomaterials. However, such often nonlinear rock behaviours result in computational challenges in geodynamic numerical modeling, especially when considering the geological time scale. Here, we aim at developing and presenting an efficient numerical code (NormaJ) that acknowledges visco-elasto-plastic rheology to study large-scale tectonic processes. We apply the finite difference method with a fully staggered Eulerian grid strategy to solve nonlinear partial differential equations, integrating the marker-in-cell technique to track deformation. The code was programmed in Julia, which provides a better computational performance and data organization than Matlab and Python, while it is more user-friendly than C and Fortran. Thus, Julia is particularly suitable for large-scale geodynamics modeling. We have benchmarked our new code NormaJ with the existing fully-vectorized Matlab code Norma (https://github.com/Norma-VEP) using a baseline example, demonstrating an increase in time efficiency by a factor of about two. Such improvement indicates a good prospect for future applications and development. NormaJ will be openly available, providing the geodynamics community an easy-to-run and accessible tool to solve tectonic evolution problems.
How to cite: Wu, H. and Ruh, J.: NormaJ: An efficient numerical code for geodynamic modeling of tectonic processes, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-9897, https://doi.org/10.5194/egusphere-egu26-9897, 2026.
