Comparison of Eulerian and Lagrangian numerical techniques for the Stokes equations in the presence of strongly varying viscosity

https://doi.org/10.1016/j.pepi.2008.06.023Get rights and content

Abstract

Numerical modeling of geodynamic problems typically requires the solution of the Stokes equations for creeping, highly viscous flows. Since material properties such as effective viscosity of rocks can vary many orders of magnitudes over small spatial scales, the Stokes solver needs to be robust even in the case of highly variable viscosity. A number of different techniques (e.g. finite difference (FD), finite element (FE) and spectral methods) are presently in use. The purpose of this study is to evaluate the accuracies of several of these techniques. Specifically, these are staggered grid, stream function and rotated staggered grid finite difference method (FDM) and an unstructured finite element method (FEM) with various arrangements of elements. Results are compared with two different analytical solutions: (1) stress distribution inside and around strong or weak viscous inclusions subjected to pure-shear and (2) density-driven flow of a simple two-layer system. The comparison shows, in case of the three FD techniques, that the manner in which viscosity parameters are defined in the numerical grid plays an important role. The application of different viscosity interpolation methods yields differences in accuracy of up to one order of magnitude. In case of the FEM, results show that the arrangement of the elements in the region of material interfaces and the way in which material properties are defined also strongly affects the accuracy.

We derived a 1D analytical solution for a simple physical model where an interface separates two domains of different viscosities. If the interface is located between two nodal points the effective viscosity for this cell is a harmonic average of the two viscosities weighted according to the fraction of each material in this cell. Numerically, fractions were evaluated by using markers bearing the material property information.

The 2D problem differs from the 1D problem in that in 2D, two viscosities are necessary for a conservative FD formulation, one in the center and one at nodal points of a cell. Harmonic (in some cases geometric) averaging of the viscosity from markers to center points and a second harmonic averaging from center to nodal points is shown to give the most accurate FD results. In presence of density variations additionally a marker based arithmetic average for density should be applied. Results obtained by an unstructured FEM with elements following exactly the material boundaries show accuracies of one to two orders of magnitude better then results produced by Eulerian FD or FE methods. If viscosities are directly sampled at integration points, however, the FEM is less accurate than FD methods. Elementwise averaging of viscosities yields similar accuracies for both methods.

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 ηexp(Q/(nRT)). Representative values are Q400 kJ/mol, R=8.3145 J/(K mol). This results in an effective viscosity contrast between near surface conditions (T400 K) and the upper mantle (T1600 K) of η/η01039 Pa s. Similarly, to solve the problem of propagation of basaltic dikes (η100 Pa s; Shaw, 1969) through the visco-elastic upper crust (η1023 Pa s; Ranalli, 1995), the numerical method should be capable of handling viscosity contrasts of 1021 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 (Q2P1) 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 PT 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):vixi=0,σijxj=ρgi,σ˜ij=2μɛ˙ij,ɛ˙ij=12vixj+vjxi,σij=δijP+σ˜ij,where vi denotes velocity, σ˜ij deviatoric stress, σij total stress, P pressure, μ Newtonian viscosity, ρ density, gi gravitational acceleration, δij 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 arevxx+vzz=0,Px+2xμvxx+zμvxz+vzx=0,Pz+2zμvzz+xμvxz+vzx=ρg.

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)

  • M. Sussman et al.

    A level set approach for computing solutions to incompressible two-phase flow

    Journal of Computational Physics

    (1994)
  • R.F. Weinberg et al.

    Polydiapirs: multiwavelength gravity structures

    Journal of Structural Geology

    (1992)
  • S. Zhong et al.

    Numerical methods in mantle convection

    Treatise in Geophysics

    (2007)
  • C. Auth et al.

    Multigrid solution of convection problems with strongly variable viscosity

    Geophysical Journal International

    (1999)
  • T. Belytschko et al.

    Arbitrary discontinuitites in finite elements

    International Journal for Numerical Methods in Engineering

    (2001)
  • M. Biot

    Theory of folding of stratified viscoelastic media and its implications in tectonics and orogenesis

    Geological Society of America Bulletin

    (1961)
  • M. Biot et al.

    Theory of gravity instability with variable overburden and compaction

    Geophysics

    (1965)
  • Buiter, S.J.H., Babeyko, A., Ellis, S., Gerya, T., Kaus, B., Kellner, A., Schreurs, G., Yamada, Y., 2006. The numerical...
  • J.P. Burg et al.

    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

    (2006)
  • J.P. Burg et al.

    Dome structures in collision orogens mechanical investigation of the gravity/compression interplay

  • J.A.D. Connolly et al.

    Decompaction weakening and channeling instability in ductile porous media: implications for asthenospheric melt segregation

    Journal of Geophysical Research

    (2007)
  • C. Cuvelier et al.

    Finite element methods and Navier–Stokes equations

    Mathematics and its Applications

    (1986)
  • Dabrowski, M., Krotkiewski, M., Schmid, D.W., 2008. MILAMIN: MATLAB-based FEM solver for large problems. Geochem....
  • H. Elman et al.

    Finite elements and fast iterative solvers with applications to incompressible fluid dynamics

    Numerical Mathematics and Scientific Computation

    (2005)
  • C. Fletcher

    Computational Techniqus for Fluid Dynamics, vol. 2

    (1991)
  • P. Fullsack

    An arbitrary Lagrangian–Eulerian formulation for creeping flows and its application in tectonic models

    Geophysical Journal International

    (1995)
  • T.V. Gerya et al.

    Robust characteristics method for modelling multiphase visco-elasto-plastic thermo-mechanical problems

    Physics of the Earth and Planetary Interiors

    (2007)
  • W. Gorczyk et al.

    Growth and mixing dynamics of mantle wedge plumes

    Geology

    (2007)
  • R.S. Huismans et al.

    Asymmetric lithospheric extension: the role of frictional plastic strain softening inferred from numerical experiments

    Geology

    (2002)
  • P. Jäger et al.

    On local tracking algorithms for the simulation of three-dimensional

    Comput. Mech.

    (2008)
  • R. Jamieson et al.

    Provenance of the greater himalayan sequence and associated rocks: predictions of channel flow models

  • R.F. Katz et al.

    The dynamics of melt and shear localization in partially molten aggregates

    Nature

    (2006)
  • B. Kaus et al.

    Eulerian spectral/finite difference method for large deformation modelling of visco-elasto-plastic geomaterials

    Bolletino Geofisico

    (2004)
  • Cited by (58)

    • Numerical Modeling of Subduction

      2023, Dynamics of Plate Tectonics and Mantle Convection
    • Stress recovery for the particle-in-cell finite element method

      2021, Physics of the Earth and Planetary Interiors
      Citation 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 Interiors
      Citation 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 Flow
      Citation 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.

    View all citing articles on Scopus
    View full text