Displays

Magma matters. From magmatism facilitating the differentiation of terrestrial planets into core, mantle and crust, to the magmatic activity that modulates plate tectonics and deep volatile cycles to maintain a habitable Earth, to volcanism that causes terrible hazards but also provides rich energy and mineral resources – igneous processes are integral to Earth and other planets. Our understanding of volcanoes and their deep magmatic roots derives from a range of disciplines including field geology, petrology and geochemistry, and geophysical imaging. Observational and experimental studies, however, are hampered by incomplete access to processes that play out across scales ranging from sub-micron size to thousands of kilometres, and from seconds to billions of years. Computational modelling provides tools for investigating igneous processes across these scales.

Over the past decade, my research has been focused on advancing the theoretical description and numerical application of multi-phase reaction–transport processes at the volcano to planetary scale. Mixture theory provides a framework to represent the spatially averaged behaviour of a large sample of microscopic phase constituents such as mineral grains, melt films, and vapour bubbles. This approach has been used successfully to model both porous flow of melt percolating through compacting partially molten rock, as well as suspension flow of crystals settling in convecting magma bodies. My recent work has introduced a new modelling framework to bridge the porous and suspension flow limits, and to extend beyound solid-liquid systems to multi-phase systems including several solid, liquid, and vapour phases. These advances provide new insights into the dynamics of crustal mush bodies, the outgassing and eruption of shallow magma reservoirs, and the generation of mineral resources by exsolution of exotic magmatic liquids.

How to cite: Keller, T.: Numerical modelling of igneous processes, EGU General Assembly 2020, Online, 4–8 May 2020, EGU2020-7447, https://doi.org/10.5194/egusphere-egu2020-7447, 2020.

We provide details on a novel formulation derived to describe the multiphysics controlling the deformation of porous rock under lithospheric conditions. The theory is developed consistent with the principles of thermodynamics and enables to capture the behaviour of porous rocks at the transition from frictional brittle behaviour to ductile viscous behaviour. It also accounts for the nonlinear feedback mechanisms derived from energetic consideration for the bi-phasic fluid-rock matrix system.

The formulation depicts a consistent, implicit visco-elasto-(visco)plastic rheology accounting for both a volumetric and a deviatoric response to applied loads, thereby avoiding the use of, the commonly assumed, plasticity limiter concept. The overstress plastic formulation introduces rate dependent mechanical behavior, an aspect that is consistent with experimental rock mechanics evidence and is also demonstrated to improve numerical stability when addressing problems related to plastic strain accumulation even in the absence of energetic feedbacks.

The introduction of a damage rheology permits to account for microstructural processes responsible for brittle-like material weakening and rate-dependent dissipative material behavior. The presence of a fluid phase is considered via a dynamic porosity, the evolution of which is demonstrated to primarily control the volumetric mechanical response of the stressed rock during incremental loading.

The above formulation has been integrated in a massively parallel, open source numerical framework with interfaces to state of the art HPC clusters. The results of a scalability and profile performance analysis on multi-core supercomputer are presented alongside with dedicated applications describing lithospheric rock deformation under different confining conditions as well as the bulk macroscopic material response recorded by laboratory experiments under triaxial conditions.

How to cite: Cacace, M. and Jacquey, A.: Thermodynamic consistent formulation for the multiphysics of a brittle ductile lithosphere - semi-brittle semi-ductile deformation and damage rheology, EGU General Assembly 2020, Online, 4–8 May 2020, EGU2020-13398, https://doi.org/10.5194/egusphere-egu2020-13398, 2020.

Scalable preconditioners for saddle point problems are essential to the solution of problems in geodynamics and beyond. Recent years have produced a wealth of research into efficient solvers for finite element methods.  These solvers are also effective, however, for orthogonal-grid finite volume discretizations of saddle point problems, also know as "staggered grid" or "marker and cell (MAC)" methods. Perhaps, ironically, due to the highly-structured nature of these discretizations, the use of advanced solvers is stymied due to the lack of a uniform topological abstraction, which is required for most scalable solvers, such as geometric multigrid.  We present new software to allow experimentation with and composition of these advanced solvers.  We focus on variable-viscosity Stokes problems with discontinuous coefficient jumps.  In particular, we attempt to demonstrate how the important know robust preconditioners may be employed, and how new variants may be experimented with.  Important solvers are compositions of block factorizations and multigrid cycles.  We demonstrate as many of these as possible, including triangular block preconditioners with nested multigrid solves, and monolithic multigrid solves with cellwise (Vanka) or field-based (Distributed Gauss-Seidel, Braess-Sarazin) smoothers.  Implementations are provided as part of the PETSc library, using the new DMStag component, and examples from the StagBL library are also shown where appropriate.  These tools are intended to help break down the barrier between cutting-edge research into advanced solvers (which is only becoming more complex, as multi-phase problems are further explored) and practical usage in geophysical research and production codes.

How to cite: Sanan, P. and May, D.: Composable, Scalable Solvers on Staggered Grids, EGU General Assembly 2020, Online, 4–8 May 2020, EGU2020-21413, https://doi.org/10.5194/egusphere-egu2020-21413, 2020.

Numerical simulations of the earthquake cycle have made great progress over the past decades to address important questions in earthquake physics and fault mechanics. However, significant challenges in bridging multiscale interactions between long-term tectonic deformation, aseismic fault slip, earthquake nucleation, and dynamic rupture still remain. In this study, we present results from GARNET, a newly-developed numerical library to simulate sequences of seismic and aseismic slip across scales. This finite difference code utilizes a fully staggered spatially adaptive rectilinear grid. Furthermore, it incorporates an automatic discretization algorithm and combines different physical ingredients, including a visco-elasto-plastic rheology and quasi- and fully dynamic formulation of inertial effects into one algorithm. While PETSc and Kokkos libraries are included for parallel computing, an adaptive time stepping is integrated into the algorithm to resolve both long- and short-time scales, ranging from years to milliseconds during the dynamic propagation of earthquake rupture.

