Geothermal energy is a crucial component of the global transition to sustainable and green energy systems due to its renewable and long-term availability. In order to study potential resources, we need to describe the subsurface by solving inverse problems. The complexity and uncertainty of these problems require the use of probabilistic inversion approaches that repeatedly solve partial differential equations over a grid of parameters describing the subsurface domain. Frequently, the high dimensionality of the parameter space to be inferred implies prohibitive computational times and reduces the sensitivity of each parameter as the grid is refined. In this work, we implement and discuss adaptive parametrization strategies in Bayesian inversions. We model the thermal conductivity structure of 2D sections of the Earth's upper mantle and perform Markov chain Monte Carlo (MCMC) inversions to recover the thermal conductivity as a probability distribution based on the likelihood of the temperature measurements. To verify the solution, we first parametrize the physical properties of the subsurface domain equal to the high-dimensional finite element grid. In order to determine the optimal metaparameters on the run we rely on adaptive MCMC techniques that accelerate the convergence and reduce the risk of getting trapped in local minima. We then use a new parametrization based on the physical structure of the geological faults of the mantle that reduces the dimensionality of the problem. By relying on transdimensional sampling through reversible-jump MCMC, we consider the number of parameters as an unknown of the inversion. In these methods, the algorithm is allowed to increase the number of parameters to invert when the solutions found are not accurate enough and to decrease it when the accuracy of the solution is not significantly affected. Our results show that we recover the thermal conductivity structure both with and without adaptive parametrization, and the performance is improved when using transdimensionality. Moreover, the proposed transdimensional inversion decreases or increases the number of parameters locally, thereby providing an efficient and robust method for addressing the often challenging lack of information on subsurface heterogeneity.

Keywords: geothermal energy; Markov chain Monte Carlo; reversible jump MCMC; transdimensional inversion; adaptive parametrization; finite elements; Poisson equation.

How to cite: Dols, A., Amaya, M., Zlotnik, S., and Díez, P.: Adaptive parameterization in Bayesian inversions using transdimensional methods, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-12436, https://doi.org/10.5194/egusphere-egu26-12436, 2026.

Magma migration is a complex natural process that controls volcanism, the formation of many types of ore deposits, the development of geothermal reservoirs and the thermal structure, and long-term evolution of the lithosphere [1-3]. Because the dynamics of magma migration are difficult to observe directly, numerical simulations provide a powerful tool to investigate magmatic systems, the coupled physiochemical processes involved, and the range of spatial and temporal scales over which these processes operate.

In this study, we present a new multi-phase numerical framework to study magma migration within the Earth, with a particular emphasis on the mechanical interactions between melt and solid. The framework is based on multiphase flow in porous media and it incorporates realistic rheological descriptions of lithospheric rocks, including visco-elasto-viscoplastic behavior, damage, strain weakening and the generation of porosity due to plastic deformation. Interaction between the fluid (magma) and solid (host rock) phases are described via a set of equations derived from a formal phase-averaging framework. An arbitrary Eulerian-Lagrangian solver is used to discretize the equations and solve the fully-coupled system. The validity of the model, and its potential to study multi-scale magmatic systems, are demonstrated using well-known benchmark tests and targeted numerical experiments.

Keywords: Dynamics of lithosphere and mantle, Mechanics, Numerical modeling, Physics of magma, Plasticity

REFERENCES

• [1] Keller, D. A. May, and B. J. Kaus, “Numerical modelling of magma dynamics coupled to tectonic deformation of lithosphere and crust,” Geophys. J. Int., Vol. 195, pp. 1406-1442, (2013).
• [2] Li, A. E. Pusok, T. Davis, D. A. May, and R. F. Katz, “Continuum approximation of dyking with a theory for poro-viscoelastic-viscoplastic deformation,” Geophys. J. Int., Vol. 234, pp. 2007-2031, (2023).
• [3] Oliveira, J. C. Afonso, S. Zlotnik, and P. Diez, “Numerical modelling of multiphase multicomponent reactive transport in the Earth’s interior,” Geophys. J. Int., Vol. 212, pp. 345-388, (2018).

