ABSTRACT. We will provide a short overview of current applications of FEniCS in geodynamics and in particular in the study of the deformation of the Earth’s mantle and lithosphere. FEniCS provides a highly adaptable and accurate tool to study the long term thermal and chemical evolution of the Earth’s mantle, the structure and dynamics of subduction zones, and the migration of melts and fluids through a deforming solid matrix. We will focus here specifically on the first two topics. The latter topic is discussed in more detail in other presentations.

ABSTRACT. Subduction zones are a challenge to model numerically. As it descends into the Earth from the surface, the cold, subducting plate drives deformation in the hot mantle that now surrounds it. Simultaneously it releases aqueous fluids, which migrate towards the surface. The fluid migration happens on time-scales much faster than the solid deformation but may still influence it. Ultimately the fluids cause melting of the mantle and, once these melts reach the surface, the arc volcanism that we observe at subduction zones around the planet.

Over long periods of time the solid deformation caused by the sinking slab can be modeled as a viscously deforming single, solid phase using coupled energy and Stokes equations with an appropriately chosen non-linear rheological model for viscosity. Modeling the fluid transport requires a multi-phase description that incorporates chemical reactions between the phases, including de- and re-hydration of the solid, melting and freezing, and fluid and solid compositional changes. This is most commonly achieved using a porous media description with Darcy's equation coupled to the underlying Stokes equation for the solid. Even ignoring the fluid transport with a simple single-phase model of the thermal structure is computationally expensive, requiring large amounts of resolution, but fully-coupled multi-phase models are additionally highly non-linear, span a large range of time-scales, and require robust solution strategies to even converge.

We will present our progress modeling this system using FEniCS through our framework TerraFERMA and discuss some of the ongoing challenges and our strategies for resolving them.

ABSTRACT. Magmatism at mid-ocean ridges generates new oceanic crust and accounts for 60% of global volcanism. The oceanic crust is emplaced in a narrow neo-volcanic region, whereas basaltic melt is generated in a wide region beneath mid-ocean ridges as suggested by a few geophysical surveys. The combined observations suggest that melt generated in a wide region at depths has to focus to a small region at the surface. We present results from a suite of two-phase models applied to the mid-ocean ridges, varying half spreading rate and intrinsic mantle permeability using a new openly available, Melt in the Mantle beneath Mid-ocean ridges (M3LT-one) models, with a goal of understanding melt focusing beneath mid-ocean ridges. M3LT-one is built using the Transparent Finite Element Rapid Model Assembler, and solves the melt migration equations that were derived in dependently by several workers. TerraFERMA leverages open source libraries, FEniCs, PETSc and SPuD to provide a common interface for building custom finite element method.

We predict ocean crust thicknesses versus spreading rates from the M3LT-one model, which fit well with fundamental observations from geophysical surveys. Three distinct melt focusing mechanisms are recognized in these models: 1) Melting pressure focusing, 2) Decompaction layers and 3) Ridge suction, of which the first two play dominant roles in focusing melt. The manifestation of these mechanisms depends largely on the choices of rheological model. The models show that regardless of spreading rates, the amount of melt and melt transport patterns are sensitive to changes in intrinsic permeability, K0: 1) Intrinsic permeability affects melt transport efficiency and therefore the amount of porosity or melt fraction, i.e. larger permeability leads to faster melt velocity, therefore lowering porosity. 2) The change in porosity leads to an increase in bulk viscosity, which depends on eta/phi. At quasi-steady state, this affects the compaction pressure and its gradient on the axial melting region, which pulls melt towards the ridge axis. In particular, higher permeability models show a greater melting pressure focusing effectsince porosities are smaller. The lack of the decompaction layers in the geophysical observations hint at the possibility that melting pressure focusing could be more significant than that presented here, which could provide constraints to mantle rheology and permeability.

ABSTRACT. Mantle peridotite reacts rapidly with the atmosphere and surface waters in regions of the Earth where it is tectonically exposed at the surface. Such reactions are thought to play an important role in a wide range of geophysical and geochemical processes (e.g., in controlling the volume of water that is carried into the Earth at subduction zones), as well as a potential sink for carbon in engineered processes for reversing anthropogenic climate change. Extensive fracture networks are frequently observed in peridotites that have undergone hydration (forming serpentine) and carbonation (forming listvenite). It is thought that solid volume changes associated with these reactions can generate stresses significant enough to fracture the rock. This "reaction-driven cracking" has been hypothesized to create a positive feedback where new fluid pathways are created and fresh reactive surfaces are exposed, offering a potential mechanism to explain the high degrees of serpentinization/carbonation observed in outcrops. To describe such processes we propose a model of small-strain poro-elasticity that includes fluid flow and reaction, which has been extended to simulate brittle fracture using the phase-field method. We will discuss the FEniCS implementation of this model, numerical issues associated with the phase-field method, and some scientific implications of the results.