Here we present results from two benchmarks (BP1 and BP3) based on the community code-verification effort for Sequences of Earthquakes and Aseismic Slip (SEAS) by the Southern California Earthquake Center (SCEC). BP1 benchmark is a 2D antiplane problem, with a 1D planar vertical strike-slip fault obeying rate-and-state friction, embedded in a 2D homogeneous, linear elastic half-space. The fault has a shallow seismogenic region with velocity-weakening friction and a deeper velocity-strengthening region, below which a relative plate motion rate is imposed. A periodic sequence of spontaneous, quasi-dynamic earthquakes and slow slip are simulated in the model. In the BP3 benchmark we consider full inertial effects during the dynamic rupture and we investigate its influence on earthquake behaviour and patterns. Results from these two benchmarks represent the first step towards more advanced seismic cycle models, which will help to enhance our understanding in earthquake physics.

How to cite: Dal Zilio, L., Li, M., van Dinther, Y., and Pranger, C.: Towards simulating sequences of seismic and aseismic slip across scales: Initial benchmarks and future directions, EGU General Assembly 2020, Online, 4–8 May 2020, EGU2020-17770, https://doi.org/10.5194/egusphere-egu2020-17770, 2020.

Ensemble based data assimilation approaches, such as the Ensemble Kalman Filter (EnKF), have been widely and successfully implemented to combine observations with dynamic models to investigate the evolution of a system’s state. Such inversions are powerful tools for providing forecasts as well as “hindcasting” events such as volcanic eruptions to investigate source parameters and triggering mechanisms. In this study, a high performance computing (HPC) adaptation of the EnKF is used to assimilate ground deformation observations from interferometric synthetic-aperture radar (InSAR) into high-fidelity, multiphysics finite element models to evaluate the prolonged unrest and June 26, 2018 eruption of Sierra Negra volcano, Galápagos. The stability of the Sierra Negra magma system is evaluated at each time step by estimating variations in reservoir overpressure, Mohr-Coulomb failure in the host rock, and tensile stress and failure along the reservoir boundary. The deformation of Sierra Negra is tracked over a decade, during which almost 5 meters of surface uplift has been recorded. The EnKF reveals that the evolution of the stress state in the host rock surrounding the Sierra Negra magma reservoir likely controlled the timing of the eruption. While increases in magma reservoir overpressure remained modest (< 10 MPa) throughout the data assimilation time period, significant Mohr-Coulomb failure is indicated in the lead up to the eruption coincident with increased seismicity along both trapdoor faults within Sierra Negra’s caldera and along the caldera’s ring faults. During the final stages of pre-eruptive unrest, the EnKF models indicate limited tensile failure, with no tensile failure along the northern portion of the magma system where the eruption commenced. Most strikingly, model calculations of significant through-going Mohr-Coulomb failure correspond in space and time with a Mw 5.4 earthquake recorded in the hours preceding the 2018 eruption. Subsequent stress modeling implicates the Mw 5.4 earthquake along the southern intra-caldera trapdoor fault as the potential catalyst for tensile failure and dike initiation along the reservoir to the north. In conclusion, the volcano EnKF approach successfully tracked the evolving stability of Sierra Negra, indicating great potential for future forecasting efforts.

How to cite: Gregg, P., Zhan, Y., Amelung, F., Albright, J., Geist, D., Mothes, P., Yunjun, Z., and Koric, S.: Linking thermomechanical models with geodetic observations to evaluate the 2018 eruption of Sierra Negra Volcano, Galápagos, EGU General Assembly 2020, Online, 4–8 May 2020, EGU2020-10835, https://doi.org/10.5194/egusphere-egu2020-10835, 2020.

This study investigates the effect of erroneous parameter values for state and parameter estimation using data assimilation. The numerical model chosen for this study solves the van der Pol equation, a second-order differential equation that can be used to simulate oscillatory processes, such as earthquakes. In the model, discrepancies in the parameter values can have a significant influence on the forecasted states of the model, which is even more significant if its behaviour is highly nonlinear. When observations of the state variables are assimilated to update the parameters along with the state variables, this improves the quality of the state forecasts. The results suggest that corrections in the model parameter not only recover the actual parameter values but also reduce state-variable errors after a certain time period. However, data assimilation that updates the state variables but not the parameter can lead to erroneous estimates as well as forecasts of the oscillation. Since the study is performed on a simplified nonlinear model framework, the consequences of these results for data assimilation in more realistic models remains to be investigated.

How to cite: Banerjee, A. and Vossepoel, F.: Considerations on Parameter and State Estimation with Ensemble Data Assimilation Methods – A Case Study with a Nonlinear Oscillator, EGU General Assembly 2020, Online, 4–8 May 2020, EGU2020-21702, https://doi.org/10.5194/egusphere-egu2020-21702, 2020.

Scientific computing related to solving partial differential equations (PDEs) frequently employ both real and complex numbers, consequently computation libraries such as PETSc provide support for these number types and their associated algebra. Other “exotic” number types such dual numbers, hyper-dual numbers and intervals are also highly desirable in the context of solving PDEs, however these are seldom available within HPC linear algebra and or discretisation libraries.

I’m this presentation I will summarise several exotic number types and discuss some of their potential uses when solving: linear PDEs, non-linear PDEs, discrete adjoint problems and PDE constrained optimisation problems. Using these as motivation, I will also describe how these exotic numbers can be supported within PETSc. The approach adopted is general (extensible), non invasive and allows users to select the particular number type representation at run-time. Moreover the general design is well suited to the characteristics of modern computational hardware and thus can efficiently exploit different forms of parallelism. Through a series of toy problems, the cost and performance of the exotic number implementations will be assessed.

