GD8.1

During the last decade, the development of numerical geodynamic tools helped the geosciences community to unravel complex thermo-mechanical processes at play during plate tectonics. Yet, the high computational cost of thermodynamic calculations, which simulates phase change, hampers our ability to integrate complex chemistry in such problems. This is particularly important for simulating magmatic processes, where the chemistry of differentiating melts can vary significantly from the mantle to the upper crust. The typical approach, currently used, is to precompute one or many phase diagrams and use them as look-up tables. For many geodynamic processes this is adequate but when the melt chemistry varies drastically it would be better to be able to do thermodynamic calculations on the fly, along with the geodynamic models.

For that, the thermodynamic computational approach must be sufficiently fast, should work fully automatically and be tuned for melting models of magmatic systems, for example by utilizing the recently developed thermodynamic melting model of Holland et al. (2018). Existing approaches do not fulfill all criteria, which is why we have developed a new computational library for this purpose. Our code is written in C, runs on massively parallel machines (MPI) and uses an adaptive mesh refinement strategy to compute phase diagrams. At the moment we have focused on the 'igneous set' of the Holland & Powell dataset (as defined in the thermocalc software) to calculate stable phase equilibria in the system K2O–Na2O–CaO–FeO–MgO–Al2O3–SiO2–H2O–TiO2–Fe2O3–Cr2O3 (KNCFMASHTOCr). The code uses pressure, temperature and bulk-rock composition as input and returns relevant petrological and geodynamic information such as (but not restricted to) stable assemblage, phase fractions and phase densities. Different than many of the existing approaches, our method can efficiently utilize initial guesses which naturally occur in geodynamic simulations where the changes in chemistry between timesteps are usually minor.

The methodology performs a Gibbs free energy minimization and involves two main steps. First, we use a combination of levelling methods (iterative change of base) to reduce the number of potential (pure and solution) phases and to bring the G-hyperplane close to solution. Second, we use a partitioning of Gibbs energy approach coupled with local minimization to satisfy the Gibbs-Duhem rule and to retrieve the final set of stable solution phases. To illustrate the efficiency of the library up to supra-solidus conditions we present a set of dry phase diagrams and compare results of our computations with thermocalc calculations.

Ongoing development includes the treatment of solvus to extend its applicability to complex wet systems involving solution phase such as amphibole.

How to cite: Riel, N., Kaus, B., Green, E., Berlie, N., and Rummel, L.: A new and efficient computational thermodynamics approach for magmatic systems, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-3910, https://doi.org/10.5194/egusphere-egu21-3910, 2021.

Thermodynamics plays an increasingly important role in computational geodynamics, as the community moves towards including chemical reactions in mechanical models of deformation. An example is mineral reactions that induce volume changes which affect the local state of stress that may trigger nonlinear feedbacks. Another example is numerical simulations of magmatic systems where the chemistry of melts changes dramatically as the melt differentiates on its way up through the lithosphere. In order to compare the results of numerical models with geochemical and petrological field data, it is crucial to predict the rock types and major element chemistry from numerical simulations of magmatic systems. Precomputing phase diagrams and including them as a lookup table, which is the current standard approach, works when the chemistry is more or less constant, but is nontrivial for magmatic systems where the chemistry changes drastically and involves a 13-dimensional system (11 oxides, plus pressure and temperature). Parameterizing the main reactions is a possibility, but this renders the comparison of simulation results with available data difficult.

In the context of magmatic systems, significant progress has been made in recent years and we now have thermodynamic melting models that are consistent with experiments and can simulate the full compositional range from mafic to felsic melt compositions. Yet, in order to be useful for geodynamic applications, we also need sufficiently fast Gibbs energy minimization software tools that can automatically determine the most stable assemblage for a given pressure, temperature and chemistry. This is a nonlinear constrained optimization problem, for which we have developed a new, parallel, solution approach. One novelty of our approach is that it can efficiently make use of good initial guesses, for example obtained from the previous timestep of a geodynamic simulation or from nearby points in pressure and temperature space. Yet, as with any gradient-based method, a risk remains to be trapped in a local minimum in the solution space that is not the overall most stable assemblage. It is thus important to explore a sufficiently broad range of starting parameters to ensure convergence toward the global minimum. Whereas this problem is well-known among users of existing thermodynamic software (such as THERMOCALC, where the starting values have to be adjusted manually), it is much less clear how nonlinear the parameter space actually is for real applications. Do we have thousands of local minima, or are there only a few (and if yes, can we precompute some of these)?