Acknowledgment

EarthSafe Doctoral Network has received funding from the European Union’s Horizon Europe research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 101120556.

How to cite: Hosseinian, N., Afonso, J. C., García-González, A., and Zlotnik, S.: Numerical modeling of magma migration in lithospheric rocks, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-5509, https://doi.org/10.5194/egusphere-egu26-5509, 2026.

A key challenge in the 21st century is the successful implementation of the energy transition, which hinges on the development of sustainable (energy) resources. In this context, hydrogen gas (H₂) generated by natural processes is a promising source of clean energy. However, we urgently need to develop the concepts and exploration strategies for this promise of natural H₂ energy to become a reality.

The most likely mechanism of large-scale natural H₂ generation in nature is the serpentinization of ultramafic mantle rocks during their chemical reaction with water. In order to predict the bulk serpentinization and natural H₂ generation that may lead to the development of exploitable H₂ deposits, we consider the following “recipe” for efficient serpentinization, which involves three main ingredients: (1) (fresh) mantle rocks that need to be at (2) optimal temperatures between ca. 200-350˚C (the serpentinization window), and (3) in contact with ample water for the reaction to take place. The serpentinization window can be expected at 8-12 kilometers below the Earth’s surface. However, mantle rocks are normally found at much greater depth; thus these rocks must be brought closer to the surface (exhumed) through geodynamic processes. Moreover, water needs to reach such depths along large faults or other structures that cut into the exhumed mantle. The challenge we are faced with is to understand where (and when) these ingredients may come together in nature, and how much natural H₂ may be generated.

Numerical geodynamic modelling is an ideal means to tackle this issue since it allows us not only to test how mantle rocks can be exhumed, but also to trace the temperature conditions and potential water availabilitiy (Zwaan et al. 2025). By combining this information, we assess favorable settings and timing of bulk natural H₂ generation in different geodynamic systems. Subsequently, we consider where the natural H₂ could be exploited. The serpentinizing mantle source rocks at 8-12 km depth cannot be directly targeted. Ideally, the natural H₂ would instead migrate and accumulate in sedimentary reservoir rocks at depths of only a couple of kilometers that are connected with the mantle source rocks via migration pathways (e.g., faults). Importantly, all key elements need to be in place for the system to work.

Our first-order modelling work and the development of natural H₂ system concepts greatly helps to direct natural H₂ resource exploration efforts, for example in the Alps and Pyrenees. Moreover, substantial opportunity lies in refining both the geodynamic modelling and natural H₂ system analysis, in field- and laboratory testing of our H₂ system concepts, and in extending such a “mineral system” modelling approach to other types of natural resources that are crucial to the energy transition. 

Reference:

Zwaan, F., Brune, S., Glerum, A.C., Vasey, D.A., Naliboff, J.B., Manatschal, G., & Gaucher, E.C. 2025: Rift-inversion orogens are potential hot spots for natural H₂ generation, Science Advances, 11, eadr3418. https://doi.org/10.1126/sciadv.adr3418

How to cite: Zwaan, F., Glerum, A. C., Brune, S., Vasey, D. A., Naliboff, J. B., Manatschal, G., and Gaucher, E. C.: Numerical geodynamic modelling for natural H2 resource exploration, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-3352, https://doi.org/10.5194/egusphere-egu26-3352, 2026.

The prevalence of viscous-dominated regimes within the Earth’s interior gives rise to Stokes-like flow systems in numerous geodynamical applications. A prominent example is sublithospheric mantle convection, which constitutes the primary driving mechanism behind the evolution of dynamic topography. In this context, numerical simulations provide more physically consistent estimates of the Lithosphere–Asthenosphere Boundary (LAB) depth than those derived from first-order isostatic approximations [1].