How to cite: May, D. and Katz, R.: Abstract numbers and algebra in PETSc, EGU General Assembly 2020, Online, 4–8 May 2020, EGU2020-22323, https://doi.org/10.5194/egusphere-egu2020-22323, 2020.

This work is part of the "Subsidence" DeepNL project which aims to identify subsurface drivers of subsidence above the Groningen (the Netherlands) gas field and to forecast future subsidence. The hydrocarbon extraction in Groningen induces a pressure reduction in the gas reservoir which triggers compaction and land subsidence. This deep-subsurface process is modeled by a disc-shaped reservoir model, which is a superposition of individual nuclei of strain based on the Geertsma's approach. We estimate the surface deformation and the strength of the disc strain using a particle method. We apply the method to one single nucleus of strain at 3 km depth and extend to a disc-shape geometry. Synthetic experiments with a single nucleus of strain and with discs of varying sizes, 2.2 km to 13.3 km diameter, at 3 km depth are performed to assess the performance of the method for an increasing degree of complexity. Sequential Importance Resampling prevents the sample degeneracy when the number of nuclei increases. Adding a jitter noise in the resampling step avoids an impoverishment of the ensemble values. The results indicate that the method estimates the surface deformation and the strength for a large number of sources and for a relatively small effective ensemble size. In further investigations, localization can provide an additional means to deal with increasing dimensions and a relatively small ensemble size.

How to cite: Kim, S., Vossepoel, F., Hanssen, R. F., Wouters, M., Govers, R., and Stouthamer, E.: A particle method strategy to estimate subsidence induced by a high-dimensional disc-strain model for reservoir compaction., EGU General Assembly 2020, Online, 4–8 May 2020, EGU2020-22117, https://doi.org/10.5194/egusphere-egu2020-22117, 2020.

Plastic strain localization in a rate-independent limit is a physical phenomenon that demands robust and reliable nonlinear solves to be modeled in a computer. It is quite common on practice that standard iterative algorithms (e.g. Newton-Raphson), being applied to this demanding problem, lead to  convergence issues or even fail. From the physical viewpoint the problem can be attributed to the difference between the dilatation and friction angles for the pressure-sensitive materials, which gets worse as this difference gets larger.

Common remedies include deriving consistent Jacobian matrix, or switching to the equivalent rate-dependent visco-plastic formulation. Both methods have their specific side effects. A very high condition number of a heavily unsymmetrical Jacobian matrix renders it nearly useless in the context of an iterative linear solver, such as multigrid. Hence, the use of the Jacobian matrix is essentially limited to a 2D formulation, for which a direct solver is practical. The visco-plastic formulation is confusing form the conceptual viewpoint. It strives to achieve the convergence by modifying the physics of the problem. Hence, the stabilization viscosity is not a pure numerical parameter that can be freely selected, but it is a physical parameter that must be determined in the laboratory. The advantages of the visco-plastic formulation vanish, if rate-independent limit is considered, or if affordable grid size is (much) larger than the intrinsic localization length-scale. The latter condition is a dominant limiting factor for a 3D model. 

In this work we share a few recipes, that can potentially improve the convergence of the rate-independent plasticity problems, without relying on the availability of a direct solver, or perturbing the physics.

How to cite: Popov, A. and Reuber, G.: How to stabilize nonlinear solvers for rate-independent plasticity problems?, EGU General Assembly 2020, Online, 4–8 May 2020, EGU2020-8803, https://doi.org/10.5194/egusphere-egu2020-8803, 2020.

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. Due to mass conservation, thinning of the crust implies horizontal spreading of the plateau towards the lower altitude surroundings. This spreading is observed in GPS velocity records. For example, around the Tibetan Plateau, horizontal surface velocities are in the order of 2 cm/yr. Crustal flow also generates differential stresses in and around the plateau. Estimating these stresses and their spatial variation in 3-D is a computational challenge and necessitates new advanced computing methods.

Here, we present a new 3-D numerical algorithm to solve the Stokes equations under gravity. The algorithm is based on an Eulerian pseudo-transient finite difference method. To test the algorithm, we consider a simplified plateau geometry and density structure. Two different simulations are performed, one pseudo-2D simulation and one full 3-D simulation considering the corner region of a plateau. The pseudo-transient method allows an explicit solution of the Stokes equations. When the pseudo-transient time derivatives approach zero, a steady-state solution for the velocity field is obtained. In the first simulations, the rheology is linear viscous. The crust-mantle interface and the interface between crust and overlying sticky air are tracked with a Cahn-Hilliard diffuse interface model. To test our results, we compare the results of the pseudo-transient model with 3-D results obtained with the implicit numerical algorithm LAMEM (bitbucket.org/bkaus/lamem) and with 2-D results obtained with a Lagrangian finite element model. We have developed two versions of the algorithm, one in Matlab (mathworks.com) and one in Julia language (julialang.org). The aim of the Julia version is to eventually utilize a parallel GPU computing environment. Furthermore, we also present first results for cylindrical coordinates for the same plateau geometries. The aim of the model with cylindrical coordinates is to quantify the impact of the Earth’s curvature on the stress state.

How to cite: Macherel, E., Podladchikov, Y., Spitz, R., and Schmalholz, S. M.: Full 3-D pseudo-transient finite difference modelling of stress distribution around continental plateaus, EGU General Assembly 2020, Online, 4–8 May 2020, EGU2020-9005, https://doi.org/10.5194/egusphere-egu2020-9005, 2020.