Here, we will discuss several examples of melting models and map out the nonlinearity of the parameter space for these cases to get better insights in how to further speed up such calculations. We also discuss how shallow and deep neural networks can be trained and implemented as part of the workflow.

How to cite: Kaus, B. and Riel, N.: Neural networks, local minima and computational thermodynamics, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-5173, https://doi.org/10.5194/egusphere-egu21-5173, 2021.

The dynamics of magmatic systems remain poorly understood, due to the lack of resolving power of geophysical methods to study active systems and the difficulty of interpreting exposed crystallized magma bodies. Numerical models are therefore helpful to connect the dots between classical geological studies, using rheological information and geometries derived from field or geophysical investigations to shed new lights on the mechanisms involved in such systems.

Taking advantage of the big CPU clusters currently available and the development of the DMStag framework as part of the PETSc infrastructure, the ERC-funded MAGMA project aims to build tools to analyse magmatic processes in the lithosphere. We developed a finite-difference staggered grid code solving the Stokes equations for visco-elasto-plastic rheologies and using analytical jacobians for linear and non-linear solvers, combined with regularized plasticity. The code is combined with both a marker and cell and semi-lagrangian advection schemes, is fully parallel and includes automated testing.

Here, we provide application examples ranging from simple benchmark validations against analytical solutions to more complex settings taking advantage of the broad rheologies and local heterogeneities permitted by high resolution settings and the finite difference method. Ongoing technical developments include adding two-phase flow and coupling to it with thermodynamic calculations to track the evolving chemistry of magmatic systems.

How to cite: Berlie, N., Kaus, B., Popov, A., Arts, M., Riel, N., and Kiss, D.: Towards modelling magmatic systems using a staggered grid finite difference method, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-4271, https://doi.org/10.5194/egusphere-egu21-4271, 2021.

It is broadly accepted that magmatism plays a key dynamic role in continental and oceanic rifting. However, these dynamics remain poorly studied, largely due to the difficulty of consistently modelling liquid/solid interaction across the lithosphere. The RIFT-O-MAT project seeks to quantify the role of magma in rifting by using models that build upon the two-phase flow theory of magma/rock interaction. A key challenge is to extend the theory to account for the non-linear rheological behaviour of the host rocks, and investigate processes such as diking, faulting and their interaction. Here we present our progress in consistent numerical modelling of poro-viscoelastic–plastic modelling of deformation with a free surface.

Failure of rocks (plasticity) is an essential ingredient in geodynamics models because Earth materials cannot sustain unbounded stresses. However, plasticity represents a non-trivial problem even for single-phase flow formulations (Spiegelman et al. 2016). The elastic deformation of rocks can also affect the propagation of internal failure. Furthermore, deformation and plastic failure drives topographic change, which imposes a significant static stress field. Robustly solving a discretised model that includes this physics presents severe challenges, and many questions remain as to effective solvers for these strongly nonlinear systems. 

We present a new finite difference staggered grid framework for solving partial differential equations (FD-PDE) for single-/two-phase flow magma dynamics (Pusok et al., 2020). Staggered grid finite-difference methods are mimetic, conservative, inf-sup stable and with small stencil — thus they are well suited to address these problems. The FD-PDE framework uses PETSc (Balay et al., 2020) and aims to separate the user input from the discretization of governing equations. The core goals for the FD-PDE framework is to allow for extensible development and implement a framework for rigorous code validation. Here, we present simplified model problems using the FD-PDE framework for two-phase flow visco-elasto-plastic models designed to characterise the solution quality and assess both the discretisation and solver robustness. We also present results obtained using the phase-field method (Sun and Beckermann, 2007) for representing the free surface. Verification of the phase-field approach will be shown via simplified problems previously examined in the geodynamics community (Crameri et al, 2012).

Balay et al. (2020), PETSc Users Manual, ANL-95/11 - Revision 3.13.

Pusok et al. (2020) https://doi.org/10.5194/egusphere-egu2020-18690 

Spiegelman et al. (2016) https://doi.org/10.1002/2015GC006228

Sun and Beckermann (2007) https://doi.org/10.1016/j.jcp.2006.05.025

