GD10.1 | Advances in Numerical Modelling of Geological Processes: Methods and Applications
Advances in Numerical Modelling of Geological Processes: Methods and Applications
Co-organized by TS8
Convener: Ludovic Räss | Co-conveners: Boris Kaus, Ivan Utkin, Dave May, Thibault Duretz
Orals
| Tue, 25 Apr, 16:15–18:00 (CEST)
 
Room N1
Posters on site
| Attendance Tue, 25 Apr, 08:30–10:15 (CEST)
 
Hall X2
Orals |
Tue, 16:15
Posters on site
Tue, 08:30
Geological, geophysical, and geotechnical data delivers insights about the physical processes governing the Earth’s evolution. Typically, the data ranges from the internal structure of the Earth, plate kinematics, composition of geomaterials, estimation of physical conditions, dating of key geological events, thermal state of the Earth to more shallow processes such as natural and “engineered” reservoir dynamics and waste sequestration in the subsurface.

Combining such data with process-based numerical models provides deeper understanding of the dynamical Earth. Process-based models are powerful tools to predict the evolution of complex natural systems resolving the feedback among various physical processes. Integrating high-quality data into numerical simulations leads to a constructive workflow to further constrain the key parameters within the models. Innovative inversion strategies, linking forward dynamic models with observables, is therefore an important research topic that will improve our knowledge of the governing physical parameters.

The complexity of geological systems arises from their multi-physics nature, as these systems combine hydrological, thermal, chemical and mechanical processes (e.g. Earth’s mantle convection). 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 tools to capture spontaneous localisation and represents a challenging though essential research direction.

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
- AI / Machine learning-based approaches
- code and methodology comparisons (“benchmarks”)
- open source implementations for the community

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

Orals: Tue, 25 Apr | Room N1

Chairpersons: Ludovic Räss, Ivan Utkin, Boris Kaus
16:15–16:18
Introduction
Methods
16:18–16:28
|
EGU23-15888
|
On-site presentation
|
­­­HPC implementation and algorithm of the Thermal Lattice Boltzmann Method for geodynamics simulations in 2D and 3D 
Peter Mora, Gabriele Morra, and Dave Yuen

The Thermal Lattice Boltzmann Method (TLBM) is a powerful numerical method for thermally driven fluid flow simulations that is starting to be applied to geodynamics research. It is based on solving the Boltzmann equations on a discrete lattice and involves two steps of movement and collision of particle number densities carrying mass density and energy density on a discrete lattice. The collision step is achieved by relaxing the distributions to the equilibrium distribution where the relaxation times relate to the kinematic viscosity and thermal diffusivity of the fluid. We present the TLBM algorithm and an optimized HPC implementation of the TLBM where the main code is written in python using MPI for python, and this code calls highly optimized c functions for the kernels which do the heavy computational work. The same code works in 2D or 3D and we calculate the optimal 2D or 3D domain decomposition at the start of each run. Edges of domains are sent and received using optimal unblocking MPI requests, with the send and receive requests and buffers initialized at the start of a run to further optimize the communication costs. We present performance results which show near linear speedup to thousands of cores provided the domain size is not too small. We achieve of order 2-3 Gflops per core which is typically over 50% of peak performance. We show 2D runs using a highly nonlinear rheology which promotes the formation of plate-tectonic like dynamics with upwelling and downwelling plumes with the horizontal motion tending to be constrained to the upper 100km of the model. We also show 2D and 3D runs with temperature dependent viscosities and power law thermal boundary layer scaling with Nusselt number. And we show runs of simulations with high Rayleigh numbers up to 10**12 and Prandtl numbers up to 10**4. The TLBM offers a means to study the effect of highly nonlinear rheologies on geodynamical processes, and may eventually lead to a more complete simulation capability for studies of planet and exoplanet evolution.

How to cite: Mora, P., Morra, G., and Yuen, D.: ­­­HPC implementation and algorithm of the Thermal Lattice Boltzmann Method for geodynamics simulations in 2D and 3D, EGU General Assembly 2023, Vienna, Austria, 23–28 Apr 2023, EGU23-15888, https://doi.org/10.5194/egusphere-egu23-15888, 2023.