Seismic anisotropy is a key property to understand the structure of the crust and mantle. In this contribution, we investigate the influence of shape (morphological) preferred orientation (SPO), crystallographic preferred orientation (CPO) and the spatial distribution of grains on seismic anisotropy in rocks (Bazargan et al., 2018). A numerical toolset has been developed with COMSOL to investigate these effects numerically, which has been benchmarked analytically and against other numerical models. Numerical samples modelled in 2D and 3D can determine anisotropy, by measurements along different sample axes, using different geometrical setups and mineral compositions. This numerical tool can include a variety of mineral arrangements and propagate P and S waves from different directions to calculate anisotropy. Current numerical results confirm directly the relations between the structural framework of the rocks (foliation, lineation) and velocity anisotropy, shear wave splitting and shear wave polarisation. This has been proven numerically with the effects of layering, which represents foliation and lineation in 2D. One of the aims of this work is to apply the fundamental results and effects of effective medium to improve our finite element method (FEM) toolbox to provide a numerical modelling tool for seismic data that have been collected in the field. Since the numerical and laboratory measurements are worked on together to verify the numerical results, to compare the models and explain the laboratory measurements have been conducted.

Here we also present laboratory measurements of directional dependence of elastic waves velocity and shear wave splitting to the internal rock structure. In the experimental part of this study, we illustrate the contribution of microstructural parameters (grain sizes, SPO and microcracks) to the elastic anisotropy of relatively similar quartzites and granites. An objective with the laboratory measurements is to investigate the effect of grain size and its possible influence on elastic wave speed and potential scattering effects due to wavelength effects. Granites are the one we use to investigate anisotropy related to SPO and CPO. Our experimental data consist of the measurements of elastic wave velocities (Vp, Vs₁ and Vs₂) at confining pressures up to 600 MPa (Bazargan et al., 2019).  numerical modelling together with laboratory measurements are used to obtain a better understanding of the role of microstructures in elastic wave propagation and its anisotropy

Bazargan, M. Almqvist, B. Hieronymus, Ch. Piazolo, S., Employing Finite Element Method using COMSOL multiphysics to predict seismic velocity and anisotropy: Application to lower crust and upper mantle rocks. EGU 2018.

Bazargan, M. Motra, H. B. Almqvist, B. G. Hieronymus, Ch. Piazolo, S., Elastic wave anisotropy in amphibolites and paragneisses from the Swedish Caledonides measured at high pressures (600 MPa) and temperatures (600 ^oC). EGU 2019.

How to cite: Bazargan, M., Bahadur Motra, H., Almqvist, B., Hieronymus, C., and Piazolo, S.: Numerical and experimental investigations of elastic wave anisotropy in monomineral and polymineral rocks, EGU General Assembly 2020, Online, 4–8 May 2020, EGU2020-10671, https://doi.org/10.5194/egusphere-egu2020-10671, 2020.

In recent years, the advent of ensemble-based methods in volcanology has greatly facilitated the use of numerical models within data assimilation frameworks that had previously been limited, either computationally or mathematically, to simpler analytical models. Because numerical models can simulate stress conditions throughout the model space, recent inversions based on assimilated volcanic deformation data are able to track not only the basic parameters of a magma reservoir, but also how those parameters affect the overall mechanical stability of the system. Although this approach has produced successful forecasts and hind-casts of volcanic eruptions, much work remains to be done in assessing its full capabilities and limitations. In particular, non-uniqueness in how source parameters are reflected in surface deformation can significantly impair the inversion’s ability to resolve the magma system’s true state and, by extension, the likelihood of eruption. While this problem is nearly intractable for deep reservoirs, for which changes in pressure and size are indistinguishable from deformation alone, preliminary synthetic tests at shallower systems have demonstrated a limited ability to resolve the main inflation mechanism. In this study, we investigate how the performance of an Ensemble Kalman Filter (EnKF) data assimilation framework varies under a wider range of experimental conditions than used in these initial investigations. In particular, we test how different mathematical implementations of the filter and how different levels of data availability affect the EnKF’s ability to distinguish inflation drivers and to accurately resolve reservoir parameters. To implement this experiment, two time series of synthetic GPS and InSAR data are generated, one in which deformation is driven by excess pressure and another in which it is driven by lateral expansion of the reservoir. For each filter implementation these datasets are down-sampled and given random noise prior to inversion, and after assimilation the resulting model is compared to the original synthetic conditions. We find that newer deterministic formulations of the EnKF are more accurate and consistent than the original stochastic implementation, although the improvement is relatively small. Moreover, some amount of parameter inflation is required to avoid model collapse, but more sophisticated adaptive inflation schemes do not produce better results than more basic formulations. Finally, we show that while increased data sampling does improve performance, this effect is subject to diminishing returns. In particular, data resolution near the center of inflation is more important than overall range of coverage. As new inversion techniques are developed or adapted from other fields, rigorous testing as demonstrated here will be a key step in being able to interpret future results and develop new forecasting frameworks for volcanic eruptions.

How to cite: Albright, J. and Gregg, P.: Optimizing Ensemble-Based Inversions for Non-unique Volcanic Systems, EGU General Assembly 2020, Online, 4–8 May 2020, EGU2020-10855, https://doi.org/10.5194/egusphere-egu2020-10855, 2020.