Crameri et al. (2012) https://doi.org/10.1111/j.1365-246X.2012.05388.x

How to cite: Li, Y., Pusok, A., May, D., and Katz, R.: Free surface models of partially molten rock with visco-elasto–plastic rheology: numerical solutions using a staggered-grid finite-difference discretisation, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-7631, https://doi.org/10.5194/egusphere-egu21-7631, 2021.

Earthquake fault zones are more complex, both geometrically and rheologically, than an idealised infinitely thin plane embedded in linear elastic material.  Field and laboratory measurements reveal complex fault zone structure involving tensile and shear fractures spanning a wide spectrum of length scales (e.g., Mitchell & Faulkner, 2009), dense seismic and geodetic recording of small and large earthquakes show hierarchical volumetric faulting patterns (e.g., Cheng et al., 2018, Ross et al., 2019) and 2D numerical models explicitly accounting for off-fault fractures demonstrate important feedback with rupture dynamics and ground motions (e.g., Thomas & Bhat 2018, Okubo et al., 2019).

Here (Gabriel et al., 2021) we adopt a diffuse crack representation to incorporate finite strain nonlinear material behaviour, natural complexities and multi-physics coupling within and outside of fault zones into dynamic earthquake rupture modeling. We use a first-order hyperbolic and thermodynamically compatible mathematical model, namely the GPR model (Godunov & Romenski, 1972; Romenski, 1988),  to describe a continuum in a gravitational field which provides a unified description of nonlinear elasto-plasticity, material damage and of viscous Newtonian flows with phase transition between solid and liquid phases.

The model shares common features with phase-field approaches but substantially extends them. Pre-damaged faults as well as dynamically induced secondary cracks are therein described via a scalar function indicating the local level of material damage (Tavelli et al., 2020); arbitrarily complex geometries are represented via a diffuse interface approach based on a solid volume fraction function (Tavelli et al., 2019). Neither of the two scalar fields needs to be mesh-aligned, allowing thus faults and cracks with complex topology and the use of adaptive Cartesian meshes (AMR). High-order accuracy and adaptive Cartesian meshes are enabled in 2D and 3D by using the extreme scale hyperbolic PDE solver ExaHyPE (Reinarz et al., 2019).

We show a wide range of numerical applications that are relevant for dynamic earthquake rupture in fault zones, including the co-seismic generation of secondary off-fault shear cracks, tensile rock fracture in the Brazilian disc test, as well as a natural convection problem in molten rock-like material. We compare diffuse interface fault models of kinematic cracks, spontaneous dynamic rupture and dynamically generated off-fault shear cracks to sharp interface reference models. To this end, we calibrate the GPR model to resemble empirical tensile and shear crack formation and friction laws. We find that the continuum model can resemble and extend classical solutions, while introducing dynamic differences (i) on the scale of pre-damaged/low-rigidity fault zone, such as out-of- plane rupture rotation; and (ii) on the scale of the intact host rock, such as conjugate shear cracking in tensile lobes. 

Our approach is part of the TEAR ERC project (www.tear-erc.eu) and will potentially allow to fully model volumetric fault zone shearing during earthquake rupture, which includes spontaneous partition of fault slip into intensely localized shear deformation within weaker (possibly cohesionless/ultracataclastic) fault-core gouge and more distributed damage within fault rocks and foliated gouges.

How to cite: Gabriel, A.-A., Li, D., Chiocchetti, S., Tavelli, M., Peshkov, I., Romenski, E., and Dumbser, M.: A unified first order hyperbolic model for nonlinear dynamic rupture processes in diffuse fracture zones, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-15237, https://doi.org/10.5194/egusphere-egu21-15237, 2021.