16:28–16:38
|
EGU23-13391
|
ECS
|
On-site presentation
High-performance Computing in modeling of Landslide Post-failure Stage using Material Point Method 
Zenan Huo, Yury Alkhimenkov, Flavio Calvo, Marc-Henri Derron, Michel Jaboyedoff, Yury Podladchikov, and Emmanuel Wyser

The post-failure of landslide is a stage where large deformations are present. It is difficult to properly resolve such large deformations using traditional mesh-based numerical methods. Meshless methods, such as the material point method (MPM), can resolve such problems by reducing the dependence on the mesh. However, the time-consuming mapping procedure between the material points and background nodes exists at each time step of MPM, consequently, one needs an efficient implementation taking advantage of modern computer hardware architectures for a high-resolution computational model. In the present study, we develop a high-performance MPM simulation package using Julia language to simulate the landslide post-failure stage. We show both the 2D and 3D computation models. The parallel algorithm on the GPU version is based on the features of MPM through CUDA.jl, a library that natively supports CUDA computing in Julia. To validate the performance of the present simulation package, we perform benchmarks on both CPU and GPU versions of the package. Furthermore, we use the uniform Generalized Interpolation MPM (uGIMP) and apply it to resolve a real problem to demonstrate the capabilities of this package.  The simulation result is in good agreement with the ground truth. HPC simulation is not only reproducing the run-out process but also provides us with a better understanding of the complex mechanisms involved in landslide movements.

How to cite: Huo, Z., Alkhimenkov, Y., Calvo, F., Derron, M.-H., Jaboyedoff, M., Podladchikov, Y., and Wyser, E.: High-performance Computing in modeling of Landslide Post-failure Stage using Material Point Method, EGU General Assembly 2023, Vienna, Austria, 23–28 Apr 2023, EGU23-13391, https://doi.org/10.5194/egusphere-egu23-13391, 2023.

16:38–16:48
|
EGU23-3521
|
ECS
|
On-site presentation
DEM crack propagation using a FEM-DEM bridging coupling 
Manon Voisin-Leprince, Joaquin Garcia-Suarez, Guillaume Anciaux, and Jean-François Molinari

The behavior of seismic faults depends on the response of the discrete microconstituents trapped in the region between continuum masses, which is usually termed “gouge”. The gouge is a particle region composed of amorphous grains. Conversely, the regions surrounding the gouge can be conceptualized as continua. The study of such system dynamics (slip) requires the understanding of several scales, from particle size to meter scale and above, to properly account for loading conditions. Our final objective in this study is to assess to what extent we can understand friction by leveraging an analogy to fracture. Dynamic friction between sliding surfaces resembles a dynamic mode-II crack, but this equivalence is brought into question when granularity at the interface is considered. Based on the theory of linear-elastic fracture mechanics (LEFM), a stress concentration should be observed at the rupture front if indeed friction can be modeled with the toolkit of LEFM.

Simulating this system numerically remains a challenge, as, in order to capture proper physics, both the continuum and discrete aspects of the system must be harmoniously incorporated and coupled into a single model. An energy-based coupling strategy between the Finite Element Method (FEM), used to resolve the continuum portions, and the Discrete Element Method (DEM), to model the granularity of the interface, is used [2]. In this exploratory study, we begin by modeling a medium with strong inter-granular cohesion [1]. The use of the coupling ensures a large enough effective domain to control nicely the crack propagation.  The linear-elastic properties of both DEM and FEM portions are therefore matched to avoid wave reflections. Both mode-I and mode-II cracks are considered.

How to cite: Voisin-Leprince, M., Garcia-Suarez, J., Anciaux, G., and Molinari, J.-F.: DEM crack propagation using a FEM-DEM bridging coupling, EGU General Assembly 2023, Vienna, Austria, 23–28 Apr 2023, EGU23-3521, https://doi.org/10.5194/egusphere-egu23-3521, 2023.

16:48–16:58
|
EGU23-1613
|
ECS
|
On-site presentation
Studying the along fault variability of slow earthquake characteristic by modeling a combined viscoelastic and damage rheology 
Sina Massoumi, Véronique Dansereau, Jérôme Weiss, and Nikolai Shapiro