In the last decade, the development of numerical geodynamic tools helped the geoscience community to explore thermo-mechanical processes at play during plate tectonics. Yet, the high computational cost of thermodynamic calculations hampers our ability to quantify multi-phase systems in which the interplay between plate-tectonics and phase transformations leads to magmatism.  Here we use the 'igneous set' of HPx-eos (thermodynamic models for minerals and geological fluids that are based on the Holland & Powell dataset and 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 calculation is performed by Gibbs free energy minimization at prescribed pressure, temperature and bulk-rock composition and is achieved in two steps. First, we employ a levelling method (iterative change of base) to reduce the number of potential stable phases. Then, the composition and proportions of stable phases at equilibrium are determined using several constrained optimization methods. We explore the computational efficiency of linear programming (e.g., Simplex), Gradient-based (e.g., SLSPQ) and Hessian-based (e.g., Newton-Raphson) methods. The accuracy and performance of tested methods are compared, and applications to geodynamic modelling are discussed.

How to cite: Riel, N., Kaus, B., Berlie, N., Rummel, L., and Green, E.: Computational thermodynamics: towards an improved Gibbs minimization tool for geodynamic modelling, EGU General Assembly 2020, Online, 4–8 May 2020, EGU2020-11861, https://doi.org/10.5194/egusphere-egu2020-11861, 2020.

How to cite: Duretz, T., de Borst, R., and Räss, L.: Modelling lithosphere dynamics with robust rheological implementations: Towards 3D, EGU General Assembly 2020, Online, 4–8 May 2020, EGU2020-13742, https://doi.org/10.5194/egusphere-egu2020-13742, 2020.

The in-depth behaviour of magmatic systems is still poorly constrained due to their lack of accessibility and the difficulty of finding good analogue representations. However, progresses in the field of geological numerical modelling can allow to better understand and interpret those constraints. The MAGMA project aims on developing tools and software for studying a range of magmatic processes in the lithosphere. On the way to building a general framework able to model the behaviour of a fluid-solid chemically coupled magmatic system, we present here the current development of a mechanical staggered-grid finite difference 2D code using robust analytical linear and non-linear solvers via the PETSc infrastructure, able to run in parallel on highly performant computers. The mesh is assembled using the recently developed DMStag framework, which is part of PETSc. This code solves the Stokes equations for elasto-visco-plastic rheologies, including tensile plasticity essential to the development of dyke structures. Here, we will present the equations and implementations used, and show initial results and benchmarks.

How to cite: Berlie, N., Kaus, B., Popov, A., and Sanan, P.: Towards a staggered-grid finite difference code for modelling magmatic systems, EGU General Assembly 2020, Online, 4–8 May 2020, EGU2020-16485, https://doi.org/10.5194/egusphere-egu2020-16485, 2020.

All divergent plate boundaries are associated with magmatism, yet its role in their dynamics and deformation is not known. The RIFT-O-MAT project seeks to understand how magmatism promotes and shapes rifts in continental and oceanic lithosphere by using models that build upon the two-phase flow theory of magma/rock interaction. Numerical models of magma segregation from partially molten rocks are usually based on a system of equations for conservation of mass, momentum and energy. One key challenge of these problems is to compute a mass-conservative flow field that is suitable for advecting thermochemically active material that feeds back on the flow. This feedback tends to destabilise the coupled mechanics+thermochemical solver. 

Staggered grid finite-volume/difference methods are: mimetic (i.e., discrete differential operators mimic the properties of the continuous differential operators); conservative by construction; inf-sup stable and "light weight" (small stencil) thus they are well suited to address these problems. We present a new framework for finite difference staggered grids for solving partial differential equations (FD-PDE) that allows testable and extensible code for single-/two-phase flow magma dynamics. We build the framework using PETSc (Balay et al., 2019) and make use of the new features for staggered grids, such as DMStag. The aim is to separate the user input from the discretization of governing equations, allow for extensible development, and implement a robust framework for testing. Any customized applications can be created easily, without interfering with previous work or tests.

Here, we present benchmark and performance results with our new FD-PDE framework. In particular, we focus on preliminary results of a two-phase flow mid-ocean ridge (MOR) model with a free surface and extensional boundary conditions. We compare flow calculations with previous work on MORs that either employed two-phase flow dynamics with kinematic boundary conditions (i.e., corner flow, Spiegelman and McKenzie, 1987), or single-phase flow dynamics with free surface (i.e., Behn and Ito, 2008). In the latter case, the effect of magma is parameterised according to a priori expectations of its role. 

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

Spiegelman and McKenzie (1987), EPSL, 83 (1-4), 137-152.

Behn and Ito (2008), Geochem. Geophys. Geosyst., 9, Q08O10.

How to cite: Pusok, A. E., May, D. A., and Katz, R. F.: Magma dynamics using FD-PDE: a new, PETSc-based, finite-difference staggered-grid framework for solving partial differential equations, EGU General Assembly 2020, Online, 4–8 May 2020, EGU2020-18690, https://doi.org/10.5194/egusphere-egu2020-18690, 2020.

Brittle-plastic flows where the yield strength is a decreasing, non-linear function of plastic strain are thought to be commonplace in the Earth, and responsible for some of its most catastrophic events. Recent work [1] has highlighted again the computational benefit of an iterative Newton-Raphson scheme that contains a linearization of the plastic flow problem that is consistent with its time discretization. However, such a consistent linearization requires a nested set of iterations to converge on a yield strength if it is governed by a law that is non-linear in strain (or strain rate).

Eckert and co-authors [2] have shown that the construction of a consistent linearization can be avoided altogether, including these inner iterations, though at the considerable cost of including the full plastic strain tensor as an objective variable alongside the displacement vector. The resulting system is therefore larger, but as it can be expressed directly, posesses the quality that it may be linearized automatically, cheaply, and accurately by finite-differencing the non-linear residual with respect to the solution variables. Their algorithm naturally incorporates predictor and corrector polynomials that are second-order accurate in time, contrasting with traditional methods that are often derived using a Backward Euler time integrator. We present a modification to this algorithm that suppresses the cost of operating it significantly by replacing the symmetric second-order plastic strain tensor with a single effective plastic strain scalar objective variable, cutting the number of unknowns by 40% (2D) and 55% (3D) This makes it computationally more on par with existing schemes that employ a consistent tangent modulus.

We demonstrate this improved algorithm with test cases of non-linear strain softening laws relevant to Earth scientists, that include regularization by both Kelvin visco-plasticity [3] and non-local measures of effective plastic strain [4]. In addition, we analyse performance of this scheme with respect to existing algorithms.

References
[1] Duretz et al. (2018). “The benefits of using a consistent tangent operator for viscoelastoplastic computations in geodynamics.” Geochemistry, Geophysics, Geosystems, 19, 4904–4924.

[2] Eckert et al. (2004). “A BDF2 integration method with step size control for elasto-plasticity.” Computational Mechanics 34.5, 377–386.

[3] Duretz et al. (2019). “Finite Thickness of Shear Bands in Frictional Viscoplasticity and Implications for Lithosphere Dynamics.” Geochemistry, Geophysics, Geosystems, 20, 5598–5616.

[4] Engelen et al. (2003). “Nonlocal implicit gradient-enhanced elasto-plasticity for the modelling of softening behaviour.” International Journal of Plasticity
19.4, 403–433.

How to cite: Pranger, C., May, D., and Le Pourhiet, L.: An improved monolithic Newton-Raphson scheme for solving plastic flow with nonlinear flow laws, EGU General Assembly 2020, Online, 4–8 May 2020, EGU2020-19046, https://doi.org/10.5194/egusphere-egu2020-19046, 2020.

Kinematic plate tectonic reconstructions are central to our understanding of the dynamics of the solid Earth over hundreds of millions of years. These reconstructions are typically based on a combination of qualitative geological and paleontological observations supplemented with quantitative geophysical data, such as paleomagnetism. As a result, reconstructing plate and continental motion generally involves largely manual processes to integrate regional tectonic histories into a geometrically self-consistent global model. Further, with this methodology it is very difficult to quantify the uncertainties of the resulting time-dependent plate configurations and motions. To overcome these difficulties, ensemble-based data assimilation methods present a promising approach to paleogeographic reconstruction as they provide a formal statistical framework to assimilate time-dependent geoscientific data of variable nature and source within a dynamical model whilst providing quantitative estimates of uncertainties on the proposed trajectory.

To develop this concept, we apply a particle filter to reconstruct the time-dependent motion of continents on Earth. We start with a kinematic-stochastic model of continental drift: the motion of continents is governed by random rigid rotations, where velocity, rate of rotation (around polygon centroid) and rate of motion change are drawn randomly. Each of the probability density functions is chosen using geometrical and geodynamical arguments. This stochastic-kinematic model is then used to compute possible continent motion trajectories. We integrate this forward model into a data assimilation framework to incorporate paleomagnetic data, starting from present day and moving backward through time. From observing system simulation experiments using this technique, results suggest the number of ensemble members needed is of the order of several thousands to obtain accurate reconstructions of continent motions. Leveraging these experiments, we present the results of applying this method to a global paleomagnetic dataset spanning the last 100 Ma.

How to cite: Bocher, M. and Tetley, M.: Reconstructing global continental motions through time using data assimilation: a proof of concept., EGU General Assembly 2020, Online, 4–8 May 2020, EGU2020-20561, https://doi.org/10.5194/egusphere-egu2020-20561, 2020.

Porous flow is of major importance for reservoir-scale processes such as waste fluid sequestration or oil and gas 
exploration. The motion of a low-viscous fluid through a high-viscous matrix (rock) can be described by the coupled 
nonlinear hydro-mechanical equations. This two-phase flow may result in the initiation of porosity waves, triggering 
high-porosity vertical pipes or chimneys. Such fluid escape features may lead to localized and fast vertical flow 
pathways that may be problematic in the case of e.g. CO₂ sequestration. Determining the porosity in such environments 
is a major challenge. Seismic imaging methods can localize the high-porosity chimneys very well in the inverted wave speed 
field but the conversion to porosity is not straightforward. 

Here, we develop an inversion framework that allows us to invert for the porosity using fluid velocities as observables 
and investigate its behavior for simple examples. We introduce the adjoint framework for the two-phase flow equations, 
which allows for efficient calculation of the pointwise gradients of the flow solution with respect to the porosity. 
These gradients are then used in a gradient descent method to invert for the pointwise porosity. Technically, the forward 
and adjoint equations are solved using a parallel iterative finite-difference pseudo-transient approach, which executes 
optimally on latest manycore hardware accelerators such as GPUs. Numerical results show that an inversion for porosity is 
feasible and that the porosity is very locally sensitive to the fluid velocity.

How to cite: Holbach, L., Reuber, G., and Räss, L.: Gradient-based inversion for subsurface porosity using the adjoint two-phase flow equations: A pseudo-transient approach, EGU General Assembly 2020, Online, 4–8 May 2020, EGU2020-21578, https://doi.org/10.5194/egusphere-egu2020-21578, 2020.

Following large earthquakes, postseismic crustal deformations are often observed for more than years. They include the afterslip and the viscoelastic deformation of the crust and the upper mantle, activated by the coseismic stress change. The viscoelastic deformation gives the stress change on the neighboring faults, hence affects the seismic activity of the surrounding area, for a long period after the large earthquake. So, estimating the viscoelastic deformation after the large earthquakes is important.

In order to estimate the time evolution of the viscoelastic deformation after a large earthquake, we also need to know the viscoelastic structure around the area. Recently, the Ensemble Kalman filter method (EnKF), a sequential data assimilation method, starts to be used for the crustal deformation data to estimate the physical variables (van Dinther et al., 2019, Hirahara and Nishikiori, 2019). With data assimilation, we get a more provable estimation by combining the data and the time-forward model than only using the data. Hirahara and Nishikiori (2019) used synthetic data and showed that EnKF could effectively estimate the frictional parameters on the SSE (slow slip event) fault, addition to the slip velocity. In the present study, I applied EnKF to estimate the viscosity and the inelastic strain after a large earthquake, both the physical property and the variables.

First, I constructed the forward model simulating the evolution of the viscoelastic deformation, following the equivalent body force method (Barbot and Fialko, 2010; Barbot et al., 2017). This method is appropriate for applying EnKF, because the ground surface deformation rate is represented by the inelastic strain at the moment, and the history of the strain is not required. Then, we applied EnKF based on the forward model and executed some numerical experiments using a synthetic postseismic crustal deformation data.

In this presentation, I show the result of a simple setting. I assumed the medium to be two layers with a homogeneous viscoelastic region underlying an elastic region. The synthetic data is made by giving a slip on a fault at time t = 0 and simulating the time evolution of the ground surface deformation. For assimilation, I assumed that the slip on the fault and the stress distribution just after the large earthquake is known. Then we executed the assimilation every 30 days after the large earthquake. I found that I can get a good estimation of the viscosity after t > 150 days.

How to cite: Ohtani, M.: EnKF estimation of the viscoelastic deformation and the viscosity, EGU General Assembly 2020, Online, 4–8 May 2020, EGU2020-20884, https://doi.org/10.5194/egusphere-egu2020-20884, 2020.

State of the art thermo-mechanical models have become very efficient at testing scenarios of tectonic evolution but uncertainties on the rheologies and the complexity of the have so far limited the potential to quantitatively predict uplift and subsidence. Coupling thermo-mechanical models to landscape evolution models remains challenging and require careful validation and better integration of field data to prevent error in interpretation.

Low temperature thermochronology has been extensively used to quantitatively constrain the thermal histories of rocks. It can provide important information on tectonic uplift (or subsidence) by measuring the erosional (or burial) response and can also map the spatial and temporal pattern of geomorphic response of a landscape.

We use the temperature evolution of our coupled thermo-mechanical models with surface processes to predict Apatite fission track data (Ages and Track lengths distributions). The aim is to provide a direct means of comparison with actual empirical thermochronometric data which will allow different model scenarios and/or model parameter choices to be robustly tested.

We present a series of 3D coupled models (Underworld / Badlands) of Rifts and the associated Apatite Fission Track predicted by the thermal evolution of the rocks exhumed to the surface. We compare models predictions to existing thermochronological transects across passive margins.

We discuss the technical challenges in obtaining sufficiently high resolution temperature field and other associated challenges that need to be addressed to satisfactory apply our model to natural examples.

How to cite: Beucher, R., Moresi, L., Brown, R., and Mallard, C.: Thermomechanical models with surface processes: Low-Temperature Thermochronology predictions for model calibration, EGU General Assembly 2020, Online, 4–8 May 2020, EGU2020-1597, https://doi.org/10.5194/egusphere-egu2020-1597, 2020.

Energy and carbon-efficient exploitation, management, and remediation of subsurface aquifers, gas and oil resources, CO₂-disposal sites, and energy storage reservoirs all require high quality modelling and simulation. The heterogeneity and anisotropy of such subsurface formations has always been a challenge to modellers, with the best current technology not being able to deal with variations at scales of less than about 30-50 m. Most formations exhibit heterogeneities and anisotropy which result in variations of the petrophysical properties controlling fluid flow down to millimetre scale and below. These variations are apparent in well-logs and core material, but cannot be characterised in the inter-well volume which makes up the great majority of the formation.

This paper describes a new fractal approach to the modelling and simulation of heterogeneous and anisotropic aquifers and reservoirs. This approach includes data at all scales such that it can represent the heterogeneity of the reservoir correctly at each scale.

Advanced Fractal Reservoir Models (AFRMs) in 3D can be produced using our code. These AFRMs can be used to model fluid flow in formations generically to understand the effects of an imposed degree of heterogeneity and anisotropy, or can be conditioned to match the characteristics of real aquifers and reservoirs. This paper will show how 3D AFRMs can be created such that they represent critical petrophysical parameters, as well as how fractal 3D porosity and permeability maps, synthetic poro-perm cross-plots, water saturation maps and relative permeability curves can all be calculated. It will also show how quantitative controlled variation of heterogeneity and anisotropy of generic models affects fluid flow. We also show how AFRMs can be conditioned to represent real reservoirs, and how they provide a much better simulated fluid flow than the current best technology.

Results of generic modelling and simulation with AFRMs show how total hydrocarbon production, hydrocarbon production rate, water cut and the time to water breakthrough all depend strongly on heterogeneity, and also depend upon anisotropy. Modelling with different degrees and directions of anisotropy shows how critical hydrocarbon production data depends on the direction of the anisotropy, and how that changes over the lifetime of the reservoir.

Advanced fractal reservoir models are of greatest utility if they can be conditioned to represent individual reservoirs. We have developed a method for matching AFRMs to aquifer and reservoir data across a wide range of scales that exceeds the range of scales currently used in such modelling. We show a hydrocarbon production case study which compares the hydrocarbon production characteristics of such an approach to a conventional krigging approach. The comparison shows that modelling of hydrocarbon production is appreciably improved when AFRMs are used, especially if heterogeneity and anisotropy are high. In this study AFRMs in moderate to high heterogeneity reservoirs always provided results within 5% of the reference case, while the conventional approach resulted in massive systematic underestimations of production rate by over 70%.

How to cite: Lorinczi, P., Glover, P., Saud, A.-Z., Sinan, S., and Daniel, G.: Modelling and Simulation of Heterogeneous and Anisotropic Formations using Advanced Fractal Reservoir Models, EGU General Assembly 2020, Online, 4–8 May 2020, EGU2020-10072, https://doi.org/10.5194/egusphere-egu2020-10072, 2020.

Multiobservable thermochemical tomography (MTT) is a recent computational approach to obtain estimates of the physical state (e.g. temperature distribution, compositional structure and rock properties) for the upper mantle [1]. It allows to jointly invert multiple independent datasets (e.g. gravity, seismic, magnetotelluric) within a thermodynamically-constrained and fully probabilistic framework. Evaluation of the plausibility of different physical states of the mantle with Markov Chain Monte Carlo (MCMC) simulations requires the solution of complex forward problems (e.g. Stokes flow, Maxwell’s equations, etc.) millions of times, making MTT computationally demanding for large-scale inverse problems. Furthermore, the number of parameters in a global study can easily reach several millions, making it increasingly difficult to 1) locate the regions of high probability and 2) sample these regions appropriately.