Rifts and rifted margins result from interaction of several physical processes, which produce a range of crustal structures, subsidence histories, and sedimentary architectures. Study of these processes in academia and industry includes kinematic modelling (i.e. cross-section restoration, backstripping) combined with simple thermomechanical models and dynamic modelling. In kinematic models, the thinning of the lower crust and mantle is kinematically imposed in the form of pure shear, which contradicts natural non-linear viscous behavior. Although, kinematic modelling can provide a crustal thinning profile, heatflow estimates, subsidence rates etc., imposed extension of the lower crust and mantle might strongly impact the result. On the other hand, a dynamic approach allows to model the whole range of possible physical processes, but it cannot be used to model particular extension histories.
Here, we show a new modelling technique, namely KineDyn, to combine the advantages of the above-mentioned approaches into a single modelling framework. Our method employs full non-linear visco-elasto-plastic rheology, surface process of erosion and sediment transport, decompression melting of the mantle, and serpentinization of mantle rocks. Faults are introduced as weak planes in the upper crust, in order to simulate faulting during the model run. In our approach, faults are initially controlled by prescribed initial locations, offsets and timings, while the rest of the model is resolved in a fully dynamic mode. Since fault planes are much weaker than the surrounding upper crust, extension of the model naturally leads to slip on the faults. We demonstrate that faults modelled this way reproduce a natural behavior, including rotation due to flexure and unloading of the fault plane.
In order to reconstruct the evolution of an existing rift or rifted margin we model extension of the lithosphere with controlled faulting. To do this we use the interpreted spatio-temporal evolution of the faulting from a seismic profile to guide the evolution of the dynamic model. After a trial-and-error process, where we correct the faults’ locations, the thicknesses of layers, surface process’s parameters, initial thermal gradient etc., we obtain the model that best fits the observations. Thus, KineDyn gives, in effect, the same results as existing section restoration techniques (i.e. the potential history of faulting) and forward modeling techniques (i.e. the likely history of sedimentation, thinning, heat flow and subsidence), while simultaneously taking into account non-linear interactions between processes occurring during rifting.
In this work we show the methodology, examples, tests and benchmarks of the technique. Finally, we present applications of KineDyn for the following rifts and rifted margins: Malawi Rift, East African Rift System, hyper-extended West Iberia Margin, and ultra-wide Santos-Benguela Rifted Margin.

How to cite: Muldashev, I., Pérez-Gussinyé, M., Neto Cavalcanti de Araújo, M., and Liu, Z.: KineDyn: Thermo-Mechanical Method for Validation of Seismic Interpretations with Dynamic Forward Models, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-8541, https://doi.org/10.5194/egusphere-egu21-8541, 2021.

In the marker-in-cell method combined with staggered finite differences, Lagrangian markers carrying information on material properties are advected with the velocity field interpolated from the staggered Eulerian velocity grid. With existing schemes, velocity interpolation from the grid points to markers violates (to some extent) mass conservation requirement that causes excess convergence/divergence of markers and opening marker gaps after significant amount of advection. This effect is especially well visible in case of diagonal simple shear deformation along planes that are oriented at 45 degrees to the grid and marker circulation through grid corners.

Here, we present a new second order velocity interpolation scheme that guarantees exact interpolation of normal strain rate components from pressure nodes (i.e. from the locations where these components are defined by solving of the mass conservation equation). This new interpolation scheme is thus applicable to both compressible and incompressible flow and is trivially expendable to 3D and to non-regular staggered grids.

Performed tests show that, compared to other velocity interpolation approaches, the new scheme has superior performance in preserving continuity of the marker field during the long-term advection including the diagonal simple shear deformation and marker circulation through grid corners. We showcase a performance-oriented implementation of the new scheme using Julia language's shared memory parallelisation features. The Julia implementation of the new advection schemes further augments the ParallelStencil.jl  related application collection with advection routines.

How to cite: Gerya, T., Duretz, T., and Ludovic, R.: New continuity-based velocity interpolation scheme for staggered grids, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-15308, https://doi.org/10.5194/egusphere-egu21-15308, 2021.

In order to treat a free surface in models of lithosphere and mantle dynamics that use a fixed Eulerian grid it is typical to use "sticky air", a layer of low-viscosity, low-density material above the solid surface (e.g. Crameri et al., 2012). This can, however, cause numerical problems, including poor solver convergence due to the huge viscosity jump and small time-steps due to high velocities in the air. Additionally, it is not completely realistic because the assumed viscosity of the air layer is typically similar to that of rock in the asthenosphere so the surface is not stress free.  