The seismo-tectonic cycle in the subduction zones is largely controlled by the level of coupling between the sliding oceanic and continental plates that strongly varies with depth. Close to the surface, at depths of a few tens of kilometers, the plate interface remains most of time locked and is occasionally broken by large earthquakes. On the other hand, the oceanic lithosphere slips into the mantle continually at large depths. Between these two zones of locked and stable slip, the transient zone is characterized by “slow earthquakes” that are mainly manifested by episodes of silent slip and tectonic tremor that are to some degree correlated in time.

 

Along-fault changes of the degree of inter-plate coupling are controlled by variations of the fault-zone rheology, which in turn is related by depth-dependent thermo-mechanical conditions and composition of rocks. The brittle-ductile transition and the slow earthquake cycle are often modeled with using the rate-and-state interface rheology. This empirical formulation represents the transition segment by assimilating brittle and frictional processes to the problem of a material interface friction. To this aim, a parametric model is obtained based on experimental studies of the frictional behavior of various materials at the laboratory scale. Although this framework reproduces the transition between a stick-slip cycle and the stable sliding behaviors it cannot represent steady-state relaxation processes and presents a limit to which it can be enriched to include the chemical, mineralogical and hydro-mechanical processes within faults.

 

To overcome these limitations, we are using a modeling approach based on a continuum volumetric rheology that allow us to model physically-based variations of parameters with depth. Namely, we use a combination of viscoelastic Maxwell and damage rheologies. The resulting model is capable to take into account the localized deformations associated with quasi-brittle processes on short time scales as well as the diffuse deformations associated with the stress relaxation in the bulk of the geophysical system over long time scales. The problem is studied in two dimensions with associated boundary conditions. Along fault variations of the important controlling parameters such as the viscosity, the cohesion coefficients, and the damage recovery time are investigated in order to understand their respective contribution in the slow earthquake cycle.

How to cite: Massoumi, S., Dansereau, V., Weiss, J., and Shapiro, N.: Studying the along fault variability of slow earthquake characteristic by modeling a combined viscoelastic and damage rheology, EGU General Assembly 2023, Vienna, Austria, 23–28 Apr 2023, EGU23-1613, https://doi.org/10.5194/egusphere-egu23-1613, 2023.

16:58–17:08
|
EGU23-6105
|
On-site presentation
Three-field finite-element modelling of coupled two-phase flow for geological problems: Towards the zero-porosity limit 
Gang Lu, Dave May, and Ritske Huismans

Two-phase flow, a system where Stokes flow and Darcy flow are coupled, is of great importance in the Earth's interior, such as in subduction zones, mid-ocean ridges, and hotspots. However, it remains challenging to solve the two-phase equations accurately in the zero-porosity limit, for example when melt is fully frozen below solidus temperature. Here we propose a new three-field formulation of the two-phase system and present a robust finite-element implementation, which can successfully solve for the system where zero and non-zero porosity domains are both present. The reformulated equations, with solid velocity (vs), total pressure (Pt), and fluid pressure (Pf) as unknowns, include penalty and regularization to avoid singularities, which exactly recover to the standard single-phase Stokes with penalty at zero porosity. The new formulation is implemented using a 2-D finite-element discretization with Q1P0Q1 elements. We demonstrate the correctness of our implementation based on benchmarks against analytical solutions, which gives expected convergence rates in both space and time. Example experiments, such as self-compaction, falling block, and mid-ocean ridge spreading, show that this formulation can robustly resolve zero- and non-zero-porosity domains simultaneously, and be used for a large range of applications in various geodynamic settings.

How to cite: Lu, G., May, D., and Huismans, R.: Three-field finite-element modelling of coupled two-phase flow for geological problems: Towards the zero-porosity limit, EGU General Assembly 2023, Vienna, Austria, 23–28 Apr 2023, EGU23-6105, https://doi.org/10.5194/egusphere-egu23-6105, 2023.

Applications
17:08–17:18
|
EGU23-17169
|
On-site presentation
Modeling focused fluid flow with geological heterogeneity 
Lawrence Hongliang Wang and Viktoriya Yarushina