However, the associated computational overburden is exceptionally high, particularly when accounting for material nonlinearities. The challenge is further complicated when attempting to incorporate them within a Markov Chain Monte Carlo (MCMC) framework that requires an exceptionally large number of evaluations [2], limiting their applicability to large-scale studies and underscores the need for novel and computationally efficient Reduced-Order Modeling (ROM) methodologies [3].

Results from linear Model Order Reduction (MOR) techniques indicate that the complexity of the problem surpasses the capabilities of projection-based ROMs designed to produce globally accurate solutions. This work introduces a localized, goal-oriented criterion to enhance linear reducibility and employs Neural Network (NN) surrogates to replace high-fidelity solver evaluations. These methodological advances jointly underpin the development of a hybrid offline–online reduction framework that efficiently reduces computational complexity while preserving the required levels of accuracy, enabling seamless model updates during parameter-space exploration.

REFERENCES

[1] Afonso, J. C., Rawlinson, N., Yang, Y., Schutt, D. L., Jones, A. G., Fullea, J., & Griffin, W. L. (2016). 3-D multiobservable probabilistic inversion for the compositional and thermal structure of the lithosphere and upper mantle: III. Thermochemical tomography in the Western-Central U.S. Journal of Geophysical Research: Solid Earth, 121(10), 7337–7370. https://doi.org/10. 1002/2016jb013049

[2] Ortega-Gelabert, O., Zlotnik, S., Afonso, J. C., & Diez, P. (2020). Fast Stokes Flow Simulations for Geophysical-Geodynamic Inverse Problems and Sensitivity Analyses Based on Reduced Order Modeling. Journal of Geophysical Research: Solid Earth, 125(3). https://doi.org/10.1029/ 2019jb018314

[3] Hesthaven, J.S., Rozza, G., Stamm, B. (2015). Certified Reduced Basis Methods for Parametrized Partial Differential Equations. SpringerBriefs in Mathematics. Springer International Publishing AG, Cham. https://doi.org/10.1007/978-3-319-22470-1

How to cite: Ramadan, M., Pichi, F., and Rozza, G.: Toward Efficient Stokes Flow Simulations in Multi-Observable Thermo-Chemical Tomography Using Model Order Reduction, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-6024, https://doi.org/10.5194/egusphere-egu26-6024, 2026.

Multi-Observable Thermochemical Tomography (MTT) is a simulation-driven, joint probabilistic inversion framework designed to estimate the thermochemical state of the Earth’s lithosphere by integrating geophysical datasets with complementary sensitivities. By jointly inverting observables such as gravity and geoid anomalies, surface heat flow, seismic dispersion, body-wave data, and magnetotelluric responses, MTT directly estimates primary thermodynamic variables, including temperature, pressure, and bulk composition, from which all secondary physical properties are derived through internally consistent thermodynamic models. This bottom-up approach provides physically-consistent constraints on lithospheric structure across regional to prospect scales.

Within this framework, MTT offers a powerful basis for characterizing lithospheric architecture and compositional domains that are commonly examined in mineral systems studies. In particular, MTT can help delineate major crustal- and lithospheric-scale structures, identify metasomatized/altered domains, and map thermochemical contrasts that serve as lithospheric-scale proxies commonly associated with specific classes of magmatic and hydrothermal mineral systems.