In order to overcome these problems, Duretz et al. (2016) introduced and tested a method for treating the free surface that instead detects and applies special conditions at the free surface. This avoids the huge viscosity jump and having to solve for velocities in the air. They applied it to a two-dimensional staggered grid finite difference / finite volume scheme, a discretization that is in common use for modelling mantle and lithosphere dynamics. Here I document the application of this approach to a three-dimensional spherical staggered grid solver in the mantle simulation code StagYY. Some adjustments had to be made to the two-dimensional scheme documented in Duretz et al. (2016) in order to avoid problems due to undefined velocities for certain boundary topographies. The approach was applied not only to the Stokes solver but also to the temperature solver, including the implementation of a mixed radiative/conductive boundary condition applicable to surface magma oceans/lakes.

References

Crameri, F., H. Schmeling, G. J. Golabek, T. Duretz, R. Orendt, S. J. H. Buiter, D. A. May, B. J. P. Kaus, T. V. Gerya, and P. J. Tackley (2012), A comparison of numerical surface topography calculations in geodynamic modelling: an evaluation of the ‘sticky air’ method, Geophysical Journal International,189(1), 38-54, doi:10.1111/j.1365-246X.2012.05388.x.

Duretz, T., D. A. May, and P. Yamato (2016), A free surface capturing discretization for the staggered grid finite difference scheme, Geophysical Journal International, 204(3), 1518-1530, doi:10.1093/gji/ggv526.

How to cite: Tackley, P.: A Three-Dimensional Version of the Free Surface Capturing Discretization for Staggered Grid Finite Difference Schemes: Implementation into StagYY, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-16228, https://doi.org/10.5194/egusphere-egu21-16228, 2021.

How to cite: Zlotnik, S., Ortega, O., Díez, P., and Afonso, J. C.: A combined Reduced Order-Bayesian scheme to drastically accelerate stochastic inversions, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-8824, https://doi.org/10.5194/egusphere-egu21-8824, 2021.

The thermal evolution of terrestrial planets depends strongly on several parameters and initial conditions that are poorly constrained. Often, direct or indirect observables from planetary missions such as elastic lithospheric thickness, crustal thickness and duration of volcanism are inverted to infer the unknown parameter values and initial conditions. The non-uniqueness and non-linearity of this inversion necessitates a probabilistic inversion framework. However, due to the expensive nature of forward dynamic simulations of thermal convection , Markov Chain Monte Carlo methods are rarely used. To address this shortcoming, some studies have recently shown the effectiveness of Mixture Density Networks (MDN) (Bishop 1995) in being able to approximate the posterior probability using only the dataset of simulations run prior to the inversion (Meier et al. 2007, de Wit et al. 2013, Käufl et al. 2016, Atkins et al. 2016).

Using MDNs, we systematically isolate the degree to which a parameter can be constrained using different “present-day” synthetic observables from 6130 simulations for a Mars-like planet. The dataset – generated using the mantle convection code GAIA (Hüttig et al. 2013)- is the same as that used by Agarwal et al. (2020) for a surrogate modelling study.

The loss function used to optimize the MDN (log-likelihood) provides a single robust quantity that can be used to measure how well a parameter can be constrained. We test different numbers and combinations of observables (heat flux at the surface and core-mantle boundary, radial contraction, melt produced, elastic lithospheric thickness, and duration of volcanism) to constrain the following parameters: reference viscosity, activation energy and activation volume of the diffusion creep rheology, an enrichment factor for radiogenic elements in the crust, and initial mantle temperature. If all observables are available, reference viscosity can be constrained to within 32% of its entire range (10¹⁹−10²² Pa s), crustal enrichment factor (1−50) to within 15%, activation energy (10⁵−5×10⁵ J mol-1 ) to within 80%, and initial mantle temperature (1600−1800K) to within 39%. The additional availability of the full present-day temperature profile or parts of it as an observable tightens the constraints further. The activation volume (4×10^-6 −10×10^-6  m³ mol^-1) cannot be constrained and requires research into new observables in space and time, as well as fields other than just temperature. Testing different levels of uncertainty (simulated using Gaussian noise) in the observables, we found that constraints on different parameters loosen at different rates, with initial temperature being the most sensitive. Finally, we present how the marginal MDN proposed by Bishop (1995) can be modified to model the joint probability for all parameters, so that  the inter-parameter correlations and the associated degeneracy can be capture, thereby providing a more comprehensive picture of all the evolution scenarios that fit given observational constraints.

How to cite: Agarwal, S., Tosi, N., Kessel, P., Padovan, S., Breuer, D., and Montavon, G.: Towards constraining Mars' thermal evolution using machine learning, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-4044, https://doi.org/10.5194/egusphere-egu21-4044, 2021.