Two-phase flow equations that couple solid deformation and fluid migration have opened new research trends in geodynamical simulations and modeling of subsurface engineering operations. A numerical model based on two-phase flow equations has been used to study the formation of focused fluid flow in ductile/plastic rocks. While the effects of material properties such as permeability, bulk viscosity, shear viscosity, and bulk moduli have been studied with simple models that contain mainly homogenous material, realistic models with geological heterogeneity are scarce. This is partly due to the physical nonlinearity of fluid-rock systems and the strong coupling between flow and deformation. Here we present numerical models with a viscoelastic approach that solves hydromechanics coupling using an efficient pseudo-transient solver, which can model focused fluid flow with sharp material boundaries. First, we study the effects of a less permeable block on the propagation of channelized fluid flow by varying the permeability factor by several orders of magnitude and block size. We found that an obstacle does not stop the propagation of the localized channels but deflects and slows them down. A wide block allows channels to pass through slowly, while a narrow block deflects the channels to the sides.  Second, we study the dynamics of fluid channels reaching a sharp geological boundary that is significantly less permeable. We also adjust the bulk viscosity and permeability exponent for different materials in our models to mimic the real geological materials. This makes it possible to consider more realistic scenarios with intraformational and top-sealing layers relevant to CO2 storage and natural fluid migration.

How to cite: Wang, L. H. and Yarushina, V.: Modeling focused fluid flow with geological heterogeneity, EGU General Assembly 2023, Vienna, Austria, 23–28 Apr 2023, EGU23-17169, https://doi.org/10.5194/egusphere-egu23-17169, 2023.

17:18–17:28
|
EGU23-11019
|
ECS
|
Virtual presentation
Numerical simulation of fluid pressure build-up below the seismic lid: Implications to ‘fault-valve’ mechanism for lode-type gold deposits 
Swarnendu Bhuyan and Mruganka Kumar Panigrahi

Orogenic gold systems are flow-controlled thermodynamic systems and typically occur in mid- to upper crustal environments where there is a strong coupling of deformation and fluid flow with attendant heat transfer and chemical reactions. Fluids are generated during metamorphic devolatilization reactions under greenschist to amphibolite/granulite facies conditions and accumulated below 15 km depth from the earth's surface due to the presence of an impermeable layer (‘seismic lid’) that prevents the upward flow of fluid. Here the fluid pressure regime ranges from hydrostatic above the seismic lid to suprahydrostatic value below the seismic lid. The formation of orogenic gold deposits is associated with fluid pressure variation and rupture of the fault. The ‘fault valve’ mechanism that operates during periodic seismic pumping is widely believed to be responsible for gold mineralization in such systems. A 2D model is generated with the help of COMSOL Multiphysics software to describe the fluctuations of fluid pressure based on fault-valve vis-à-vis seismic pumping mechanism taking various assumptions, standard physical and lithological parameters, and governing equations related to Darcy’s law, storage coefficient equation. The rectangular cross-section covers a region of 50 km long and 25 km deep. The 200 m width seismic lid is located at 15 km depth and a fault of infinite length along a strike of 300m width cutting through the seismic lid. Specific sense of movement on the fault and the inferences of tectonic movements are not considered for this study. The intact rocks have low porosity and low permeability and a fixed heat flux is assigned for the bottom boundary and other boundaries are thermally insulated. Based on results obtained from the numerical simulation, the followings can be concluded. (1) The fluctuation of the fluid pressure shows a larger variation below the seismic lid and the zone where the fault penetrates the seismic lid (Fig 1). (2) A high-angle fault seems favourable for fluid flow and may not give rise buildup of supralithostatic fluid pressure that is essential for fault-valve process to operate. On the other hand, orogenic gold deposits are hosted in high angle reverse faults/ shear zones. Therefore the operation of Sibson’s cycle for the origin of lode-type gold deposit is need to be more critically evaluated.

                                                                                                 

                                                      

 

                                                                               

How to cite: Bhuyan, S. and Panigrahi, M. K.: Numerical simulation of fluid pressure build-up below the seismic lid: Implications to ‘fault-valve’ mechanism for lode-type gold deposits, EGU General Assembly 2023, Vienna, Austria, 23–28 Apr 2023, EGU23-11019, https://doi.org/10.5194/egusphere-egu23-11019, 2023.