ABSTRACT. Subduction zone models (Wilson et al., this meeting) and "Reactive Cracking" models (Evans et al, this meeting) are prime examples of a general class of problems that can be describe as reactive fluid flow in strongly deformable, permeable media. A range of models built on FEniCS and PETSc (many through our framework TerraFERMA) have made progress in understanding the coupled fluid-solid mechanics of these systems. However to properly describe melting, hydration and carbonation reactions in open systems, along with accurate computation of material properties such as density, latent and sensible heat etc., requires integrating internally consistent thermodynamic models into the dynamic models. This talk will report on progress of the ENKI project, a NSF-SI^2 project to develop software for flexible construction and calibration of custom thermodynamic models for a range of Earth materials, as well as their integration into general multi-phase flow simulations.

Key components of the ENKI framework include python based tools for describing and recording models of the Gibbs free energy (or other Thermodynamic potentials) of a range of liquid and mineral phases. These models are described using SymPy which is used for automatic differentiation and code-generation of optimized C/C++ libraries (together with python bindings) for consistent thermodynamic properties of phases as a function of T,P and composition. These phase libraries are designed to be used in a range of applications including Bayesian calibration software (Wolf), optimization code for calculating thermodynamic equilibrium and phase diagrams (Ghiorso) , as well as kinetic reaction modeling for disequilibrium transport. Here we will primarily describe the kinetics modeling framework we have developed. This software allows for flexible description of reaction kinetics between phases and again uses SymPy to generate efficient C++ libraries and python bindings for direct inclusion in dynamic modeling software. We will demonstrate its use and behaviour in TerraFERMA/FEniCS models for a range of Earth Science applications and discuss steps required for integration into large scale subduction zone and carbon sequestration models.

ABSTRACT. Many real world applications occur in heterogeneous media, for example, subsurface flow and processes in composite materials. For accurate numerical solution of such problems, we should use a very fine grid that resolve all small scale heterogeneities which leads to the computationally expensive algorithms. To reduce size of the discrete systems, we construct a coarse grid approximation and calculate upscaled parameters of the model. For problems in stochastic media, we should recalculate upscaled parameters for each realization of the heterogeneous properties.

In this work, we present a machine learning algorithm for fast prediction of the upscaled properties of the model. The proposed method based on the construction of the convolutional neural network and fast training on the GPU. We train neural networks on the set of the selected realizations of the local microscale stochastic fields to learn a map between stochastic fields and upscaled properties. Numerical results are presented for two and three-dimensional model problems and show that proposed method provide fast and accurate upscaled property predictions.

ABSTRACT. In this work, a single fluid flow poroelastic model for simulating triaxial cell tests implemented in FeniCS project is presented. The mathematical model follows the Biot theory of poroelasticity considering material transverse isotropy, as well as the compressibility of the pore fluid and the solid constituents of the rock skeleton and includes geomechanics and flow phenomena such as displacements and pore fluid pressure [1]. For the numerical solution of the equation system, a finite element method (FEM) in space and a backward Euler finite difference method in time are applied, resulting in a fully implicit scheme. Its computational implementation was carried out in the programing language Python using the open-source computing platform FeniCS project [2]. From the methodological point of view, each stage of model development (conceptual, mathematical, numerical and computational) is described. The resulting model is applied to the classical Mandel problem [3] using data from a case study given in [4]. The numerical solutions are compared with the analytical ones for validation purposes of the FEniCS code.

References

[1] Biot, M. A. (1955). Theory of elasticity and consolidation for a porous anisotropic solid. J.Appl. Phys., 26:182 – 185.

[2] Logg, A., Mardal, K.-A., and Wells, G. N. (2012). Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book. Springer-Verlag, Berlin.

[3] Abousleiman, Y., Cheng, A. H.-D., Cui, L., Detournay, E., and Roegiers, J.-C. (1996). Mandel’s problem revisited. Geotechnique, 46(2):187–195.

[4] Sangnimnuan, A., Li, J., and Wu, K. (2018). Development of efficiently coupled fluid- flow/geomechanics model to predict stress evolution in unconventional reservoirs with complex-fracture geometry. SPE Journal, pages 640–660.

ABSTRACT. Plastic deformation is the crucial features of material science which is occurred by the motion of dislocations. During this process, dislocations may group together into patterns which are spontaneous emergence of low and high dislocation density. In this presentation, within continuum dislocation dynamics theory we investigate how to deal with dislocations and their motion to obtain such patterns.