In order to overcome these limitations, we have combined and implemented a number of techniques, such as reduced-order modelling and efficient parallelization of both the forward problems and the MCMC algorithms, which dramatically accelerate the solution of the forward problems. Our software, LitMod1D_4INV and LitMod3D_4INV, allow to compute a proposal in less than 1 second, even when solving multiple complex forward problems together. We develop a multi-level parallel MPI driver for a collection of advanced MCMC sampling strategies to locate and sample high-probability regions efficiently. The massive amounts of data generated by large-scale MTT inversions need to be managed efficiently. We output results to open-source freeware formats, such as HDF5, TileDB, designed for big data problems. We emphasize that our methods and approaches are not only useful for MTT, but for any demanding inverse problem.

In this contribution, we will present applications of our software to complex, large-scale MTT problems and discuss its benefits, limitations and future improvements.

[1] J.C. Afonso, N. Rawlinson, Y. Yang, D. L. Schutt, A. G. Jones, J. Fullea, W. L. Griffin, 3‐D multiobservable probabilistic inversion for the compositional and thermal structure of the lithosphere and upper mantle: III. Thermochemical tomography in the Western‐Central U.S., Journal of Geophysical Research, 121, doi:10.1002/2016JB013049, 2016

How to cite: Fomin, I. and Afonso, J.: Advances on multiobservable thermochemical tomography for the physical state of the upper mantle, EGU General Assembly 2020, Online, 4–8 May 2020, EGU2020-13423, https://doi.org/10.5194/egusphere-egu2020-13423, 2020.