17:28–17:38
|
EGU23-12057
|
ECS
|
On-site presentation
Self-consistent Reconstructions of the Earth's Mantle in Space and Time using Nonlinear Rheologies 
Sia Ghelichkhan and Rhodri Davies

Reconstructing the spatial and temporal evolution of Earth’s mantle through the recent geological past stands as one of the grand challenges in Geodynamics. One method to invert for the mantle’s evolution is to reformulate mantle flow as an optimisation problem using the adjoint method, where uncertain properties, such as the mantle’s previous thermo-chemical states, are found by minimising a misfit functional that represents the difference between model predictions and geodynamic inferences from various disciplines, including seismology, geodesy, and geochemistry. While the rapid growth in high-performance computing capacities has underpinned an ever-growing number of such reconstruction models, they often make several simplifying physical assumptions, or are limited in the number of assimilated datasets, thus limiting their applicability.

Here we present our latest attempts at reconstructing the evolution of Earth’s mantle using complex non-linear rheologies. Our approach builds upon a novel algorithmic differentiation method as implemented in dolfin-adjoint, together with state-of-the-art optimisation methods, developed using the Rapid Optimisation Library. Using analytical and synthetic examples, we show that the self-consistent derivation of the adjoint equations in our approach provides a pathway for accurate inversions for past-mantle flow.

How to cite: Ghelichkhan, S. and Davies, R.: Self-consistent Reconstructions of the Earth's Mantle in Space and Time using Nonlinear Rheologies, EGU General Assembly 2023, Vienna, Austria, 23–28 Apr 2023, EGU23-12057, https://doi.org/10.5194/egusphere-egu23-12057, 2023.

17:38–17:48
|
EGU23-12654
|
ECS
|
On-site presentation
Shell vs. plate tectonics: numerical stress quantification in a shortening lithosphere with strain localization 
William R. Halter, Roman Kulakov, Thibault Duretz, and Stefan M. Schmalholz

The mechanical characteristics of a shell, having a double curvature, are fundamentally different to the characteristics of a plate, having no curvature in its undeformed state. Geometrically, the Earth’s lithosphere is a shell rather than a plate. However, most geodynamic numerical models applied to study the deformation of the lithosphere do not consider this curvature. It is currently unclear whether the shell-type geometry of the lithosphere has a significant impact on lithosphere deformation on the scale of few 1000 kilometers. This study investigates the importance of considering lithospheric shells and compares numerical results of a shortening shell-type and plate-type lithosphere. We apply the two-dimensional state-of-the-art thermo-mechanical code MDoodz (Duretz et al. 2021). We consider a shortening lithosphere in an initially curved and in an initially rectangular geometry and calculate the spatio-temporal stress distribution inside the deforming lithosphere. We further present preliminary results on the effects and relative importance of various softening mechanism, leading to strain localization and subduction initiation, such as thermal softening, grain size reduction, or anisotropy generation due to fabric development.

 

REFERENCES

Duretz T., R. de Borst and P. Yamato (2021), Modeling Lithospheric Deformation Using a Compressible Visco-Elasto-Viscoplastic Rheology and the Effective Viscosity Approach, Geochemistry, Geophysics, Geosystems, Vol. 22 (8), e2021GC009675

How to cite: Halter, W. R., Kulakov, R., Duretz, T., and Schmalholz, S. M.: Shell vs. plate tectonics: numerical stress quantification in a shortening lithosphere with strain localization, EGU General Assembly 2023, Vienna, Austria, 23–28 Apr 2023, EGU23-12654, https://doi.org/10.5194/egusphere-egu23-12654, 2023.

17:48–17:58
|
EGU23-11502
|
ECS
|
On-site presentation
Quantifying Geodynamical Influences through Physics-Based Machine Learning: A Case Study from the Alpine Region 
Denise Degen, Ajay Kumar, Mauro Cacace, Magdalena Scheck-Wenderoth, and Florian Wellmann

