Power-law viscous flow describes the first-order features of long-term lithosphere deformation. Due to the ellipticity of the Earth, the lithosphere is mechanically analogous to a shell, characterized by a double curvature. The mechanical characteristics of a shell are fundamentally different to the characteristics of plates, having no curvature in their undeformed state. The systematic quantification of the magnitude and the spatiotemporal distribution of strain, strain-rate and stress inside a deforming lithospheric shell is thus of major importance: stress is, for example, a key physical quantity that controls geodynamic processes such as metamorphic reactions, decompression melting, lithospheric flexure, subduction initiation or earthquakes. Calculating these stresses in a three-dimensional (3-D), geometrically and mechanically heterogeneous lithosphere requires high-resolution and high-performance computing.

Here, we present numerical simulations of 3-D power-law viscous flow. We employ the pseudo-transient finite difference (PTFD) method, which enables efficient simulations of high-resolution 3-D deformation processes by implementing an iterative implicit solution strategy of the governing equations. The main challenges for the PTFD method are to guarantee convergence, minimize the required iteration count and speed-up the iterations. We implemented the PTFD algorithm using the Julia language (julialang.org) to enable optimal parallel execution on multiple CPUs and GPUs using the ParallelStencil.jl module (https://github.com/omlins/ParallelStencil.jl). ParallelStencil.jl enables execution on multi-threaded CPUs and Nvidia GPUs using a single switch.

We present PTFD simulations of mechanically heterogeneous (weak and less dense spherical inclusion), incompressible 3-D power-law viscous flow under gravity in cartesian, cylindrical and spherical coordinates systems. The viscous flow is described by a linear combination of a linear viscous and a power-law viscous flow law, representing diffusion and dislocation creep, respectively. The iterative solution strategy builds upon pseudo-viscoelastic behavior to minimize the iteration count by exploiting the fundamental characteristics of viscoelastic wave propagation. We performed systematic numerical simulations to investigate the impact of (i) buoyancy versus shear forces and (ii) linear versus power-law viscous flow on the vertical velocity of the spherical inclusion under bulk strike-slip shearing. We report the systematic results using the controlling dimensionless numbers and compare the numerical results with analytical predictions for buoyancy-driven flow of inclusions in a power-law matrix. We also aim to unveil preliminary results for a vertically and locally loaded power-law viscous lithosphere showing the impact of different lithosphere curvatures on the resulting stress field.

How to cite: Macherel, E., Podladchikov, Y., Räss, L., and Schmalholz, S. M.: GPU-based pseudo-transient finite difference solution for 3-D gravity- and shear-driven power-law viscous flow, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-8816, https://doi.org/10.5194/egusphere-egu22-8816, 2022.

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 (Keller et al., 2013). Here we present our progress in consistent numerical modelling of poro-viscoelastic-viscoplastic (VEVP) flow. We show that a VEVP model with a new, hyperbolic yield surface can help to robustly simulate both shear and tensile modes of plastic failure in a two-phase system. 

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 with shear failure only. In two-phase systems, tensile failure of rocks can also occur due to an overpressured liquid phase. 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.

An appropriate rheological model is required to meet this challenge. The most straightforward choice is a Maxwell visco-elasto-plastic model, but this leads to grid-scale localisation and hence mesh-dependence. To obtain mesh-independent shear localisation, we employ the visco-elasto-viscoplastic model by introducing a viscous dashpot in parallel to the plasticity element. Whilst this formulation has shown promise in regularising shear failures in a single-phase flow model (de Borst and Duretz, 2020), its incorporation within two-phase systems has not been examined. We will show that the linear Griffith criteria for the tensile failure can lead to convergence issues whereas a new, hyperbolic yield surface is proposed to resolve these numerical issues. This yield surface provides a smooth transition between the two modes of failure.

The underlying PDEs are discretised using a conservative, finite-difference, staggered-grid framework implemented with PETSc (FD-PDE) that supports single-/two-phase flow magma dynamics. Here, we present simplified model problems using the FD-PDE framework for poro-viscoelastic-viscoplastic models designed to characterise the solution quality and assess both the discretisation and solver robustness. It has been observed that employing the hyperbolic yield surface improved the robustness in simulating plastic failures in both modes.

References

Keller, T., May, D. A., & Kaus, B. J. P., (2013). Numerical modelling of magma dynamics coupled to tectonic deformation of lithosphere and crust, Geophysical Journal International, v195, 1406-1442, https://doi.org/10.1093/gji/ggt306.

de Borst, R., Duretz, T., (2020). On viscoplastic regularisation of strain-softening rocks and soils. International Journal for Numerical and Analytical Methods in Geomechanics, v44, 890-903. https://doi.org/10.1002/nag.3046.

How to cite: Li, Y., Pusok, A., May, D., and Katz, R.: Simulation of partially molten rocks with visco-elasto-viscoplastic rheology and a hyperbolic yield surface for plasticity, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-5594, https://doi.org/10.5194/egusphere-egu22-5594, 2022.

The formation and maintenance of narrow, lithospheric shear zones and their importance in plate-tectonics remain one of the major problems in geodynamics. While the cause and consequence of strain localization and weakening within the lithosphere remain debated, it is clear that these processes play an essential role in lithospheric deformation across a wide range of spatio-temporal scales. Here, we analyze the efficiency of strain localization in a 2-D visco-elasto-plastic medium for a strain-dependent weakening and healing (SDWH) rheology using 2-D numerical, thermo-mechanical experiments with kinematic boundary conditions. Such a parameterized rheology successfully mimics more complex transient weakening and healing processes, akin to a grain-size sensitive composite (diffusion and dislocation creep) rheology. In addition, the SDWH rheology allows for memory of deformation. This enables self-consistent formation and reactivation of inherited weak zones within the lithosphere and sustains those weak zones over an extended period of time. We further analyze the resulting shear zone patterns and seek to answer the questions: What is the typical, effective intensity of strain localization? What are the dimensions of the resulting shear zones? Are such shear zones mesh-dependent in numerical models and, if so, can we exploit existing regularization approaches for the SDWH rheology?

How to cite: Fuchs, L., Duretz, T., and Becker, T. W.: Strain localization in a visco-elasto-plastic medium using strain-dependent weakening and healing rheology, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-11469, https://doi.org/10.5194/egusphere-egu22-11469, 2022.

One of the great challenges involved in modelling the lithosphere is its plastic behaviour, especially when dealing with compressible materials. Shear fractures are designated as mode 2 and 3 and can be described using a Linear Mohr Coulomb envelope  or a simplification of it like Drucker-Prager. Meanwhile, mode 1 fractures are created when the normal stresses become tensile  and require another yield function, such as the Griffith criterion or a tension cap function.

While the governing equations are well known and widely employed in engineering codes, they are usually expressed with a displacement formulation. Most geodynamic codes, on the other hand, use pressure and velocity as their primary variables. A numerically robust method that takes all plasticity modes into account in a staggered finite difference discretization remains an open task. Here we present a composite yield function implemented with pressure-velocity formulation, capable of producing produce shear and tensile failure.

We have implemented this in a new code that employed PETSc through the recently updated PETSc.jl Julia interface, while utilizing the automatic differentiation tools in julia. We found this workflow to significantly reduce the development time of complex nonlinear coupled codes.  

We will describe the implementation, propose regularization schemes and discuss benchmark cases and simple applications. We demonstrate Newton convergence for most cases and will discuss different methods to combine multiple plastic flow laws.

How to cite: Berlie, N., Kaus, B., Popov, A., Kiss, D., and Riel, N.: How to break the lithosphere: a compressible pressure-velocity formulation for elasto-visco-plastic rheologies that includes shear and tensile failure with dilation, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-6566, https://doi.org/10.5194/egusphere-egu22-6566, 2022.

The rate- and state-dependent friction (RSF) laws (Dieterich, 1979; Ruina, 1983) have been widely successful in capturing the behavior of sliding surfaces in laboratory settings, as well as reproducing a range of natural fault slip phenomena in numerical models.

Studies of exhumed fault zones make it clear that faults are not two-dimensional features at the junction of two distinct bodies of rock, but instead evolve into complex damage zones that show clear signs of multi-scale fracturing, grain diminution, hydro-thermal effects and chemical and petrological changes. Many of these observed factors have been experimentally verified, and several studies have furthered our theoretical understanding of earthquakes and other seismic phenomena as volumetric, bulk-rock processes, including Sleep (1995, 1997), Lyakhovsky and Ben-Zion et al. (2011, 2014a,b, 2016), Niemeijer and Spiers et al. (2007, 2016, 2018), Roubicek (2014), and Barbot (2019).

While the established numerical modeling approach of simulating faults as planar features undergoing friction can be a useful and powerful homogenization of small-scale volumetric processes, there are also cases where this practice falls short -- most notably when studying faults that grow and evolve in response to a changing tectonic environment. This is mainly due to the computational challenges associated with automating the construction of a fault-resolving conformal mesh.

Motivated by this issue, we formulate a generalization of RSF as a plastic or viscous flow law with generation, diffusion, and healing of damage that gives rise to mathematically and numerically well-behaved finite shear bands that closely mimic the behavior of the original laboratory-derived formulation (Pranger et al., submitted). The proposed formulation includes the well-known RSF laws for an infinitely thin fault as a limit case as the damage diffusion length scale tends to zero. We will show the behavior of this new bulk RSF formulation with results of high-resolution 1D and 2D numerical simulations.

Dieterich, J.H. (1979), J. Geophys. Res., 84 (B5), 2161.
Ruina, A. (1983), JGR: Solid Earth, 88 (B12), 10359–10370.
Sleep, N.H. (1995), JGR, 100 (B7), 13065–13080.
Sleep, N.H. (1997), JGR: Solid Earth 102 (B2), 2875–2895.
Roubíček, T. (2014), GJI 199.1, 286–295.
Lyakhovsky, Hamiel and Ben-Zion (2011), J. Mech. Phys. Solids, 59, 1752-1776.
Lyakhovsky and Ben-Zion (2014a), PAGeoph 171.11, 3099–3123.
Lyakhovsky and Ben-Zion (2014b), J. Mech. Phys. Solids 64, 184–197.
Lyakhovsky, Ben-Zion et al. (2016), GJI 206.2, 1126–1143.
Barbot (2019), Tectonophysics 765, 129–145.
Niemeijer and Spiers (2007), JGR 112, B10405,
Chen and Spiers (2016), JGR: Solid Earth 121, 8642–8665.
van den Ende, Chen et al. (2018), Tectonophysics 733, 273-295.
Pranger et al. (202X), ESSOAr (https://www.essoar.org/doi/10.1002/essoar.10508569.1)

How to cite: Pranger, C., Sanan, P., May, D., Le Pourhiet, L., Räss, L., and Gabriel, A.: Rate and state friction on spontaneously evolving faults, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-10412, https://doi.org/10.5194/egusphere-egu22-10412, 2022.

The development of highly efficient, robust, and scalable numerical algorithms lags behind the rapid increase in massive parallelism of modern hardware. In this work, we address this challenge with the accelerated pseudo-transient iterative method. This method is motivated by the physical analogy between numerical iterations and transient processes converging to a steady state.

We analytically determine optimal iteration parameters for a variety of basic physical processes such as diffusion, diffusion-reaction and non-inertial viscous fluid flow featuring Maxwell viscoelastic rheology. We further confirm the validity of theoretical predictions with numerical experiments.

We provide an efficient numerical implementation of various pseudo-transient solvers on graphical processing units (GPUs) using the Julia language. We achieve a parallel efficiency over 96% on 2197 GPUs in distributed memory parallelisation weak scaling benchmarks. 2197 GPUs allow for unprecedented terascale solutions of 3D variable viscosity Stokes flow involving over 1.2 trillion degrees of freedom.

We verify the robustness of the method by handling contrasts up to 9 orders of magnitude in material parameters such as viscosity, and arbitrary distribution of viscous inclusions for different flow configurations. Moreover, we show that this method is well suited to tackle strongly nonlinear problems such as shear-banding in a visco-elasto-plastic medium.

We additionally motivate the accessibility of the method by its conciseness, flexibility, physically motivated derivation, and ease of implementation. This solution strategy has thus a great potential for future high-performance computing applications, and for paving the road to exascale in the geosciences and beyond.

How to cite: Utkin, I., Rass, L., Duretz, T., Omlin, S., and Podladchikov, Y.: Assessing the robustness and scalability of the accelerated pseudo-transient method towards exascale computing, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-9815, https://doi.org/10.5194/egusphere-egu22-9815, 2022.

Transient superstructures in mantle convection whose life and morphology vary with Rayleigh and Prandtl number have recently been demonstrated. These superstructures appear to be a two-scale phenomenon where smaller scale rolls organize into larger scale convection cells. Simulation of such superstructures requires the ability to model 3D convection in box with very large width/height ratio of order greater than 10, and with resolution to resolve the thermal boundary layer at Rayleigh numbers of 10⁸ to 10¹⁰, respectively at least 100 height levels and 200 height levels. We achieve this with an efficient parallel implementation of the Lattice Boltzmann Method using Python which operates with high efficiency and linear speedup on thousands of cores. We present simulations with Rayleigh numbers of up to 10¹⁰ and Prandtl numbers from 1 to 100 to illustrate covering regimes from a magma ocean to solid mantle convection. We further present simulations using the LBM to model variable viscosity – specifically, temperature dependent– and illustrate the existence of pulsating plumes. We further demonstrate power law scaling between Nusselt number and Rayleigh number Nu  ~ Ra^g, which to first order is consistent with the Grossmann and Lohse theory.

How to cite: Mora, P., Morra, G., and Yuen, D.: Simulation of 3D transient superstructures in mantle convection and variable viscosity via the Lattice Boltzmann Method, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-9069, https://doi.org/10.5194/egusphere-egu22-9069, 2022.

The core-mantle boundary (CMB) is the most prominent compositional boundary inside the Earth. Its topography provides insight on lower mantle flow and the thermochemical structure above the CMB. Yet, CMB topography remains challenging to observe and estimates from seismology vary substantially. Numerical models of mantle convection provide complementary means to estimate CMB topography. Classically, topography is determined from the normal stresses acting on the CMB. However, this is known to face severe complications when applied to the surface boundary of the mantle, leading to non-Earth-like topographic scales and a different style of subduction. A (quasi-)free surface yields more Earth-like predictions, but for the CMB this comparison has never been made.

Here, we compare CMB topography predicted from mantle convection modelling using different treatments of the CMB. Specifically, we test the role of a ‘sticky core’, a quasi-fluid approximation the core. We compare results predicted by different codes (with either sticky core or true free base) and compare to a simple analytical case. Also, we simulate the evolution of subduction and deep thermochemical provinces to compare the topography of the (quasi-)free CMB and the free-slip approach. Initial results indicate that the sticky core approach can reproduce CMB topography reasonably well, but has rather high computational cost (grid resolution, number of particles). In analogy to the sticky air at the surface, the viscosity contrast of the sticky core layer determines the quality of predicted topography, with larger contrasts (≥10³) leading to acceptable levels of artificial CMB topography. In dynamic flow cases with vigorous mantle convection, entrainment by plumes further complicates application of the sticky core, but can be tackled with an unmixing procedure. A true free base tends to better accuracy than the sticky core approach and avoids the problem with entrainment, but it also comes with additional computational costs as various forces at the CMB have to be taken into account.

How to cite: Rolf, T., Crameri, F., Heyn, B. H., and Thielmann, M.: Testing a (quasi-)free base for modelling core-mantle boundary topography, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-9133, https://doi.org/10.5194/egusphere-egu22-9133, 2022.

Motions within the Earth mantle and tectonics constitute a single self-organized system which is cooling the planet over its geological history. Since the end of the XXth century, models of mantle convection self-generating plate tectonic behavior have progressed to a state that makes them applicable to global tectonic problems. The possibility of combining geological and geophysical data with dynamic models to retrieve the recent history of mantle flow and tectonics becomes realistic. Therefore, it is a challenge to build inverse methods to study inverse and sensitivity problems in the Earth's mantle convection. We have automatically generated the tangent-linear and the adjoint source code from the StaggYY code (Tackley, Phys. Earth Planet. Int. 171, 7-18, 2008). The Fortran code of the model was translated to the corresponding derivative codes using TAF (Transformation of Algorithms in Fortran), source-to-source translator. All codes run in parallel mode, using MPI (Message Passing Interface). The economic taping strategy of TAF, including re-computations, and checkpointing, helps to keep the memory footprint of the adjoint code low and the performance high. We highlight some key features of the automatic differentiation, evaluate the performance of the adjoint code, and show first results from 2D and 3D sensitivity fields, focusing on the relationships between temperature in the mantle and tectonics. Ultimately the addjoint code shall be applied to inversion and assimilation problems using a bayesian framework.

How to cite: Coltice, N., Blessing, S., Giering, R., and Tackley, P.: Automatic generation of the adjoint of the StagYY mantle convection model, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-9232, https://doi.org/10.5194/egusphere-egu22-9232, 2022.

Density-driven flows in porous media are frequently encountered in natural systems and arise from the gravitational instabilities introduced by fluid density gradients. They have significant economic and environmental impacts, and numerical modelling is often used to predict the behaviour of these flows for risk assessment, reservoir characterisation or management. However, modelling density-driven flow in porous media is very challenging due to the nonlinear coupling between flow and transport equations, the large domains of interest and the wide range of time and space scales involved. Solving this type of problem numerically using a fixed mesh can be prohibitively expensive.  Here, we apply a dynamic mesh optimisation (DMO) technique along with a control-volume-finite element method to simulate density-driven flows. DMO allows the mesh resolution and geometry to vary during a simulation to minimize an error metric for one or more solution fields of interest, refining where needed and coarsening elsewhere. We apply DMO to the Elder problem for several Rayleigh numbers. We demonstrate that DMO accurately reproduces the unique two-dimensional (2D) solutions for low Rayleigh number cases at significantly lower computational cost compared to an equivalent fixed mesh, with speedup of order x16. For unstable high Rayleigh number cases, multiple steady-state solutions exist, and we show that they are all captured by our approach with high accuracy and significantly reduced computational cost, with speedup of order x6. The lower computational cost of simulations using DMO allows extension of the high Rayleigh number case to a three-dimensional (3D) configuration and we demonstrate new steady-state solutions that have not been observed previously. Early-time, transient 3D patterns represent combinations of the previously observed, steady-state 2D solutions, but all evolve to a single, steady-state finger in the late time limit.

How to cite: Bahlali, M. L., Salinas, P., and Jackson, M. D.: Dynamic mesh optimisation for efficient numerical simulation of density-driven flows: Application to the 2- and 3-D Elder problem, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-12156, https://doi.org/10.5194/egusphere-egu22-12156, 2022.