Despite recent advances incorporating ray-based seismic tomography solvers (Fomin, I., Afonso, J. C., Gorbatov, A., Salajegheh, F., Dave, R., Darbyshire, F. A., et al. (2026). Multi-observable thermochemical tomography: New advances and applications to the superior and North Australian cratons. Journal of Geophysical Research: Solid Earth, 131, e2025JB031939. https://doi.org/10.1029/2025JB031939 ), the integration of full-waveform seismic solvers within the MTT framework has not yet been achieved. Full-waveform inversion (FWI) offers enhanced sensitivity to both seismic velocity and density and the potential for improved spatial resolution relative to traditional tomography approaches. However, the computational cost of FWI remains prohibitive, particularly in probabilistic or ensemble-based inversion settings required for uncertainty quantification.

This contribution presents a computational strategy aimed at reducing the cost of full wavefield simulations to enable probabilistic seismic FWI within the MTT framework. We extend reduced-order modeling (ROM) techniques to the spectral element method (SEM), which is widely used for accurate time-domain seismic wave propagation in complex geological settings. Specifically, we consider projection (Galerkin)–based ROMs in which the SEM wavefield is approximated in a low-dimensional reduced basis constructed from representative high-fidelity solutions. While ROM approaches are well established for simpler formulations, their application to SEM-based elastic wave simulations remains challenging due to the method’s high dimensionality and complex operator structure. Beyond MTT, such reductions are also relevant to SEM-based workflows that require large numbers of forward simulations, including ground motion studies and FWI with many sources at regional-to-global scales.

We develop and test a reduced-order SEM formulation using synthetic benchmark models relevant to lithospheric-scale imaging. Results demonstrate computational speed-ups of up to two orders of magnitude relative to full SEM simulations, while retaining sufficient accuracy in simulated wavefields for inversion purposes. These results represent a first proof of concept toward incorporating probabilistic FWI into multi-observable thermochemical tomography and reducing a key computational barrier to uncertainty-aware, physics-based lithospheric imaging.

How to cite: Jamasb, A., Afonso, J.-C., Garcia Gonzalez, A., Rozza, G., Pichi, F., Zlotnik, S., Meijde, M. V. D., and Fadel, I.: Enabling Probabilistic Full Waveform Inversion in Multi-Observable Thermochemical Tomography through Reduced-Order Spectral Element Modeling, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-13298, https://doi.org/10.5194/egusphere-egu26-13298, 2026.

Mantle convection plays a fundamental role in governing the thermal and dynamical
evolution of terrestrial planets, yet its numerical simulation remains computationally ex-
pensive due to strong nonlinearities, high Rayleigh numbers, and the presence of thin
thermal boundary layers. In this work, we present a non-intrusive reduced-order modeling
(ROM) framework for two-dimensional mantle convection based on Proper Orthogonal
Decomposition combined with Radial Basis Function interpolation (POD–RBF).
High-fidelity full-order model (FOM) simulations are first performed using a finite-
volume discretization of the incompressible Boussinesq equations under the infinite-Prandtl-
number approximation. The FOM is carefully validated across a wide range of Rayleigh
numbers. Particular attention is devoted to high-Rayleigh-number regimes, where mesh
refinement studies are conducted to improve accuracy and ensure reliable reference solu-
tions.
The ROM is constructed from snapshot data of velocity and temperature fields. POD
analysis reveals a rapid decay of singular values, indicating a low-dimensional structure
of the solution manifold. The parametric dependence of the reduced coefficients is recon-
structed using RBF interpolation, yielding a fully data-driven and non-intrusive ROM.
To rigorously assess predictive capability, the ROM is validated using test points ex-
cluded from the training dataset. Leave-One-Out cross-validation demonstrates that the
ROM accurately predicts unseen solutions across the parameter space, with low relative
L2 errors for both velocity and temperature fields. Field-level comparisons confirm that
the dominant flow structures and thermal patterns are faithfully reproduced.
The framework is further extended to transient simulations, where both time and
Rayleigh number are treated as parameters. This two-dimensional parametric unsteady
ROM successfully captures time-dependent dynamics while providing significant compu-
tational speed-up. The proposed approach offers a robust and efficient tool for parametric
mantle convection modeling and provides a solid basis for future extensions toward three-
dimensional configurations and uncertainty quantification.

How to cite: Haider, Q., Tonicello, N., Girfoglio, M., and Rozza, G.: Non-Intrusive POD–RBF Reduced OrderModeling for Parametric and Transient MantleConvection, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-13256, https://doi.org/10.5194/egusphere-egu26-13256, 2026.

Geophysical inverse problems are inherently ill-posed due to sparse, noisy, and indirect observations, making Uncertainty Quantification (UQ) a fundamental requirement for reliable subsurface characterization. Bayesian inversion provides a comprehensive probabilistic framework for inferring subsurface parameters by coherently combining prior knowledge with observational data through the likelihood function. However, the practical deployment of Bayesian methods in large-scale geophysical settings is often hampered by the prohibitive computational cost of repeated forward model evaluations. In this context, uncertainty is often not solely driven by observational noise; a substantial and sometimes dominant contribution arises from model error, resulting from simplified physical descriptions, numerical discretization, and uncertain boundary conditions. When these sources of uncertainty are neglected or inadequately represented, Bayesian inversions may yield biased posterior estimates and unrealistically narrow uncertainty bounds. These limitations are particularly acute in deep Earth applications, where complex rheologies, poorly constrained geometries, and computationally intensive forward models coexist.

A key challenge is the accurate delineation of the Lithosphere–Asthenosphere Boundary (LAB), which plays a central role in controlling mantle dynamics, lithospheric deformation, and deep geothermal processes. Despite the necessity of relying on Bayesian approaches to estimate the LAB and its associated uncertainties, the high computational cost of repeated evaluations of the forward solver makes this unfeasible within realistic time frames [1]. To address these limitations, this work investigates Reduced-Order Modeling (ROM) techniques to enable efficient Bayesian inversion of LAB geometry in geodynamical Stokes flow models. ROMs construct low-dimensional surrogates of high-fidelity solvers, allowing rapid forward simulations while preserving the dominant physical behavior of mantle flow. By integrating ROMs with Bayesian inference, the proposed framework enables effective and reliable UQ for LAB characterization.
Keywords: Geophysical inverse problems; Bayesian inversion; Uncertainty Quantification; Reduced-Order Modeling; Lithosphere–Asthenosphere Boundary

Acknowledgement This research was conducted within the EarthSafe Doctoral Network and has received funding from the European Union’s Horizon Europe research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 101120556.

References [1] Olga Ortega-Gelabert, Sergio Zlotnik, Juan Carlos Afonso, and Pedro 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, 125(3):e2019JB018314, 2020.

How to cite: Shahzaib, M., Díez, P., Zlotnik, S., Muixí, A., and Amaya, M.: Coupling Bayesian Inversion and Reduced-Order Modeling: Application to Lithosphere–Asthenosphere Boundary Estimation, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-8341, https://doi.org/10.5194/egusphere-egu26-8341, 2026.

Three-dimensional (3D) Magnetotelluric (MT) probabilistic inversion remains rare in real-world applications because it requires solving the forward problem thousands to millions of times, often making the computational cost prohibitive. Since the total duration of an inversion is directly controlled by the performance of the forward solver, the high computational overhead of 3D MT modeling remains a significant challenge, particularly for large-scale problems requiring high mesh resolutions. To address the poor scaling of existing strategies, we introduce DD–POD, a hybrid framework that integrates Domain Decomposition (DD) with Proper Orthogonal Decomposition (POD). The DD formulation partitions the global problem into subdomains, bypassing the memory limitations of traditional direct solvers and enabling simulations with substantially finer discretizations. Implementing this distributed architecture alone yields simulations that are at least 50% faster than global full-order approaches. Building on this foundation, the integration of POD eliminates the need for repeated large-scale linear system solves within the iterative DD process, delivering total forward-solver speed-ups exceeding 90%. Benchmark experiments and a real-world case study demonstrate that DD–POD consistently outperforms standard global POD strategies in computational efficiency with an acceptable trade-off in numerical accuracy.

(This work was supported by the Marie Sklodowska-Curie Actions (Doctoral Network with Grant agreement No. 101120556))

How to cite: Tao, L., Zlotnik, S., Muixí, A., Zyserman, F. I., Afonso, J. C., and Diez, P.: Reducing Computational Costs in 3D Magnetotelluric Simulations via Domain Decomposition and Reduced-Order Modeling, EGU General Assembly 2026, Vienna, Austria, 3–8 May 2026, EGU26-10570, https://doi.org/10.5194/egusphere-egu26-10570, 2026.