Estimation of subglacial Geothermal Heat Flux (GHF) is of paramount importance to better understand the dynamics  of cryosphere and ice flow of the East Antarctica Ice Sheet (EAIS). Unfortunately, the GHF of East Antarctica is still poorly known and constrained, and direct measurements are still challenging. The EIAS is underlain by major subglacial mountain ranges and basins resulting from distinct geodynamic domains. These include Northern Victoria Land-Ross Sea, the Transantarctic Mountains, the Wilkes Subglacial Basin, the Gamburtsev Subglacial Mountains, the East Antarctic System and a major transpressional fault zone in between (e.g. Cianfarra & Maggi, 2017), which hosts clusters of subglacial lakes. The distribution of sedimentary basins and tectonic structures may affect the GHF in that it exhibits strong regional variations as testified by the presence of subglacial lakes at bedrock topographic elevation/depth with a range exceeding 1500 m, from deep subglacial basins to the flanking highlands. In the framework of the G-IDEA (Geo Ice Dynamics of East Antarctica) project, heat flow from the basement is quantified in key areas of East Antarctica between 60°E and 180°E, by an innovative application of the HCA (Hybrid Cellular Automata) method: the description of stationary conditions of the temperature field is used to replicate the observed distribution of wet vs dry ice-rock contacts in an ice-flowing environment. Evaluation of the geothermal flux is performed in key areas based on the numeric modeling of the ice-rock interaction, which can replicate the spatial distribution of wet contacts and subglacial lakes and is related to local dynamics of the ice sheet and its interaction with the atmosphere. The model takes into account the spatial distribution of the Curie temperature depth as derived from literature. The heat flux is estimated by modeling the stationary state of the ice-rock system with the HCA numerical method, and by its discretization into a large number of cells. Each cell is characterized by physical parameters such as density, enthalpy, thermal capacity and conductivity. By their interaction it is possible to compute their temperature evolution and to replicate the heat diffusion by conduction and convection (the ice movement) in the interfaces ice-rock and ice-atmosphere. The final resolution of the model is about 100 m. The presence of possible anomalous heath flow in the bedrock are identified by a stochastic approach that allow the estimation of the error in the computed heath flow values.