Advanced numerical methods and increasing computational power have enabled us to incorporate numerical forward models into geodynamic inverse frameworks. We now have several strategies to constrain the rheological properties of the crust and lithosphere. Yet, the initial geometry of geological formations (e.g., salt bodies, magma bodies, subducting slabs) and associated uncertainties are, in most cases, excluded from the inverse problem and assumed to be part of the a priori knowledge. Usually, geometrical properties remain constant, or we employ simplified bodies like planes, spheres or ellipsoids for their parameterization.

Here, we present a simple method to parameterize complex three-dimensional bodies and incorporate them into geodynamic inverse problems. Our approach enables us to automatically create an entire ensemble of initial geometries within the uncertainty limits of geophysical imaging data. This not only allows us to account for geometric uncertainties, but also enables us to investigate the sensitivity of geophysical data to the geometrical properties of the geological structures.

We present 3 areas of application for our method, covering salt diapirs, magmatic systems and subduction zones, using both synthetic and real data.

How to cite: Spang, A., Baumann, T. S., and Kaus, B. J. P.: Geodynamic inversion with uncertain initial geometries, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-5680, https://doi.org/10.5194/egusphere-egu21-5680, 2021.

Commonly, the physical properties of the Earth (e.g., velocity, density) are parameterized as continuous fields. The most popular representation are grids and basis functions like spherical harmonics or splines. In an inversion context it is quite common that not all the parameters are fully constrained by the available inputdata. This relates to the common issues of insufficient resolution, incomplete coverage, and trade-offs due tonon-uniqueness. By applying some form of regularization to the inverse problem, a well-behaved and unique solution can be obtained, but this solution depends on the details of the chosen regularization.

Transdimensional approaches address the regularization problem by using a model representation with a variable number of parameters. The number of parameters is adjusted according to the requirements of the input data using the reversible jump Monte Carlo Markov Chain (rj-MCMC) algorithm. The output is an ensemble of variable resolution models that provides insight into the required model complexity and trade-offbetween parameters.

Here, I present synthetic tests from a joint inversion of satellite gravity gradients and normal modes for the Earth's velocity and density structure. The mantle's seismic velocity and density inside a 2-D spherical annulus are described by a variable number of discrete anomalous volumes, each with a variable size, shape, location and strength of velocity and density anomaly. The discrete anomalies are adjusted using the transdimensional approach in order to fit the gravity and normal mode data. This synthetic example shows promising results, because the synthetic model can recovered reasonably well.

How to cite: Szwillus, W.: Possibilities of transdimensional inversion for estimating deep Earth velocity and mantle structure, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-14782, https://doi.org/10.5194/egusphere-egu21-14782, 2021.

Realistic modelling of tightly coupled hydro-geomechanical processes is relevant for the assessment of many hydrological and geotechnical applications. Such processes occur in geological formations and are influenced by natural heterogeneities. Current numerical libraries offer capabilities and physics coupling, that have proven to be valuable in simulating various applications in geotechnical fields such as underground gas storage, rock fracturing, land subsidence and Earth resources extraction. However, implementation and verification of full heterogeneity of subsurface properties using high resolution field data in coupled simulations has not been done yet. Hence, we develop, verify and document RHEA (Real HEterogeneity App), an open-source fully coupled finite element application capable of including node-resolution hydro-geomechanical properties in coupled simulations. We propose a simple, yet powerful workflow to allow the integration of fully distributed hydro-geomechanical properties. We then verify the code with analytical solutions in one and two dimensions and propose a benchmark semi-analytical problem to verify heterogeneous systems with sharp gradients. Finally, we exemplify RHEA's capabilities with a comprehensive example integrating realistic properties. With this we demonstrate that RHEA is a verified open-source application able to integrate complex geology to perform scalable fully coupled hydro-geomechanical simulations. Our work is a valuable tool to assess real world hydro-geomechanical challenging systems that may include different levels of complexity like heterogeneous geology with several time and spatial scales and sharp gradients produced by contrasting subsurface properties.

How to cite: Bastias, J., Rau, G., Wilkins, A., and Blum, P.: RHEA: A verified simulator for hydro-geomechanical heterogeneity, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-3440, https://doi.org/10.5194/egusphere-egu21-3440, 2021.