Characterizing the influence of geodynamical models is important to improve our understanding of the development and current state of subsurface properties. Which are, in turn, of great societal relevance, for questions such as renewable energy. However, enabling a quantifiable characterization is a major challenge in Geodynamics, due to the high computational cost associated with both the model and the analysis for characterizing the influential parameters. The high cost of the model is caused by a high dimensionality in space, and time and a large number of input parameters. The cost of the probabilistic analyses is related to the large number of individual model solves required for performing the characterization.

To address this computational challenge, we employ the non-intrusive RB method, which combines advanced mathematical algorithms and novel machine learning methods. The method produces models that considerably reduce the dimensionality, yielding an acceleration of several orders of magnitude while maintaining the physical principles. In contrast, to other machine learning methods, the non-intrusive RB method produces explainable models, which is a crucial property for later analyses and predictions.

In this work, we demonstrate how the methodology can be beneficially used for the construction of reliable surrogate models of large-scale geodynamical applications without impacting the underlying physics. Furthermore, we show the benefits of global variance-based sensitivity analysis to quantifiable characterize the influence of the densities and viscosities on both the topography and velocity for the designated case study of the Alpine Region. We employ a global sensitivity analysis to account for possible parameter correlations and nonlinearities.

How to cite: Degen, D., Kumar, A., Cacace, M., Scheck-Wenderoth, M., and Wellmann, F.: Quantifying Geodynamical Influences through Physics-Based Machine Learning: A Case Study from the Alpine Region, EGU General Assembly 2023, Vienna, Austria, 23–28 Apr 2023, EGU23-11502, https://doi.org/10.5194/egusphere-egu23-11502, 2023.

17:58–18:00
Outlook

Posters on site: Tue, 25 Apr, 08:30–10:15 | Hall X2

Chairpersons: Dave May, Thibault Duretz, Boris Kaus
X2.224
|
EGU23-9242
|
solicited
|
Highlight
Using Julia for the next generation of HPC software for geodynamic modelling 
Albert de Montserrat Navarro, Boris Kaus, Ludovic Räss, Ivan Utkin, and Paul Tackley

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

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

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

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

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

X2.225
|
EGU23-1010
|
ECS
|
Global scale numerical geodynamic modelling with a free surface using a volume of fluid method 
Timothy Gray, Paul Tackley, and Taras Gerya

The study of coupled Earth systems, and in particular the coupled interactions between the lithosphere, atmosphere, and biosphere, have received greater attention in recent years (Gerya et al. 2020). Interactions between these systems occur primarily at the surface, and are driven on the large scale by topographic and bathymetric evolution controlled by deep mantle processes. However, due to the large difference in length scales between the mantle and the surface, it is difficult to capture topographic evolution to a high degree of accuracy in existing global mantle convection models including a free surface boundary condition.

Global mantle convection models incorporating a free surface often employ a marker-in-cell technique with a layer of “sticky air” (i.e. material with the density of the air and sufficiently low viscosity, which is still much higher than that of real air) to characterise the surface. However, accurate topographic evolution using this method requires a high density of markers near the surface. This need for additional computational resources motivates alternative methods of tracking the interface between the air and rock layers, as is done frequently in existing multiphase fluid flow codes. A volume of fluid method with piecewise-linear interface reconstruction provides a suitable method for tracking a surface in a performant way with the sub-grid level topographic resolution that is necessary for coupling global scale geodynamic models to models of other Earth systems.

We demonstrate benchmarks of an implementation of a volume of fluid method within the existing advanced mantle convection code StagYY (Tackley, 2008). Our method is applicable to both 2D and 3D geometries, and on both Cartesian and non-Cartesian grids. Models of global scale topography and evolution produced using StagYY may later be used as a tool for further studies on the coupling of mantle dynamics with modelling of the landscape, and the evolution of the atmosphere and biosphere.

How to cite: Gray, T., Tackley, P., and Gerya, T.: Global scale numerical geodynamic modelling with a free surface using a volume of fluid method, EGU General Assembly 2023, Vienna, Austria, 23–28 Apr 2023, EGU23-1010, https://doi.org/10.5194/egusphere-egu23-1010, 2023.