How to cite: Salvini, F., Cianfarra, P., Capponi, G., Crispini, L., Federico, L., and Rossi, C.: Geothermal Heat Flux in East Antarctica from HCA numerical modeling between 60-180°E Longitude, EGU General Assembly 2020, Online, 4–8 May 2020, EGU2020-19852, https://doi.org/10.5194/egusphere-egu2020-19852, 2020.

Advances in numerical modelling of geological processes are based upon, and driven by, diagnosing models. Such model diagnostics are often performed by hand, by eye, or else, by individually written routines that are neither tested or testable, nor reproducible.

Collecting geodynamic diagnostics, automating them in a robust manner to be applied to the multitude of different geodynamic models and codes, and providing them back to the community can foster additional progress within the modelling community.

In this presentation, I introduce the latest version of StagLab (Crameri 2018; www.fabiocrameri.ch/StagLab; currently version 5.0), which is a growing resource of geodynamic diagnostics, openly available, and easy to use. StagLab works seamlessly with StagYY (Tackley 2008) and can be made compatible with any other mantle convection code, if the respective developers start to provide machine-readable and documented output. Moreover, StagLab represents model data fairly to its users and to the readers of their papers. StagLab allows its users, whether professional or beginner, to produce state-of-the-art post-processing of geodynamic models, and publication-ready figures and movies, in a blink of an eye; all fully tested, testable and reproducible.

Crameri (2018), Geodynamic diagnostics, scientific visualisation and StagLab 3.0, Geosci. Model Dev., http://dx.doi.org/10.5194/gmd-11-2541-2018

Tackley (2008), Modelling compressible mantle convection with large viscosity contrasts in a three-dimensional spherical shell using the Yin-Yang grid, PEPI, http://dx.doi.org/10.1016/j.pepi.2008.08.005.

How to cite: Crameri, F.: Towards ready-to-use open source automated geodynamic diagnostics and fair representation of numerical models, EGU General Assembly 2020, Online, 4–8 May 2020, EGU2020-18384, https://doi.org/10.5194/egusphere-egu2020-18384, 2020.