Comparison of Eulerian and Lagrangian numerical techniques for the Stokes equations in the presence of strongly varying viscosity
Introduction
Numerical modeling of geodynamic processes on geological timescales requires solving of Stokes-like equations regardless of whether the rheology is viscous, non-Newtonian, visco-elastic or visco-elasto-plastic (e.g. Gerya and Yuen, 2003, Gerya and Yuen, 2007, Braun and Sambridge, 1994, Schmeling and Jacoby, 1981, Moresi et al., 2002, Moresi et al., 2003, Moresi et al., 2007, Zhong et al., 2007, Schmalholz et al., 2001, Vasilyev et al., 2001, Fullsack, 1995, Kaus et al., 2004). A major complication in numerically solving such problems is that the effective material properties of rocks vary many orders of magnitudes due to either changes in temperatures or differences in composition. For example, the effective viscosity of olivine has an Arrhenius-type temperature dependence of the form . Representative values are kJ/mol, J/(K mol). This results in an effective viscosity contrast between near surface conditions ( K) and the upper mantle ( K) of Pa s. Similarly, to solve the problem of propagation of basaltic dikes ( Pa s; Shaw, 1969) through the visco-elastic upper crust ( Pa s; Ranalli, 1995), the numerical method should be capable of handling viscosity contrasts of over length scales of tens of meters. Similar difficulties, although slightly less dramatic, occur in visco-plastic problems, in which the effective viscosity inside a shear zone typically is several orders of magnitude smaller than outside the shear zone (with a width of a few numerical grid points).
The Stokes problem, in the case of constant viscosity, is one of the classical problems in numerical mathematics, as it is the prototype equation for saddle-point problems. Therefore, numerous methods have been developed to solve the equations efficiently (see e.g. Elman et al., 2005, for an overview). Unfortunately these methods are insufficiently studied for cases where viscosity is strongly varying as the Stokes equations become very ill-conditioned. This makes it difficult to directly apply those methods and so far almost no theoretical work seems to have been performed for geodynamically relevant cases.
Moresi and Solomatov (1995) developed a multigrid method for 2D convection problems capable of handling viscosity variations of up to 14 orders of magnitude. Although this is an impressive result, the reason why a solution could be found with iterative methods could be related to the relatively smooth variations in viscosity (which was also found by Auth and Harder, 1999, Kameyama et al., 2005, Choblet, 2005). Sharp gradients in effective viscosities cause a deterioration in the performance of numerical methods, and may also lead to non-convergence of solutions (e.g. Tackley, 2000, Choblet, 2005, Albers, 2000, Trompert and Hansen, 1996). Indeed, it was shown by Moresi et al. (1996), in which direct and multigrid methods were compared with analytical solutions, that the accuracy of the solution depends on the local (elementwise) variation of viscosity. If the jump in material properties occurs exactly at the element boundary, the errors in velocity and pressure are relatively independent of the total viscosity contrast. If the jump in material properties occurs inside an element, however, errors increased by more than two orders of magnitude. Finally, they demonstrated that the errors in velocity are one to two orders of magnitudes smaller than the errors in pressure.
This was also found in Schmid (2002), in which the dynamic pressure solution around circular and elliptic inclusions was compared with a numerical FE solution. Schmid found that it was necessary to use an unstructured FEM, together with higher order () shape functions (instead of the linear shape functions employed by, e.g. Moresi et al., 1996) to obtain accurate pressure solutions.
If one is interested in accurate velocity solutions only, most methods seem to be reliable provided no strong viscosity contrasts occur within an element. If one, however, is also interested in accurate dynamic pressures one has to consider possible effects of the definition of material properties in numerical grids.
Motivation for the development of numerical techniques that accurately resolve the pressure field during deformation comes from a variety of sources:
- •
It has become more common to track the pressure–temperature evolution of rocks in numerical codes of lithospheric deformation, for the purpose of comparison with observed data (e.g. Toussaint et al., 2004, Jamieson et al., 2006, Burg and Gerya, 2006). Since tectonic stress can be significant under certain circumstances (Petrini and Podladchikov, 2000), it is crucial to use dynamic pressure instead of lithostatic load for reconstructing P–T paths.
- •
Treatment of plasticity requires knowledge of the absolute values of stress. Using dynamic pressure for this improves the sharpness of shear bands (Gerya and Yuen, 2007). Moreover, one of the interesting features of Mohr–Coulomb plasticity is that it lowers the pressure inside the shear-band, which makes this rheology different than that of ductile shear zones Vermeer, 1990, Mancktelow, 2006.
- •
Fluid flows down pressure gradients. Lowering the solid pressure locally thus focuses fluids into brittle shear zones, a mechanism that is often invoked to explain strain weakening (Huismans and Beaumont, 2002). Self-consistent modeling of such strain-weakening processes requires coupling of solid deformation with Darcy-type flow, and hence requires accurate dynamic pressure solutions (Morency et al., 2007).
- •
Melt migration through compacting, two-phase flow materials requires solving equations for the fluid and the solid matrix (e.g. McKenzie, 1984, Scott and Stevenson, 1986). With few notable exceptions Katz et al., 2006, Scott, 1988, Spiegelman, 2003 most formulations restrict their analysis to cases in which the solid matrix is rigid. Whereas this greatly simplifies the analysis, a full description of the problem seems necessary for certain problems such as shear-localization in partially molten rocks (Katz et al., 2006), and requires accurate stresses in the solid matrix.
These applications demonstrate that it is important to understand the accuracy of numerical methods for Stokes flow in the presence of large variations in material properties. As it is described above, the studies of Moresi et al. (1996) and Schmid (2002) show that the FEM would be an appropriate tool to solve such problems. Its disadvantage is that it requires exact knowledge of locations of material interfaces; there are many situations in geodynamic modeling in which this is not the case, or in which it would be prohibitively expensive to dynamically generate such a mesh. Examples are cases in which melt spontaneously localizes into small zones (e.g. Connolly and Podladchikov, 2007), models with shear-localization (e.g. Kaus and Podladchikov, 2006, Buiter et al., 2006) and geodynamic models of multi-phase lithosperic-scale deformation (e.g. Gorczyk et al., 2007, Fullsack, 1995). The most practical methods for such problems, in particular in 3D, seem to be Eulerian or arbitrary Lagrangian–Eulerian methods with quasi-structured grids as they are routinely employed in the mantle convection community.
Among the described methods there are other methods, which might be potentially powerful in tracking free surfaces such as an interface between two different materials. Some examples are Level-set, boundary integral or volume-of-fluid methods Mühlhaus et al., 2007, Sussman et al., 1994. Each of this methods has its own advantages and disadvantages. The volume-of-fluid method has difficulties to handle merging and folding interfaces. It would require reordering of the interface points and would result in complex programming. Boundary integral methods require the integral description of the equations and might become difficult to implement if viscosity is temperature dependent or if phase transitions occur. While all these interface tracking methods work very well for simple problems with only few different materials and interfaces, its implementation for multi-phase problems is insufficiently studied.
Thus the objective of this study is to obtain a better insight in the accuracies of the dynamic pressure and velocity solution for various 2D FDM and FEM of the Stokes equations. It has been demonstrated that FD and finite volume methods are very efficient in solving multi-component geodynamic problems (e.g. Gerya and Yuen, 2007, Schott et al., 2000, Weinberg and Schmeling, 1992, Katz et al., 2006). The main aims of the present work are
- •
To demonstrate that, with relatively simple suggestions, the accuracy of the pressure solution of the Stokes solver can be improved. Many authors use FDM to model geodynamical problems. Models using FDM require the application of interpolation schemes; if the interpolation is carefully chosen the accuracy of solution can be improved by almost one order of magnitude.
- •
To compare the accuracy of FEM with those of FDM. Of interest is also how accurate FEM soltions are if phase boundaries are not exactly followed by the mesh (i.e. Eulerian FEM).
We employ three different FD formulations: (1) a staggered grid velocity–pressure FDM, (2) a stream function FDM, and (3) a rotated staggered grid velocity–pressure FDM, which was recently demonstrated to be effective for simulating wave propagation in heterogeneous rocks and an unstructured FEM using different element setups. Numerical solutions are compared with analytical solutions for the flow around circular inclusions of different viscosity and with analytical solutions of density-driven flow. The staggered grid and stream function FDM require viscosities to be defined both at centers and at nodal points of a cell, while the rotated staggered grid FDM only requires viscosity defined at center points. We demonstrate that the manner in which viscosities are defined at these locations is of extreme importance for the accuracy of the overall solution. In case of the FEM the accuracy is affected by the type of mesh as well as by the averaging method of viscosity.
In the following, we begin by describing the governing equations. We than illustrate the problem of viscosity interpolation with a 1D problem for which an analytical solution can be found. We continue with the description of the numerical methods, definition of benchmark studies and a discussion of our results and conclusions.
Section snippets
Governing equations
The governing continuum mechanics equations for slowly moving Newtonian viscous, incompressible flow are the Stokes equations (e.g. Turcotte and Schubert, 1982, Schubert et al., 2001):where denotes velocity, deviatoric stress, total stress, P pressure, Newtonian viscosity, density, gravitational acceleration, the Kronecker delta, and the Einstein summation convention is assumed.
1D analytical solution and its numerical application
To obtain insight in the accuracy of numerical solutions in the presence of strong viscosity variations, we derive an analytical solution for a quasi 1D case, which we employ to estimate the accuracy of numerical solutions. We consider a thin bar of length L that is composed of two regions with different constant viscosity separated by an interface. The system is subjected to stress boundary conditions (Fig. 1).
2D governing equations
The Stokes equations (Eqs. (1), (2), (3), (4), (5)) in 2D are
Eq. (36) represents conservation of mass, Eqs. (37), (38) represent conservation of momentum in x- and in z-direction.
In presence of highly varying viscosity, the Stokes solver is critical for accurate results. Therefore, accuracy tests for pressure and velocity solutions are performed for three different FDM and an FEM using different arrangements of
Viscosity interpolation methods for finite difference methods
Discretization of the Stokes equations for the staggered grid and stream function FDM requires viscosity to be defined at both nodal and center points of a control volume (Fig. 6(a) and (b)). The 1D analytical solution (Section 3) showed that an effective viscosity should be calculated if an interface intersects a cell. Markers can be used to evaluate a more exact location of an interface intersecting a cell (marker-in-cell method employed in Gerya and Yuen, 2003, Weinberg and Schmeling, 1992,
Model setups
The accuracy of the three different FD formulations was compared to analytical solutions in two different cases: (1) a circular inclusion test and (2) a Rayleigh–Taylor instability test. The FEM results were only compared with the circular inclusion test.
2D pressure distribution
In Schmid and Podladchikov (2003) a 2D analytical solution is derived for the setup given in Section 6.1. Pressure is defined throughout the whole domain and is directly used to evaluate the accuracy of the numerical schemes. For the previously described pure shear boundary conditions, the analytical solution has zero pressure inside the strong clast and highest and lowest pressures directly at the clast–matrix boundary (Fig. 10(a)).
Fig. 10(b) shows the numerically computed pressure
Discussion
The analysis presented here demonstrates how material properties are defined in different FD formulations in presence of strongly varying material properties, in particular viscosity. Quantitatively it was shown that the accuracy of the pressure (circular inclusion test) or velocity solution (circular inclusion and Rayleigh–Taylor instability test) is critically dependent on the averaging scheme. We also demonstrated, that certain arrangements of elements in FEM show large differences in
Conclusions
In this paper we have analyzed the effects of viscosity averaging and numerical methodology on the accuracy of the Stokes equations in the presence of viscosity variations.
In 1D analytical and numerical tests, we demonstrated that the optimal viscosity in a cell, where an interface is separating two domains of different viscosities, is a harmonic average weighted according to the fractions of the two viscosities. The harmonically marker ratio average could represent the solution very well,
Acknowledgements
This work was supported by the Swiss National Science Foundation grant 200021-107889. Helpful comments on the manuscript of James Connolly and Stefan Schmalholz are gratefully acknowledged. We thank Taras Gerya and Erik Saenger for stimulating discussions. We further thank the reviewers R. Katz and M. Krotkiewski and the editor D. Schmid for insightful and helpful reviews.
References (66)
A local mesh refinement multigrid method for 3D convection problems with strongly variable viscosity
Journal of Compuational Physics
(2000)- et al.
Dynamic Lagrangian remeshing (DLR): a new algorithm for solving large strain deformation problems and its application to fault-propagation folding
Earth and Planetary Science Letters
(1994) Modelling thermal convection with large viscosity gradients in one block of the “cubed sphere”
Journal of Computational Physics
(2005)- et al.
Characteristics-based marker-in-cell method with conservative finite-differences schemes for modeling geological flows with strongly variable transport properties
Physics of the Earth and Planetary Interiors
(2003) - et al.
Multigrid iterative algorithm using pseudo-compressibility for three-dimensional mantle convection with strongly variable viscosity
Journal of Computational Physics
(2005) - et al.
Numerical simulation of geodynamic processes with the portable extensible toolkit for scientific Computation
Physics of th Earth and Planetary Interiors
(2007) - et al.
A Lagrangian integration point finite element method for large deformation modeling of viscoelastic geomaterials
Journal of Computational Physics
(2003) - et al.
The accuracy of finite element solutions of Stokes’s flow with strongly varying viscosity
Physics of the Earth and Planetary Interiors
(1996) - et al.
Modeling the propagation of elastic waves using a modified finite-difference grid
Wave Motion
(2000) - et al.
The diversity of tectonics from fluid-dynamical modeling of the lithosphere–mantle system
Tectonophysics
(2000)
A level set approach for computing solutions to incompressible two-phase flow
Journal of Computational Physics
Polydiapirs: multiwavelength gravity structures
Journal of Structural Geology
Numerical methods in mantle convection
Treatise in Geophysics
Multigrid solution of convection problems with strongly variable viscosity
Geophysical Journal International
Arbitrary discontinuitites in finite elements
International Journal for Numerical Methods in Engineering
Theory of folding of stratified viscoelastic media and its implications in tectonics and orogenesis
Geological Society of America Bulletin
Theory of gravity instability with variable overburden and compaction
Geophysics
The role of viscous heating in barrovian metamorphism of collisional orogens: thermomechanical models and application to the lepontine dome in the central alps
Journal of Metamorphic Geology
Dome structures in collision orogens mechanical investigation of the gravity/compression interplay
Decompaction weakening and channeling instability in ductile porous media: implications for asthenospheric melt segregation
Journal of Geophysical Research
Finite element methods and Navier–Stokes equations
Mathematics and its Applications
Finite elements and fast iterative solvers with applications to incompressible fluid dynamics
Numerical Mathematics and Scientific Computation
Computational Techniqus for Fluid Dynamics, vol. 2
An arbitrary Lagrangian–Eulerian formulation for creeping flows and its application in tectonic models
Geophysical Journal International
Robust characteristics method for modelling multiphase visco-elasto-plastic thermo-mechanical problems
Physics of the Earth and Planetary Interiors
Growth and mixing dynamics of mantle wedge plumes
Geology
Asymmetric lithospheric extension: the role of frictional plastic strain softening inferred from numerical experiments
Geology
On local tracking algorithms for the simulation of three-dimensional
Comput. Mech.
Provenance of the greater himalayan sequence and associated rocks: predictions of channel flow models
The dynamics of melt and shear localization in partially molten aggregates
Nature
Eulerian spectral/finite difference method for large deformation modelling of visco-elasto-plastic geomaterials
Bolletino Geofisico
Cited by (58)
Numerical Modeling of Subduction
2023, Dynamics of Plate Tectonics and Mantle ConvectionStress recovery for the particle-in-cell finite element method
2021, Physics of the Earth and Planetary InteriorsCitation Excerpt :As subdivision of the mesh by searching for material interfaces increases the computational complexity (Puckett et al., 2018), smoothing the discontinuities by averaging particle properties with the arithmetic or harmonic mean (Deubelbeiss and Kaus, 2008; Schmeling et al., 2008) or using a smooth basis function in the FEM (Bardenhagen and Kober, 2004) can be a simpler approach. Averaging of the viscosity in cells with multi-phase materials in the FD method was systematically studied by Deubelbeiss and Kaus (2008). However, they did not make full use of the possibilities offered by the FEM, as they use a single value for all quadrature points in one element.
Modelling the interplate domain in thermo-mechanical simulations of subduction: Critical effects of resolution and rheology, and consequences on wet mantle melting
2017, Physics of the Earth and Planetary InteriorsCitation Excerpt :On the contrary, a harmonic averaging simulates simple shear (at constant stress) parallel to the channel surface. While arithmetic averaging yields an artificial strengthening, the harmonic one reproduces well a thin channel low viscosity, and is more appropriate to simulate a subduction channel (Schmeling et al., 2008; Deubelbeiss and Kaus, 2008). As both kinds of shear are expected in a subducting slab, the averaging scheme should change depending on deformation pattern, which is not easy to handle numerically (Schmeling et al., 2008).
Direct numerical simulations of gas–solid–liquid interactions in dilute fluids
2017, International Journal of Multiphase FlowCitation Excerpt :We argue that this hybrid approach reduces computational cost without a significant loss in accuracy (Section 4.1). Most importantly, we avoid spurious pressure oscillations in the vicinity of the interface (Fig. 6), which are a common problem in previous approaches (Deubelbeiss and Kaus, 2008). Tracking the gas phase through a level-set function enables us to model strongly deforming interfaces that may exhibit topological changes like bubble coalescence and bursting.