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

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.