X2.226
|
EGU23-1845
|
ECS
|
On the statement and numerical solution of the thermal problem within inversion methods for the study of lithospheric structure. 
Mariano Tomás Fernandez, Sergio Zlotnik, and Pedro Díez

One of the goals of geophysicists is mapping and understanding the current structure of the Earth including its variations in composition, temperature and dynamical state. This structure is only accessible via indirect observations and, therefore, the mathematical problem to be solved is of an inverse kind. Within the inverse solver, many forward problems will be tested until finding a configuration compatible with the observations. This work deals with the problem statement and numerical solution of the forward thermal problem that arises from an inverse solver. In this case, we will use a simple parameterization of the Lithosphere-Asthenosphere Boundary (LAB), but the results are useful for other parametric description (e.g. one parameter per each cell). 
A simplified model is used to show the ill-posedness of the mathematical problem arising when the LAB --an isotherm whose location is determined by the input parameters-- is imposed within the domain, over-constraining the forward problem. This is well-known in the community and several authors have proposed different approaches to circumvent it. Nevertheless, the strategies used in practice usually involve some non-physical procedures such as transitional regions where two different temperature fields are made compatible by smearing out differences. Generally, the solution in these regions does not comply with the governing equation and exhibits a non-physical behaviour. 
In this work, we propose a specific problem statement for the temperature with interior essential conditions. The resulting problem is mathematically sound and results in a two-step numerical solver. This guarantees a self-consistent temperature field, in the sense that it respects the thermal governing equations everywhere. 
The numerical domain is divided into two subdomains (lithosphere and asthenosphere) that are solved separately in the same mesh, using an unfitted mesh methodology. First, the temperature of the lithosphere is computed using the essential condition on the LAB. Second, the temperature in the mantle is obtained by minimizing a residual that measures the compatibility between the two subdomains in terms of LAB temperatures and across-LAB fluxes. This is done by adjusting the proper fluxes at the bottom of the numerical domain. 
Several examples are presented showing that the obtained temperature fields are stable and oscillation-free. Moreover, the resulting fluxes at the bottom of the domain are reasonable and compatible with the expected values.

How to cite: Fernandez, M. T., Zlotnik, S., and Díez, P.: On the statement and numerical solution of the thermal problem within inversion methods for the study of lithospheric structure., EGU General Assembly 2023, Vienna, Austria, 23–28 Apr 2023, EGU23-1845, https://doi.org/10.5194/egusphere-egu23-1845, 2023.

X2.227
|
EGU23-9515
|
ECS
Seismic response to volcanic processes at Mount Etna: coupling thermomechanical simulations with seismic wave-equation modeling 
Michelle Bensing, Sergio Vinciguerra, and Luca De Siena

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

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

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

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

X2.228
|
EGU23-11594
|
ECS
GPU-based numerical models of rapid ductile strain localization due to thermal runaway 
Arne Spang, Marcel Thielmann, and Daniel Kiss

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

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

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

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

X2.229
|
EGU23-12450
|
ECS
GPU-based finite-difference solution of 3-D stress distribution around continental plateaus in spherical coordinates 
Emilie Macherel, Yuri Podladchikov, Ludovic Räss, and Stefan M. Schmalholz

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

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

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

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

X2.230
|
EGU23-13828
An unconstrained formulation for thermodynamic complex solution phase minimization 
Nicolas Riel, Boris Kaus, and Eleanor Green

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

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

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

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

X2.231
|
EGU23-14864
|
ECS
Generalization of the Nitsche method to apply oblique boundary conditions in regional geodynamic models 
Anthony Jourdon, Dave A. May, and Alice A. Gabriel

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

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

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

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

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

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

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

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

 

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

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

X2.232
|
EGU23-17307
Modelling lithosphere deformation with non-linear anisotropic constitutive models 
Roman Kulakov, William Halter, Stefan Schmalholz, and Thibault Duretz

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

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

X2.233
|
EGU23-14034
Towards integrated numerical models of lithospheric-scale magmatic systems 
Boris Kaus, Daniel Kiss, Albert de Montserrat, Nicolas Riel, Nicolas Berlie, and Arne Spang

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

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

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

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

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