DOUAR: A new three-dimensional creeping flow numerical model for the solution of geological problems

Abstract

We present a new finite element code for the solution of the Stokes and energy (or heat transport) equations that has been purposely designed to address crustal-scale to mantle-scale flow problems in three dimensions. Although it is based on an Eulerian description of deformation and flow, the code, which we named DOUAR (‘Earth’ in Breton language), has the ability to track interfaces and, in particular, the free surface, by using a dual representation based on a set of particles placed on the interface and the computation of a level set function on the nodes of the finite element grid, thus ensuring accuracy and efficiency. The code also makes use of a new method to compute the dynamic Delaunay triangulation connecting the particles based on non-Euclidian, curvilinear measure of distance, ensuring that the density of particles remains uniform and/or dynamically adapted to the curvature of the interface. The finite element discretization is based on a non-uniform, yet regular octree division of space within a unit cube that allows efficient adaptation of the finite element discretization, i.e. in regions of strong velocity gradient or high interface curvature. The finite elements are cubes (the leaves of the octree) in which a $q_1$–$p_0$ interpolation scheme is used. Nodal incompatibilities across faces separating elements of differing size are dealt with by introducing linear constraints among nodal degrees of freedom. Discontinuities in material properties across the interfaces are accommodated by the use of a novel method (which we called divFEM) to integrate the finite element equations in which the elemental volume is divided by a local octree to an appropriate depth (resolution). A variety of rheologies have been implemented including linear, non-linear and thermally activated creep and brittle (or plastic) frictional deformation. A simple smoothing operator has been defined to avoid checkerboard oscillations in pressure that tend to develop when using a highly irregular octree discretization and the tri-linear (or $q_1$–$p_0$) finite element. A three-dimensional cloud of particles is used to track material properties that depend on the integrated history of deformation (the integrated strain, for example); its density is variable and dynamically adapted to the computed flow. The large system of algebraic equations that results from the finite element discretization and linearization of the basic partial differential equations is solved using a multi-frontal massively parallel direct solver that can efficiently factorize poorly conditioned systems resulting from the highly non-linear rheology and the presence of the free surface. The code is almost entirely parallelized. We present example results including the onset of a Rayleigh–Taylor instability, the indentation of a rigid-plastic material and the formation of a fold beneath a free eroding surface, that demonstrate the accuracy, efficiency and appropriateness of the new code to solve complex geodynamical problems in three dimensions.

Introduction

In recent years, modelling the Earth’s crust and upper mantle deformation has led to increased insight concerning the way the Earth responds to both tectonic and erosional forciong and indirectly to the climate Willett et al., 1993, Batt and Braun, 1997, Beaumont et al., 2001. Apart from a few exceptions Braun, 1993, Braun, 1994, Braun and Beaumont, 1995, such modelling has been limited to two-dimensional analysis, mostly for computational efficiency reasons. Convection or mantle-scale flow calculations have been more easily performed in three dimensions, in part due to the relatively simple geometry of the problem and the lack of interactions with a free or eroding upper surface (Houseman, 1988, Albers, 2000, Tackley, 1998, among many others).

The need for a three-dimensional model capable of taking into account the large stresses arising from surface topography gradient is however growing (Braun, 2006). Key questions regarding the potential couplings and feedbacks between tectonics, erosion and climate can only be properly addressed using a plan-view surface processes model and, thus a full three-dimensional representation of deformation in the underlying crust (Stolar et al., 2005).

Three dimensional calculations are inherently computationally costly. Using a uniform spatial discretization, most three-dimensional convection models are limited in their spatial resolution to meshes of the order of $100^3$ elements or nodes (Tackley, 1998, for example). Owing to the non-linear and localizing nature of lithospheric rheologies, deformation gradients are much greater in the lithosphere than in the convecting parts of the Earth’s mantle and a finer spatial discretization is required to capture them. This is the reason why complex meshing algorithms have been developed and are commonly used for the solution of lithospheric-scale deformation or flow problems in two dimensions (Braun and Sambridge, 1994). These are, however, difficult to generalize to three dimensions. Furthermore, the evolving nature of the deformation or velocity field (in part due to the formation of shear zones or faults) requires the use of adaptative meshing techniques in which the numerical spatial discretization evolves with the flow.

Two contrasting methodologies have commonly been used to solve lithospheric-scale deformation problems. Explicit time-stepping methods are based on the dynamical force balance equation ($\text{F}=m\gamma$) in which a pseudo-mass (m) has been introduced to damp numerical oscillations. These methods require relatively few operations per time steps but a large number of time steps. There have been several implementations of these implicit algorithms to solve problems of lithospheric deformation in two dimensions Poliakov et al., 1993, Hassani et al., 1997. Remarkably, none of these have been so far ported to three-dimensions. Implicit methods in which the equations of static equilibrium have been linearized to form a large system of algebraic equations require less time steps but become computationally expensive in three dimensions. Multigrid iterative methods are commonly used to solve these large systems of equations but their convergence is poor when dealing with highly non-linear problems or those involving a free surface (Moresi and Solomatov, 1998).

Here we present a newly developed finite element code to solve the three dimensional Navier–Stokes equations that we purposely developed to address $Re=0$, high $Ra$ and infinite $Pr$ flows characterized by a free and/or eroding surface. The new model, that we called DOUAR, is in principle capable of tracking any interface, such as the Moho or a stratigraphic marker, deforming with the flow. It is based on a multi-scale octree-based discretization method and uses a fast, yet accurate direct solver for the solution of the large system of algebraic equations resulting from the implicit time-stepping, finite element discretization of the static force balance equations. In this paper, we present in detail the various components of this new code that are based on existing and novel algorithms. We also present the results of selected computations that demonstrate the usefulness and accuracy of the methodology.