The complexity in the physics of geological processes arises from their multi-physics nature, as they combine hydrological, thermal, chemical and mechanical processes. Multi-physics couplings are prone to nonlinear interactions ultimately leading to spontaneous localisation of flow and deformation. Understanding the couplings among those processes therefore requires the development of appropriate numerical tools.
Integrating high-quality data into physics-based predictive numerical simulations may lead to 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 parameters.
We invite contributions from the following two complementary themes:
#1 Computational advances associated with
- Alternative spatial and/or temporal discretisation for existing forward/inverse models
- Scalable HPC implementations of new and existing methodologies (GPUs / multi-core)
- Solver and preconditioner developments
- Combining PDEs with AI / Machine learning-based approaches (physics-informed ML)
- Automatic differentiation (AD) and differentiable programming
- Code and methodology comparisons (benchmarks)
#2 Physics advances associated with
- Development of partial differential equations to describe geological processes
- Inversion strategies and adjoint-based modelling
- Numerical model validation through comparison with observables (data)
- Scientific discovery enabled by 2D and 3D modelling
- Utilisation of coupled models to explore nonlinear interactions
The research output presented in this session can be submitted to the ongoing Special Issue (SI) in the EGU journal of Geoscientific Model Development (GMD): https://www.geoscientific-model-development.net/articles_and_preprints/scheduled_sis.html
Modeling the evolution of glaciers and ice sheets over glacial cycle timescales is critical for understanding landscape transformation through glacial erosion, predicting future changes, and assessing their impacts on sea-level rise and water availability. However, solving the partial differential equations (PDEs) governing thermomechanical ice flow at the high spatial and temporal resolutions required for these timescales is computationally prohibitive using traditional CPU-based solvers. GPU-accelerated methods offer a promising pathway to overcome these challenges.
In this study, we present a physics-informed deep learning approach leveraging GPUs, which integrates traditional numerical approximation with deep learning techniques. Using a regular grid and finite difference methods for spatial discretization, we train a Convolutional Neural Network (CNN) to minimize the energy associated with high-order ice flow equations -- a non-linear elliptic problem -- within the iterative time-stepping of a glacier evolution model. The resulting CNN, which is similar to a Variational Physics-Informed Neural Network, delivers multiple benefits: computational efficiency optimized for GPU usage, high fidelity to the original model, unsupervised training that eliminates the need for pre-generated datasets, and relatively simple implementation. Additionally, the emulator incorporates memory of prior solutions, reducing the computational cost of training -- a memory-intensive task.
Embedded within the "Instructed Glacier Model" (IGM) framework, the emulator's capabilities are demonstrated through high-resolution, large-scale simulations of glaciated landscape formation over extended timescales. This work underscores the potential of combining deep learning with physical modeling to develop scalable, efficient tools for simulating complex glaciological processes.
How to cite: Jouvet, G. and Cordonnier, G.: Emulating Viscous Ice Flow Dynamics with Physics-Informed Deep Learning, EGU General Assembly 2025, Vienna, Austria, 27 Apr–2 May 2025, EGU25-7270, https://doi.org/10.5194/egusphere-egu25-7270, 2025.
Simulating subsurface geological processes, such as compaction-driven fluid flow and rock deformation, is essential for understanding natural phenomena and addressing challenges in energy production and resource extraction. Traditional numerical methods, while effective, are computationally expensive and struggle to efficiently model large-scale dynamical problem in subsurface systems. This creates a bottleneck for large-scale simulations and real-time decision-making. Recent advances in machine learning (ML) offer promising solutions to enhance simulation efficiency. Neural operators, which learn mappings between function spaces, provide a flexible, scalable approach to modeling complex systems. Unlike traditional methods, neural operators can generalize across varying inputs and geometries, offering a more efficient and versatile alternative. This study explores the potential of advanced machine learning techniques, specifically Fourier Neural Operators (FNO) and Physics-Informed Neural Operators (PINO), to model two critical subsurface geomechanical processes: compaction-driven fluid flow and elastic stress analysis for tunnelling.
For the first case, numerical simulations were conducted to generate a dataset of up to ~10,000 samples, derived from ~1,000 different initial porosity conditions represented by randomly generated polygons. The FNO model was trained using Nvidia A100 GPUs (80G). Training loss decreased rapidly during early epochs and stabilized below 0.02 after approximately 50,000 epochs. Models trained with larger datasets (e.g., 9,753 samples) demonstrated significantly improved validation performance, achieving a validation loss of ~0.06. In contrast, models trained on smaller datasets exhibited overfitting, with validation losses exceeding 0.3. The trends in validation loss, evaluated using 60 test cases with elliptical initial conditions excluded from the training data, underscored the importance of dataset size in enhancing model generalization for machine learning-based geological simulations. The validation results demonstrated high predictive accuracy, with maximum errors below 10%. Models trained on larger datasets achieved superior performance, particularly for cases with sharper structural features. These findings highlight the capability of FNO models to effectively generalize and reproduce the dynamics of complex fluid flow in subsurface environments. The second case focuses on elastic stress analysis for tunneling, where stresses and deformations around underground excavations are critical to ensuring structural stability. Preliminary efforts have been directed toward generating numerical datasets to train FNO and PINO models, with the goal of capturing stress distribution and deformation patterns under diverse geological and engineering conditions. While results are still emerging, early indications suggest that PINO may provide additional advantages by incorporating physical laws directly into the training process, potentially reducing the amount of data required and improving computational efficiency.
This work demonstrates the transformative potential of neural operators in addressing computational challenges associated with subsurface geomechanical modeling. By combining the flexibility of data-driven methods with the robustness of physics-informed approaches, FNO and PINO offer scalable and efficient alternatives to traditional numerical methods.
How to cite: Wang, L. H. and Yarushina, V.: Harnessing Neural Network and Operators to Simulate Subsurface Geomechanical Processes , EGU General Assembly 2025, Vienna, Austria, 27 Apr–2 May 2025, EGU25-13428, https://doi.org/10.5194/egusphere-egu25-13428, 2025.
Our previous studies revealed that the impact of small-scale capillary heterogeneities is crucial for accurately predicting carbon dioxide (CO2) plume migration in geological CO2 storage sites, such as the Endurance site in the UK. While high-fidelity dynamic modelling would require excessive computational resources, conventional upscaling methods of geological models often result in dynamic models underestimating lateral CO2 migrations.
A novel capillary-limit steady-state upscaling approach is based on macroscopic invasion percolation and addresses this accuracy-feasibility trade-off. It incorporates small-scale capillary effects into upscaled local water/gas saturation functions: capillary pressure and phase permeabilities. We are developing an open-source algorithm implementation to encourage industrial adoption of the approach [1]. We present the latest advances in the library's features and performance and numerical experiments on our public set of dynamic models of real CO2 storage sites in the North Sea [2].
The latest library advances include substantial parallelisation, single-core optimisations, an optional hydrostatic term, support for anisotropic fine-scale permeability, a stochastic re-upscaling approach for porosity and permeability fields upscaled by averaging, library infrastructure, and more. Numerical experiments aim to assess the impact of upscaling under uncertainties in rock and multiphase flow properties. We also attempt to downscale CO2 plumes simulated at coarse scales using data from the algorithm's percolation step, providing an estimate for fine-scale dynamics.
[1] https://github.com/ImperialCollegeLondon/StrataTrapper
[2] https://github.com/ImperialCollegeLondon/StrataTrapper-models
How to cite: Elizarev, M., Krevor, S., and Muggeridge, A.: Capillary heterogeneity upscaling using macroscopic percolation: code advances and field-scale dynamic CO2 storage simulations, EGU General Assembly 2025, Vienna, Austria, 27 Apr–2 May 2025, EGU25-4162, https://doi.org/10.5194/egusphere-egu25-4162, 2025.
The flow of fluids within porous rocks is an important process with numerous applications in Earth sciences. Modeling the compaction-driven fluid flow requires the solution of coupled nonlinear partial differential equations that account for the fluid flow and the solid deformation within the porous medium. Despite the nonlinear relation of porosity and permeability that is commonly encountered, natural data show evidence of channelized fluid flow in rocks that have an overall layered structure. Layers of different rock types routinely have discontinuous hydraulic and mechanical properties.
We present numerical results [1] obtained by a novel space-time method [2] based on a fixed-point scheme inspired by the mathematical analysis [3], combined with a space-time least-squares formulation. This approach can handle discontinuous initial porosity (and hence permeability) distributions. It furthermore exhibits optimal convergence independently of the discontinuities, while standard approximations, as e.g. finite differences, tend to show lower order convergence in discontinuous regimes.
The space-time method enables a straightforward coupling to models of mass transport for trace elements. Our results show the influence of different kinds of layering in the development of fluid-rich channels and mass transport [1].
References
[1] Fluid flow channeling and mass transport with discontinuous porosity distribution, S. Boisserée, E. Moulas and M. Bachmayr, arXiv Preprint (2024), https://doi.org/10.48550/arXiv.2411.14211.
[2] An adaptive space-time method for nonlinear poroviscoelastic flows with discontinuous porosities, M. Bachmayr and S. Boisserée, arXiv Preprint (2024), https://doi.org/10.48550/arXiv.2409.13420.
[3] Analysis of nonlinear poroviscoelastic flows with discontinuous porosities, M. Bachmayr, S. Boisserée and L. M. Kreusser, Nonlinearity (2023), https://doi.org/10.1088/1361-6544/ad0871.
How to cite: Boisserée, S., Moulas, E., and Bachmayr, M.: Fluid flow channeling and mass transport with discontinuous porosity distribution, EGU General Assembly 2025, Vienna, Austria, 27 Apr–2 May 2025, EGU25-1710, https://doi.org/10.5194/egusphere-egu25-1710, 2025.
Compositional concentration profiles across individual crystals or diffusion couples are largely determined by diffusion and growth processes. These two processes are particularly important during the formation of high-temperature rocks such as igneous and metamorphic rocks. The numerical simulation of concentration profiles in crystals is a widely used technique in various fields such as geospeedometry or diffusion chronometry. Compared to single crystals, coupled diffusion pairs yield tighter constraints on the experienced temperature and pressure ranges and thus provide additional information for our models. However, the numerical description of concentration profiles within diffusion couples is challenging due to the sharp compositional gradients. Discontinuities in concentration, which are related to the different mineral properties, commonly occur at the interface of two minerals and lead to technical implementation issues.
To address these issues, we have developed a Finite Element (FE) package in Julia that can calculate the evolution of concentration profiles in diffusion couples with moving interfaces. Both growth and diffusion processes are considered. An adaptive grid enables the accurate reproduction of rapid concentration changes and discontinuities. Our code can be applied to various examples of single crystals or diffusion couples, integrating any combination of growth, diffusion, and temperature dependency. Additionally, it is possible to calculate concentration profiles based on the thermodynamically-constrained, Stefan-Interface condition. Results from our models can be used in petrology and geodynamic applications to provide tighter constraints concerning in the pressure and temperature evolution of magmatic and metamorphic mineral assemblages.
How to cite: Stroh, A., Aellig, P., and Moulas, E.: Numerical modeling of simultaneous diffusion and mineral growth, EGU General Assembly 2025, Vienna, Austria, 27 Apr–2 May 2025, EGU25-4080, https://doi.org/10.5194/egusphere-egu25-4080, 2025.
Finely resolved Geodynamic mantle convection models are crucial to understand
the detailed physics governing major geological process and to infer
the mineralogical state of the Earth. Given the scale needed to globally
resolve features of the mantle such as unstable boundary layers
and asthenospheric flows where viscosity can change by around four orders of
magnitude within 50 km, traditional sparse matrix methods become unsuitable
due to their immense memory requirements. TerraNeo builds upon the
massively parallel matrix-free finite element framework HyTeG which uses
geometric multigrid solvers on block-structured hybrid tetrahedral grids.
The extreme scalability of HyTeG has been demonstrated previously, solving
Stokes problems with trillions (1e12) of unknowns, corresponding to ≃ 1
km global resolution of Earth’s mantle.
We discuss the features and scalability of TerraNeo with respect to the
numerical treatment of the truncated anelastic liquid approximation as a
model for mantle convection. The compute kernels are evaluated and verified
through geophysical applications, convergence studies, and community
benchmarks covering sharp viscosity variations, nonlinear rheologies and
mixed boundary conditions.
How to cite: Ilangovan, P., Kohl, N., Mohr, M., Brown, H., D'Ascoli, E., Papanagnou, I., Vilacis, B., and Bunge, H.-P.: Extreme-Scale Geodynamic Modelling with TerraNeo, EGU General Assembly 2025, Vienna, Austria, 27 Apr–2 May 2025, EGU25-14096, https://doi.org/10.5194/egusphere-egu25-14096, 2025.
The behavior of interfaces in seismic faults, and tribological systems in general, is governed by the interaction of discrete microconstituents trapped between contacting surfaces, often referred to as the gouge or the third-body. This is an amorphous wear particle agglomeration undergoing significant deformation, while the surrounding regions experience comparatively moderate strain. Understanding the dynamics of such systems, and in particular the evolution of the third body, can rely on particle-based numerical models such as the Discrete Element Method (DEM). The predicted behavior of the gouge can be sensitive to the boundary conditions, and therefore to the system size. However, the important computational costs prevent handling arbitrarily large domain sizes, which calls for cheaper Discrete-Continuum coupled methods. To accurately capture the behavior of continuum (long-range boundary) and discrete regions (gouge), we employed a hybrid modeling strategy combining the Finite Element Method (FEM) for continuum regions and the Discrete Element Method (DEM) [1,2].
To accommodate the evolving nature of the third body and prevent limitations imposed by the size of the discrete region, we will introduce in this presentation an adaptive coupling. This approach allows FEM regions to transition dynamically into DEM regions when a sufficient deformation criterion is met. Such a condition is evaluated within the third body near the coupling region. The adaptive framework supports large-scale simulations, and it will be demonstrated to support amorphous and ordered (crystalline) material setups for a gouge. Finally, the adaptive coupling is used to model the evolution of a third body comprising elliptical rigid bodies, which will be shown to impact the gouge evolution in certain conditions. Our findings underscore the importance of coupling techniques in modeling the complex, multiscale nature of frictional interfaces and contribute to a deeper understanding of the role of granularity in dynamic friction and third-body evolution.
[1] Xiao, S. P. and T. Belytschko (2004). “A bridging domain method for coupling continua with molecular dynamics”. Computer Methods in Applied Mechanics and Engineering. doi: 10.1016/j.cma.2003.12.053
[2] Voisin-Leprince, M., J. Garcia-Suarez, G. Anciaux, and J.-F. Molinari (2022). “Finite element method–discrete element method bridging coupling for the modeling of gouge”. International Journal for Numerical Methods in Engineering. doi: 10.1002/nme.7171
How to cite: Anciaux, G., Voisin-Leprince, M., and Molinari, J.-F.: Adaptive FEM-DEM bridging coupling to study third-body/gouge evolution, EGU General Assembly 2025, Vienna, Austria, 27 Apr–2 May 2025, EGU25-15996, https://doi.org/10.5194/egusphere-egu25-15996, 2025.
Modern geodynamic models have become increasingly complex, coupling detailed numerical approximations to many physical processes with large observational datasets. This coupling creates unique challenges for modern research software such as how to combine complex rheologies utilizing multiple flow mechanisms with the simultaneous modeling of mineral microstructure; how to model realistic geometries and evolving surface topography while simultaneously including large observational datasets of topography and surface deformation; and how to utilize highly-optimized and scalable numerical solvers while keeping up with changing high-performance computing architectures.
We here present our approach to reconciling these challenges: The next major release of ASPECT - The Advanced Solver for Planetary Evolution, Convection, and Tectonics. During the six years since our last major release, we have implemented many new features and improvements. Here we report on a new major release that highlights ASPECT's increased flexibility in modeling complex tectonic and convection problems. New features we will present at the workshop are in particular:
- A new default Stokes solver utilizing a matrix-free geometric multigrid preconditioner
- Complex rheologies like visco-elasto-plasticity including Peierls-, dislocation-, and diffusion-creep
- Models of pinned grain-size evolution in a two-mineral assemblage
- Evolution of crystal-preferred orientation using DREX like algorithms
- Utilizing modern external libraries for the accurate solution of ordinary differential equations
- Extended support for efficiently including large-scale datasets in parallel models
- Interfaces to surface evolution modeling software like Fastscape and others
- Optimizing finite element type, degree, and advection method for different compositions
- Major improvements to the structure of the code base, plugin systems, and user interface
As usual the release is open-source and freely available at https://aspect.geodynamics.org/. We hope that providing well-documented, flexible, and tested geodynamic research software provides the community with the necessary tools to tackle the geodynamic research questions of the next decade.
How to cite: Gassmöller, R., Heister, T., Bangerth, W., Dannberg, J., Fraters, M., Glerum, A., Myhill, R., and Naliboff, J.: ASPECT 3.0: The Advanced Solver for Planetary Evolution, Convection, and Tectonics, EGU General Assembly 2025, Vienna, Austria, 27 Apr–2 May 2025, EGU25-16478, https://doi.org/10.5194/egusphere-egu25-16478, 2025.
The movement of the Earth's crust and mantle in geodynamics is typically modeled as the flow of a viscous fluid governed by the Stokes equations. Incorporating plasticity into material rheology often results in mesh-dependent behavior, which poses challenges for accurate numerical simulations. Several approaches have been proposed to mitigate mesh dependence and develop solvers that decouple solution errors from viscosity and mesh size.
Traditionally, Taylor-Hood (TH) and Crouzeix-Raviart (CR) elements of order 2 are used for geodynamics simulations. In this study, we examine the numerical solution of variable-viscosity Stokes flow with plasticity and Drucker-Prager type yielding using Scott-Vogelius (SV) compatible finite elements in combination with pseudo-Jacobian and augmented Lagrange methods. The Scott-Vogelius element is unique among finite elements for the mixed formulation of Stokes flow, as it has an associated De Rham complex. This theoretically ensures a divergence-free velocity field. We investigate the degree of decoupling between velocity errors, pressure errors, and viscosity-induced errors in a viscoelasto-plastic case study.
Our results show that Taylor-Hood elements (CG2 × CG1 for velocity and pressure) fail to provide accurate solutions in such cases. While the low-order CR elements perform better, the higher-order SV elements (CG4 × DG3) yield the best results.
We conclude that due to the inherent mesh-dependent behavior and viscosity dependent errors in TH elements, CR or SV elements should be preferred for geodynamics simulations.
How to cite: Norouzi-Moghanjoghi, K., García-Pintado, J., and Perez-Gussinye, M.: On the Application of Compatible Finite Elements for Divergence-Free and Mesh-Independent Viscoelastic-Plastic Rheology in Geodynamics Simulations, EGU General Assembly 2025, Vienna, Austria, 27 Apr–2 May 2025, EGU25-17228, https://doi.org/10.5194/egusphere-egu25-17228, 2025.
The transition to exascale (>1000 Petaflops) computing necessitates the adaptation of numerical modeling tools to efficiently utilize emerging high-performance computing architectures. Within the ChEESE-2P project, the further development of LaMEM (Lithosphere and Mantle Evolution Model) focuses on achieving scalable performance on advanced systems, including the EuroHPC supercomputer LUMI (currently #3 in Europe). Taking advantage of using the PETSc library, LaMEM demonstrates strong and weak scalability, achieving linear performance up to 512 compute nodes and supporting high-resolution simulations with grids up to 10243
In response to the increasing emphasis on GPU-based computing, ongoing efforts are directed towards optimizing LaMEM for GPU architectures, including both NVIDIA and AMD systems. Preliminary results highlight significant progress in enabling GPU-accelerated runs and improving resource utilization. This work highlights LaMEM's ability to perform large-scale geodynamic simulations, contributing to the broader goal of integrating physics-based models with available data.
How to cite: Ibragimov, I., Kaus, B., and Popov, A.: Scaling staggered grid code on pre-exascale machines, EGU General Assembly 2025, Vienna, Austria, 27 Apr–2 May 2025, EGU25-9151, https://doi.org/10.5194/egusphere-egu25-9151, 2025.
From complex magmatic systems to geothermal reservoirs, fluid-rock dynamics have posed immense modeling hurdles in the Geosciences for decades. These systems are influenced by interactions between magmatic heat sources, fluid flow, host rock deformation, and chemical heterogeneities. To enhance our understanding and predictive capabilities of these intricate systems, we develop a novel forward and inverse modeling code designed to simulate fluid migration within a deforming, porous host-rock. We employ the Julia programming language, chosen for its differentiability and efficiency, facilitating the simple integration of various composable packages.
Forward simulations are performed in the advanced automatic differentiation (AD) framework of Julia, allowing for flexible adjustments of the underlying coupled system of equations. We utilize a staggered-grid, implicit finite-difference solver, along with the GeoParams.jl package to implement visco-elasto-plastic rheologies and solve the coupled fluid-rock interactions under non-linear Darcy and incompressible Stokes-flow regimes.
The AD framework of Julia allows for the application of the adjoint method in parameter sensitivity analysis, significantly reducing computational demands compared to traditional inversion techniques. This framework lays the foundation for adjoint inversions as the gradients calculated for the sensitivities are needed for gradient descent algorithms. The effectiveness of our framework is demonstrated through representative case studies, illustrating its applicability to understanding the dynamic behavior of two-phase systems influenced by both thermal and mechanical processes.
How to cite: Frasunkiewicz, J., Kaus, B., Popov, A., Schuler, C., and Riel, N.: Leveraging Differentiable Programming in Julia: Forward and Inversion Modeling of Two-Phase Systems, EGU General Assembly 2025, Vienna, Austria, 27 Apr–2 May 2025, EGU25-10269, https://doi.org/10.5194/egusphere-egu25-10269, 2025.
The adjoint method for the Stokes equations provides a versatile and highly efficient approach to investigate the underlying physics of geodynamic processes. Reuber et al. (2018) demonstrated that adjoint sensitivities can be used to develop scaling laws for processes like folding and subduction dynamics. The gradients derived using the adjoint method can also directly be used in inversions in geodynamic applications. However, previous implementations of the adjoint method have typically been highly problem-dependent and often limited to viscous rheologies. Extending it to other nonlinear rheologies typically required substantial additional work, which is likely one of the reasons that the method has not yet been widely adopted in solid Earth geosciences.
To overcome this problem, we use automatic differentiation (AD) to compute the gradients needed to develop an adjoint solver for the Stokes equations. The gradients are computed using the Julia package Enzyme.jl. The adjoint solver is designed to be problem-agnostic, where the gradients are automatically computed for any user-defined rheology, from a simple linear viscous model to a complex visco-elasto-viscoplastic composite rheology. This functionality is added to the JustRelax.jl thermo-mechanical solver, where we use the same pseudo-transient solver strategy to solve both the forward and adjoint problems. This approach ensures that the adjoint solver remains consistent and fully generic.
The method is applied to analyse horizontal plate motion around subduction zones. For different material parameters, it is possible to calculate sensitivity kernels that show, for each location in the numerical domain, how much these parameters influence the horizontal plate motion (e.g. Reuber et al (2020)). The scaling of sensitivities for different parameters is discussed to enable a quantitative comparison. This approach is then used to identify the most influential factors affecting plate motion.
Reuber, G. S., Popov, A. A., & Kaus, B. J. (2018). Deriving scaling laws in geodynamics using adjoint gradients. Tectonophysics, 746, 352-363.
Reuber, G. S., Holbach, L., Popov, A. A., Hanke, M., & Kaus, B. J. (2020). Inferring rheology and geometry of subsurface structures by adjoint-based inversion of principal stress directions. Geophysical Journal International, 223(2), 851-861.
How to cite: Schuler, C., de Montserrat, A., Frasunkiewicz, J., Aellig, P., and Kaus, B.: Parameter Sensitivity Analysis of Plate Motion using the Adjoint Method and Automatic Differentiation, EGU General Assembly 2025, Vienna, Austria, 27 Apr–2 May 2025, EGU25-10823, https://doi.org/10.5194/egusphere-egu25-10823, 2025.
Pegmatites and rare metal granites are granitic igneous rocks distinguished by their texture, which is dominated by crystal growth. These rocks are frequently enriched in rare elements (e.g., Li, Cs, Be, Nb, Ta) and represent economically significant deposits, classified among the critical raw materials identified by the European Commission. Our objective is to better constrain the tectonic parameters that govern the emplacement of pegmatites within the continental crust, including their migration rates and durations, with a particular emphasis on the role of temperature in these crustal migration processes.
To model those fluid migrations, two-phase flow in Julia, based on porosity waves with compressible fluid is used. The porosity is interpreted as melt [1]. To improve the yet existing codes [2], we implement temperature in our two phase flow formulation from energy conservation. Temperature allow a thermomechanical coupling with rock viscosity. A first equation with a simple exponential coupling is used to understand thermal implication on viscosities.
The model represents a continental crust with partial melting occurring within the lowermost 5 km, where temperature is maintained at 750°C due to underplating. A constant geothermal gradient is applied from this depth to the surface. A 10% porosity anomaly is introduced in the partially molten zone, while a baseline porosity of 1% is applied throughout the model to ensure numerical stability. The fluid viscosity is set at 10^4 Pa.s, while at 750°C, the rock viscosity is 10^16 Pa.s. A constant permeability of 10^-11 m² is applied throughout the model. Thermomechanical couplings of varying strength are implemented to assess their impact on migration processes. Accordingly, the rock viscosity at 450°C is varied between 10^16Pa.s and 10^21 Pa.s.
Models reveals two distinct mechanisms that halt migration. The first occurs when the thermomechanical coupling is low (soft and hot crust). Allowing rock viscosity to remain low, so melt migration can outpace thermal diffusivity. Enabling the melt to be in an undercooling state. This means that the magma can migrate beyond the point at which the surrounding rock reaches the crystallization temperature of the melt, a necessary condition for pegmatite formation. The second case arises when thermomechanical coupling is strong, causing the surrounding rock's viscosity to become too high for the magma to reach undercooling condition. In this scenario, the magma crystallizes as soon as it reaches the surrounding rock at its crystallization temperature, potentially becoming trapped by a viscous layer (of different or colder nature).
The use of a geothermal gradient more representative of a metamorphic core complex, along with an improved thermomechanical coupling, should refine the estimates of migration time and distance. Similarly, the introduction of viscous heterogeneity would help highlight the geological structures that may or may not facilitate the migration of these magmas to shallower levels of the crust.
References:
[1] L. Räss, T. Duretz & Y.Y. Podladchikov (2020). https://doi.org/10.1093/gji/ggz239
[2] A. Plunder & al. (2022). https://doi.org/10.1016/j.lithos.2022.106652
How to cite: Bergogne, M., Le Pourhiet, L., Räss, L., and Plunder, A.: 1D modelling of pegmatite migration, EGU General Assembly 2025, Vienna, Austria, 27 Apr–2 May 2025, EGU25-11537, https://doi.org/10.5194/egusphere-egu25-11537, 2025.
The staggered grid finite difference is a robust discretization method for the high-resolution 3D geodynamic simulations that involve heterogeneous material parameters. Achieving its scalability on parallel machines inevitably requires the application of multigrid solvers. In particular, a coupled velocity-pressure geometric multigrid preconditioner based on Galerkin coarsening scheme demonstrates very good results. However, this method relies on assembled matrices which have a significant memory imprint and prohibits achieving peak performance due to suboptimal use of the limited memory bandwidth. A geometric multigrid method, based on the re-discretization of linear operators on the coarser levels, converges generally slower, but can be implemented in a completely matrix-free manner. It poses a valuable alternative to the Galerkin method, since an increased number of iterations can be compensated by a greater performance of the matrix-vector products computed on the fly without storing matrices in the memory.
Here we present a hybrid approach that allows optimal combination between various types of coarsening techniques for staggered grid discretizations. This work is performed within the framework of ChEESE-2p project (Centre of Excellence for Exascale in Solid Earth) and involves the flagship code LaMEM (Lithosphere and Mantle Evolution Model), which is based on Portable Extensible Toolkit for Scientific computation (PETSc), following an approach suggested for finite element discretizations by May et al. (2015). Here, we extend it to the staggered grid finite difference, discuss the optimal solver parameter selection, and document performance gains that can be achieved by using the matrix-free operators.
We typically start with a few levels of re-discretized matrix-free operators, followed by Galerkin geometric coarsening operating on assembled matrices. This approach ensures that most of the optimization and memory saving is already obtained at the top levels, whereas more robust Galerkin coarsening can be used at coarser levels without compromising the convergence. At the coarse grid level, we either utilize a parallel sparse direct solver or a black-box algebraic multigrid method. The number of processors participating in a coarse grid solve can be optimally selected via PETSc sub-communicator framework (Telescope).
D. A. May, J. Brown, L. Le Pourhiet, 2015. A scalable matrix-free multigrid preconditioner for finite element discretizations of heterogeneous Stokes flow, Comput. Methods Appl. Mech. Engrg., 290, 496–523.
How to cite: Popov, A. A., Kaus, B. J. P., and Ibragimov, I.: Scalable hybrid multigrid for staggered grid discretizations in geodynamics, EGU General Assembly 2025, Vienna, Austria, 27 Apr–2 May 2025, EGU25-11825, https://doi.org/10.5194/egusphere-egu25-11825, 2025.
In geodynamic models, the incorporation of a free surface boundary condition is crucial to better understand and resolve, for example, the coupled interaction between the lithosphere, deep mantle, and surface processes. A free surface is similarly important in ice flow models to capture the geometry of the ice sheet and its temporal evolution.
While the implementation of a free surface is relatively straightforward in finite element methods due to their capacity to manage complex and deformable geometries, and boundary conditions, the treatment of a free surface boundary (or any other boundary) that is not aligned with the staggered grid of a finite difference (FD) scheme poses a significant challenge. A common approach adopted by FD geodynamic codes, such as I2/3VIS (Gerya and Yuen, 2007) or LAMEM (Popov and Kaus, 2016), involves the incorporation of an additional rheology layer above the surface, simulating the presence of air. However, the of Stokes solvers is constrained by the viscosity contrast occurring within the domain, typically in the range of 6 to 7 orders of magnitude. Consequently, the viscosity of the air layer is limited to values within the 1e16-1e18 Pa*s range. This approach is not only physically unrealistic, but also leads to an inaccurate topography evolution, and introduces a very strong and sharp viscosity contrast at the rock and air interface, which hinders the convergence of iterative solvers particularly the Accelerated Pseudo-Transient (APT) method employed in this study.
To address these limitations, we propose the implementation of a variational Stokes approach (Larionov et al. 2017), which allows for the incorporation of both real free surface and solid wall boundary conditions. This approach is then combined with either a marker chain or a level set to track evolution of the surface. We demonstrate that this approach greatly improves the convergence rate of the iterative APT solver, as well as demonstrate its accuracy and applicability to geodynamic and ice flow simulations with a set of benchmarks and toy codes.
Gerya, Taras V., and David A. Yuen. "Robust characteristics method for modelling multiphase visco-elasto-plastic thermo-mechanical problems." Physics of the Earth and Planetary Interiors 163.1-4 (2007): 83-105.
Popov, Anton, and Boris Kaus. "3D modelling of non-linear visco-elasto-plastic crustal and lithospheric processes using LaMEM." EGU General Assembly Conference Abstracts. 2016.
Larionov, Egor, Christopher Batty, and Robert Bridson. "Variational stokes: A unified pressure-viscosity solver for accurate viscous liquids." ACM Transactions on Graphics (TOG) 36.4 (2017): 1-11.
How to cite: de Montserrat Navarro, A., Aellig, P., Gray, T., Räss, L., and Utkin, I.: Variational Stokes: free surface for staggered grid finite differences schemes for geodynamic and ice flow modelling, EGU General Assembly 2025, Vienna, Austria, 27 Apr–2 May 2025, EGU25-16111, https://doi.org/10.5194/egusphere-egu25-16111, 2025.
Although geodynamics and ice flow dynamics address distinct physical systems, they share significant computational and modelling challenges. Both require vast, data-intensive simulations on next-generation high-performance computing (HPC) platforms. With limited observational data, these models must be rigorously constrained to improve their predictive power. Our work focuses on differentiable modelling of Earth’s largest ice sheets and high-resolution 3D geodynamic processes, such as magmatic systems and the formation of the Alps.
We are developing differentiable multi-physics solvers for extreme-scale geophysical simulations on GPUs - ∂GPU4GEO. These high-performance, scalable tools leverage advanced programming techniques, particularly automatic differentiation (AD) within the Julia programming language. Using Enzyme.jl, an AD tool integrated with the LLVM compiler, we combine differentiation with compiler optimisations. This approach enables highly efficient reverse-mode AD, achieving near-theoretical peak performance.
Building on the GPU4GEO PASC project (2020–2024), we are extending pseudo-transient solvers with differentiable modelling capabilities. The modular GPU4GEO software stack, composed of specialised Julia packages, provides solvers for diverse physical systems and customisable building blocks. By integrating Enzyme.jl into the entire stack, we enable high-performance AD on GPUs while maintaining support for distributed-memory parallelism via MPI. These developments ensure scalability on flagship supercomputers and facilitate efficient exploration of geophysical processes.
This collaborative effort targets applications requiring large-scale simulations to address critical scientific challenges. The resulting computational tools are optimised for next-generation GPU architectures, offering transformative potential for geodynamics and glaciology research.
How to cite: Räss, L., Utkin, I., De Montserrat, A., Kaus, B., Tackley, P., Moses, W., and Duretz, T.: Differentiable multi-physics solvers for extreme-scale geophysics simulations on GPUs, EGU General Assembly 2025, Vienna, Austria, 27 Apr–2 May 2025, EGU25-16577, https://doi.org/10.5194/egusphere-egu25-16577, 2025.
Reconstructing the evolution of a landscape provides insights into its geological and/or climatic history and into processes shaping the earth's surface: what was the configuration of the drainage network before a specific geological or climatic event, what are the areas that are currently most sensitive to fluvial incision or hillslope processes, or to which extent lithological contrasts influenced landscape evolution are frequent questions. Most of the time, these questions are addressed with forward models in which only a small part of the parameter space is explored.
Landscape evolution can be simulated using a diffusion-advection equation where the diffusive term represents hillslope erosion-deposition processes and advection simulates river incision. In this case, the advection velocity can be calculated from drainage area and erodibility parameters of the Stream Power Law. The model can also include a source term, which simulates tectonic uplift. This approach permits to solve a PDE and formulate an adjoint model that can be used for parameter inversion and sensitivity analysis. In this study, we use the Firedrake package which includes automatic differentiation procedures for building the adjoint model. Our results illustrate different model sensitivities to diffusion and erodibility coefficients, and show that it is able to reconstruct spatial variations of these coefficients.
We then apply the adjoint model to sensitivity analysis and to parameter inversion in real-world cases. The first one is the southeastern border of the French Massif Central, for which we seek at understanding what was the topography prior to a major incision by tributaries of the Rhone River. The second case is the footwall topography along a segment of the Wasatch normal fault, USA, for which we aim at quantifying temporal uplift rate variations.
In the French Massif Central, inversion of the initial conditions reveals that the pre-incision topography consists of a relatively smooth and flat footwall delimited by a well-defined and linear fault escarpment that corresponds to a Mesozoic normal fault system. In the Wasatch range, the model indicates a significant increase in the uplift rate of the Wasatch Range, from 0.2 to 1 mm.yr-1, since approximately 2 Ma, aligning well with recently published estimates.
How to cite: Petit, C., Jourdon, A., and Coltice, N.: Reconstructing landscapes: an adjoint model of the stream power and diffusion equation., EGU General Assembly 2025, Vienna, Austria, 27 Apr–2 May 2025, EGU25-17495, https://doi.org/10.5194/egusphere-egu25-17495, 2025.
Glacier evolution models (GEMs) are important tools in glaciology to predict future glacier response from climate forcing. However, reconstructing past climates requires inversion tools to infer the climate forcing that explains paleoglacier extents documented through historical records or geomorphic features.
Such "inverted" GEMs are less common compared to forward GEMs but are important to better constrain climate from the past.
For instance, Visnjevic et. al proposed a model to estimate the Equilibrium Line Altitude (ELA) from reconstructed paleoglacier extents.
However, their approach assumes stationary glaciers, neglecting temporal dynamics, and employs a heuristic inversion technique.
Recent implementations of automatic differentiation (AD), coupled with Graphic Processing Unit (GPU) performance improvements, provide a promising pathway to develop fully differentiable and computationally efficient GEMs. Here, we introduce an Invertible Glacier Evolution Model (IGEM), a new framework designed to overcome the limitations of existing inversion methods.
Our IGEM relies on the Shallow Ice Approximation (SIA) for the ice flow, and surface mass balance is computed with the Positive Degree Day (PDD) .
The model’s tensor-based architecture leverages GPU acceleration and enables efficient computation of gradients with respect to input climate variables, such as temperature and precipitation, which are used for PDD calculations. The gradient-descent inversion scheme employed in our IGEM converges more rapidly, delivers more accurate solutions, and offers greater generality (e.g., it is not constrained by the stationary assumption) compared to heuristic inversion methods.
The main challenge here is due to the fact that one forward GEM simulation requires thousands of iterations to model a glaciation spell.
To compute the gradient of the cost function with respect to climate forcing, a chain derivation of all operations within the forward GEM is necessary, which is memory-challenging, especially on GPUs.
To address this, our IGEM selectively recomputes a subset of intermediate operations during gradient computation. Instead of storing all operations, only those essential for computing gradients are cached, while others are recomputed during the "backward" pass. This approach reduces memory usage at the cost of increased computation time, enabling the methodology to handle large-scale problems effectively.
To illustrate the feasibility of our approach, we apply it to the inference of climate proxies at the Aletsch Glacier for the period 1880–2010. We leverage sequentially dated observations of the glacier geometry during the same timeframe. Given the nonuniqueness of the problem, the method permits the derivation of a set of compatible temperature and precipitation proxies, which are evaluated against weather station data near the glacier.
This proof-of-concept shows that our IGEM approach enables the extraction of compatible climate proxies, such as temperature and precipitation, provided documented glacier former extents. By bridging the gap between glacier dynamics and climate reconstruction, our IGEM has the potential to advance our understanding of past climates in formerly glaciated regions.
How to cite: Lleshi, K., Jouvet, G., and Herman, F.: Retrieving climate proxies from an invertible glacier evolution model., EGU General Assembly 2025, Vienna, Austria, 27 Apr–2 May 2025, EGU25-17922, https://doi.org/10.5194/egusphere-egu25-17922, 2025.
Geodynamic modeling has become a crucial tool for investigating the dynamics of Earth deformation across various scales. This approach involves solving mechanical problems characterized by significant property variations (e.g., viscosity, shear modulus, conductivity) under nearly incompressible conditions. Recent advancements in technology have facilitated the development of iterative pseudo-transient solvers, which require minimal global communication and enable near-optimal parallel scaling on supercomputers. However, selecting numerical parameters that ensure both robust and rapid convergence remains a challenging task. In this contribution, we explore potential strategies to address these challenges and provide application examples using both finite difference and face-centered finite volume methods.
How to cite: Duretz, T., de Montserrat, A., Sevilla, R., Räss, L., and Utkin, I.: Automatic tuning of iterative pseudo-transient solvers for geodynamic modelling, EGU General Assembly 2025, Vienna, Austria, 27 Apr–2 May 2025, EGU25-19536, https://doi.org/10.5194/egusphere-egu25-19536, 2025.
-
In this contribution, we generalise the optimally accurate operators proposed and used in the series of studies on the simulation of seismic wave propagation, especially based on Geller & Takeuchi (1995). Although the operators have been mathematically proven more accurate than conventional methods, the demonstration has been made without a recipe ready for different configurations and the theory is complicated using normal-mode theory, which prevents other physics from applying the methods. Here we show that the operators can be systematically obtained for any form of partial differential equations and we show several applications with numerical examples.
How to cite: Fuji, N. and Duretz, T.: Optimally accurate operators for partial differential equations, EGU General Assembly 2025, Vienna, Austria, 27 Apr–2 May 2025, EGU25-2822, https://doi.org/10.5194/egusphere-egu25-2822, 2025.
In this work we describe dimensionality reduction methods to solve parametric problems governed by partial differential equations (PDE). The goal of these methodologies is to elucidate and reduce the dimensionality of the manifold containing the family of solutions of a parametric problem. This leads to a reduced system that can often be solved in real time. These techniques are being successfully applied in many fields in science and engineering, e.g. [1,2] as well as in geodynamics [3,4].
The application of these techniques involves two phases: i) creation of a reduced space, often done via a sampling of the parametric space and a singular value decomposition, and ii) use of the reduced space to find a new solution within the family.
Here we want to compare methodologies that, in their second step, include or neglect the physics described by the PDEs. We call these, physics-based and physics-agnostic approaches. Reduced Basis methods being an example of the first, and surrogate modelling one of the second.
Applications of flow in porous media are used as examples to test the strengths and weaknesses of the different approaches. These methodologies are a potential tool to be used in situations where it is unaffordable to obtain a very large training set (big data).
REFERENCES
[1] Rocas M., A. García-González, S. Zlotnik, X. Larráyoz and P. Díez. Nonintrusive Uncertainty Quantification for automotive crash problems with VPS/Pamcrash. Finite Elements in Analysis & Design, Vol. 193, doi:10.1016/j.finel.2021.103556, 2021.
[2] Muixí A., S. Zlotnik, P. Calvet, M. Espanol, I. Lodoso-Torrecilla, M.P. Ginebra, P. Díez and A. García-González. A multiparametric advection-diffusion reduced-order model for molecular transport in scaffolds for osteoinduction. Biomechanics and Modeling in Mechanobiology, doi:10.1007/s10237-022-01577-2, 2022.
[3] Ortega O., S. Zlotnik, J.C. Afonso and P. Díez. Fast Stokes flow simulations for geophysical-geodynamic inverse problems and sensitivity analyses based on reduced order modeling. Journal of Geophysical Research: Solid Earth, Vol. 125, 1–25, doi:10.1029/2019JB018314, 2020.
[4] Manassero M., J.C. Afonso, F. Zyserman, A. Jones, S. Zlotnik and I. Fomin. A Reduced Order Approach for Probabilistic Inversions of 3D Magnetotelluric Data II: Joint inver- sion of MT and Surface-Wave Data. Journal of Geophysical Research - Solid Earth, doi:10.1029/2021JB021962, 2021.
How to cite: García-González, A. and Zlotnik, S.: Physics-based and physics-agnostic reduced order modeling, EGU General Assembly 2025, Vienna, Austria, 27 Apr–2 May 2025, EGU25-11462, https://doi.org/10.5194/egusphere-egu25-11462, 2025.
I present a set of numerical experiments, which show the formation, evolution, and influence of small-scale rolls at the lithosphere-asthenosphere boundary (LAB). The rolls originate from lateral gradients of temperature and are not related to classical large-scale Rayleigh–Bénard convection. Rather, they represent a form of horizontal convection, arising from the circumstance that the lower part of the lithosphere (both oceanic and continental) can contribute to the advection of material (due to a relatively low viscosity) but is characterized by a non-adiabatic thermal regime. The formation of convection rolls indicates that the process is relatively steady, with a relaxation time of several hundreds kyrs. Although the LAB geometry influences the formation of convective cells, these features form even in presence of a flat LAB surface, whenever there is a lateral thermal change within the lithosphere, for example at the transition between oceanic and continental lithosphere along continental margins. An important observation is that the thermal structure of the oceanic lithosphere close to a spreading center induces an ascending flow even in absence of extension. Consequently, an active component of spreading exists regardless of whether two plates are moving apart. In these experiments, the active component of spreading induces a velocity between 0.6 and 1.2 cm/yr, which adds to the velocity imposed with boundary conditions. Such active component develops even in the case of closed systems and determines a state of compressional stress within the lithosphere. The structure of the ascending flow in the melting regime below a spreading center suggests that it results from the superposition of two small-scale rolls with opposite polarity, associated with horizontal convection.
How to cite: Schettino, A.: The role of small-scale horizontal convection in lithosphere-asthenosphere interaction, EGU General Assembly 2025, Vienna, Austria, 27 Apr–2 May 2025, EGU25-12536, https://doi.org/10.5194/egusphere-egu25-12536, 2025.
The Thermal Lattice Boltzmann Method (TLBM) offers an alternative to classical PDE based solvers for planetary dynamics and is based on solving the Boltzmann Equation on a discrete lattice which maps perfectly onto parallel computers. We present strong scaling performance runs on the Shaheen III HPC cluster at KAUST using up to 300K cores for a 3D whole mantle model at a 3km grid resolution. Based on the throughput performance achieved, the TLBM can model mantle simulation in 3D for one convection cycle's worth of physical time in less than a day of CPU at a 3km resolution. We present 2D performance results which indicate that the TLBM can model ultra-high Rayleigh numbers in 2D well into the ultimate turbulent regime up to Ra = 1018 on 300K cores in a 1x1 aspect ratio model. We also present an update of the status of TLBM developments and example runs and run times of: (1) the transition to plate tectonics in the Archean using a temperature dependent viscosity, a yield stress formulation of the rheology and partial melting, (2) high Rayleigh number simulations in the turbulent regime up to Ra = 1015 in 2D, (3) 2D whole mantle modelling in a circular annulus and accuracy benchmarks against ASPECT, and (4) 3D simulations of whole mantle convection at a resolution of 30 km on just 96 cores, and (5) iron droplets from an impactor on a turbulent magma ocean settling to form the iron core using a combined TLBM and multiphase LBM. We believe that the TLBM and multiphase TLBM have the potential to lead to new insights in the dynamics and evolution of the Earth and exoplanets from the early lava world stage onwards including plate tectonics due to their high throughput performance and near linear scaling to 100s of K cores. These capabilities enable geodynamical modelling with never-before-seen resolutions in 2D and 3D, high Rayleigh numbers well into the ultimate turbulent regime, studies of turbulent magma oceans and core formation, and phase space studies of planetary dynamics.
How to cite: Mora, P., Morra, G., O'Neill, C., Honarbakhsh, L., Kuang, J., and Barges, A.: The Thermal Lattice Boltzmann Method: new developments, strong scaling to 300K cores, and potential to yield major advancements in geodynamics, EGU General Assembly 2025, Vienna, Austria, 27 Apr–2 May 2025, EGU25-20801, https://doi.org/10.5194/egusphere-egu25-20801, 2025.
Neutrino oscillation tomography is potentially a method for probing the properties of Earth's deep interior, complementing classical geophysical and geochemical methods. It relies on the detection of neutrinos, subatomic particles that interact weakly with matter and can traverse the Earth’s interior essentially unimpeded. Neutrinos exist in three types, called "flavors": electron, muon, and tau. As they propagate, they can change from one flavor to another, a phenomenon known as neutrino oscillation. Oscillation probabilities are influenced by the electron density profile along the neutrino’s path, determined by the matter density and the proton-to-nucleon ratio (Z/A) distribution. By measuring neutrino oscillations, it is thus possible to retrieve information about the composition and density variations in the Earth’s interior.
In this work, we present sensitivity kernels from neutrino oscillation tomography for a spherically symmetric Earth model. Our goal is to identify which depth ranges can be effectively studied using this technique. To understand the constraints that neutrino oscillation tomography can provide on Earth's structure, we first model the sensitivity of neutrino tomography to the planet's composition and density assuming an ideal neutrino detector. Then, to derive realistic sensitivities, we apply the detector’s response (i.e., resolution) of next-generation neutrino telescopes. We show that an ideal detector is most sensitive to the outer core, while realistic detectors with lower resolution but large detection volumes shift the sensitivity focus to shallower depths. Finally, we discuss how this method could provide complementary insights into the structure of large low velocity provinces (LLVPs) at the base of mantle and the water content in the mantle transition zone (MTZ).
How to cite: Goos, I. A., Coelho, J. A. B., Deniz Hernandez, Y. A., Durand, S., Fuji, N., Mittelstaedt, E. L., Pestes, R., and Van Elewyck, V.: Probing Earth’s interior with neutrinos: sensitivity kernels for a 1-dimensional Earth model, EGU General Assembly 2025, Vienna, Austria, 27 Apr–2 May 2025, EGU25-20446, https://doi.org/10.5194/egusphere-egu25-20446, 2025.
