Viscoelastic formulation for modelling of plate tectonics

Louis Moresi, Frédéric Dufour, Hans Mühlhaus

► To cite this version:

Louis Moresi, Frédéric Dufour, Hans Mühlhaus. Viscoelastic formulation for modelling of plate tectonics. Bifurcation and Localisation Theory in Geomechanics, 2001, Perth, Australia. 10.1201/9781003210931-46. hal-01009106

HAL Id: hal-01009106

https://hal.science/hal-01009106

Submitted on 19 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Viscoelastic formulation for modeling of plate tectonics.

Louis Moresi, Frédéric Dufour, Hans Mühlhaus

CSIRO Exploration and Mining, Perth, 6009, AUSTRALIA

**ABSTRACT:** The Earth's tectonic plates are strong, viscoelastic shells which make up the outermost part of a thermally convecting, predominantly viscous layer. In order to build a more realistic simulation of the planet's evolution, the complete viscoelastic convection system must be included. A particle-in-cell finite element method is demonstrated which can simulate very large deformation viscoelasticity. This is applied to a plate-deformation problem. Numerical accuracy is demonstrated relative to analytic benchmarks, and the characteristics of the method are discussed.

## 1 INTRODUCTION

Underneath the lithospheric plates of the Earth lies the mantle (Figure 1). Approximately 3000km deep, it is composed of solid rock that is warm enough to deform like a viscous fluid, albeit at incredibly slow speeds of a few centimetres per year. The plates move because the mantle is forever stirring as heat generated by natural radioactive decay struggles to escape via thermal convection. The plates which form the ocean floors are part of this circulation and are sucked down when they become old, cold and dense. The continental crust is formed by lower density rock which remains buoyant despite being cold. In the lithosphere the rocks are significantly cooler and behave as a viscoelastic, brittle solid. In regions of high stress, brittle failure gives rise to earthquakes.

This picture of the Earth's interior is widely accepted by geophysicists. It clearly indicates that the fundamental process is thermal convection; plate tectonics is the manner in which the system organizes. Therefore, a consistent model of plate behaviour must contain a description of the convection system of which the plate is a part.

There are some fundamental problems which need to be addressed before the routine application of engineering principles to the lithosphere. The principle issues is that plate tectonics is itself only a kinematic description of the observations: a fully consistent dynamic description of the motion of the plates is still sought.

There have been some major steps towards the simulation of plate tectonics in recent years by solving brittle/viscous fluid flow equations (e.g. Tackley 1998, 2000, Moresi & Solomatov, 1998). However, incorporating these developments into an engineering description of the lithosphere has proved to be difficult because of the absence of elastic stresses from the simulations. In the past, viscoelastic con-

vection simulations have been limited to models with explicit layering in which a non-convecting viscoelastic layer is coupled to a viscous convecting domain (Podladchikov et al, 1993). Models of subduction zones which incorporate viscoelasticity, faulting, and free-surface behaviour have generally been limited to modest evolution times, after which further deformation produces severe remeshing problems (e.g. Melosh, 1978, Gurnis et al, 1996). The Natural Element Method of Braun & Sambridge (1995) is capable of handling very large deformation viscoelasticity, but with considerable complexity, particularly in the extension to 3D. Having identified the need for efficient, large-scale convection simulations with elastic effects in an evolving cool lithosphere, we present a method for simulating viscoelastic-brittle materials in extreme deformation.

## 2 MATHEMATICAL MODEL

We begin our analysis in a general way with the classical momentum conservation equation:

<figure><img src="image_2.png" /><figcaption>Figure 1: A simplified cross section of the Earth with major layerings shown to scale except for the upper boundary layer which is exaggerated in thickness by a factor of roughly two.</figcaption></figure>
$$ \nabla \cdot \boldsymbol{\sigma} = \boldsymbol{f} \qquad (1) $$

where $\boldsymbol{\sigma}$ is the stress tensor and $\boldsymbol{f}$ a force term. As we are interested only in very slow deformations of highly viscous materials, (infinite Prandlt number) we have neglected all intertial terms in (1). It is convenient to split the stress into a deviatoric part, $\boldsymbol{\tau}$, and an isotropic pressure, $\boldsymbol{p}$,

$$ \boldsymbol{\sigma} = \boldsymbol{\tau} - \boldsymbol{p}\boldsymbol{I} \qquad (2) $$

where $\boldsymbol{I}$ is the identity tensor.

## 2.1 Viscoelasticity

There are a number of different viscoelastic models, we will use the Maxwell model which has been used in previous studies of lithospheric deformation where viscous and elastic effects are important such as post-glacial rebound (Peltier, 1974). This model assumes that the strain rate tensor, $\boldsymbol{D}$, defined as:

$$ D_{ij} = \frac{1}{2} \left( \frac{\partial V_i}{\partial x_j} + \frac{\partial V_j}{\partial x_i} \right) \qquad (3) $$

is the sum of an elastic strain rate tensor $D_e$ and a viscous strain rate tensor $D_v$. The velocity vector, $\boldsymbol{V}$, is the fundamental unknown of our problem and all these entities are expressed in the fixed reference frame $\boldsymbol{x}_i$. Now we decompose each strain rate tensor

$$ D_e = \frac{1}{3} \mathrm{tr}(D_e)\boldsymbol{I} + \hat{D}_e \qquad (4) $$

and

$$ D_v = \frac{1}{3} \mathrm{tr}(D_v)\boldsymbol{I} + \hat{D}_v \qquad (5) $$

where $\hat{D}$ is the deviatoric part of $\boldsymbol{D}$ and $\mathrm{tr}(\boldsymbol{D})$ represents the trace of the tensor.

Individually we express each deformation tensor as a function of the deviatoric stress tensor $\boldsymbol{\tau}$ and pressure $\boldsymbol{P}$, which finally gives a tensorial equation:

$$ \frac{\nabla \boldsymbol{\tau}}{2\mu} + \frac{\boldsymbol{\tau}}{2\eta} = \hat{D}_e + \hat{D}_v = \hat{D} \qquad (6) $$

where $\frac{\nabla \boldsymbol{\tau}}{2\mu}$ is the Jaumann corotational stress rate for an element of the continuum, $\mu$ is the shear modulus and $\eta$ is shear viscosity. The isotropic part gives a scalar equation for the pressure:

$$ \frac{\dot{p}}{K_e} + \frac{p}{\xi} = -\mathrm{tr}(\boldsymbol{D}) \qquad (7) $$

where $K_e$ is the bulk modulus and $\xi$ is the bulk viscosity. $\dot{p} \equiv \dot{p}$ as it $p$ is a scalar.

$$ \frac{\dot{\boldsymbol{\tau}}}{\boldsymbol{\tau}} = \dot{\boldsymbol{\tau}} + \boldsymbol{\tau}\boldsymbol{W} - \boldsymbol{W}\dot{\boldsymbol{\tau}} \qquad (8) $$

where $\boldsymbol{W}$ is the material spin tensor,

$$ W_{ij} = \frac{1}{2} \left( \frac{\partial V_i}{\partial x_j} - \frac{\partial V_j}{\partial x_i} \right) \qquad (9) $$

The $\boldsymbol{W}$ terms account for material spin during advection which reorients the elastic stored-stress tensor.

We note that the form of equation (7) is unsuited to conventional fluids as the material has no long term resistance to compression. This behaviour is, however, relevant to the simulation of the coupled porous-flow, matrix deformation problem. Here it is common to ascribe an apparent bulk viscosity to the matrix material in order to model compaction effects (e.g. McKenzie 1984), particularly for large scale geophysical systems where the details of the pore network cannot be measured directly.

## 2.2 Numerical implementation

As we are interested in solutions where very large deformations may occur — including thermally driven fluid convection, we would like to work with a fluid-like system of equations. Hence we obtain a stress / strain-rate relation from (6) by expressing the Jaumann stress-rate in a difference form:

$$ \frac{\dot{\boldsymbol{\tau}}}{\boldsymbol{\tau}} \approx \frac{\boldsymbol{\tau}^{t+\Delta t} - \boldsymbol{\tau}^t}{\Delta t} - \boldsymbol{W}^t \boldsymbol{\tau}^t + \boldsymbol{\tau}^t \boldsymbol{W}^t \qquad (10) $$

where the superscripts $t$, $t + \Delta t$ indicate values at the current and future timestep respectively. (6) and (7) become respectively

$$ \boldsymbol{\tau}^{t+\Delta t} = \frac{\eta \Delta t}{\alpha + \Delta t} \hat{D}^{t+\Delta t} + \frac{\alpha}{\alpha + \Delta t} \boldsymbol{\tau}^t + \frac{\alpha \Delta t}{\Delta t + \alpha} (\boldsymbol{W}^t \boldsymbol{\tau}^t - \boldsymbol{\tau}^t \boldsymbol{W}^t) \qquad (11) $$

and

$$ p^{t+\Delta t} = -\frac{\xi \Delta t}{\beta + \Delta t} D_{kk}^{t+\Delta t} + \frac{\beta}{\beta + \Delta t} p^t \qquad (12) $$

where $\alpha = \eta/\mu$ is the shear relaxation time and $\beta = \xi/K_e$ is the bulk relaxation time. We can simplify the above equations by defining an effective viscosity $\eta_{\mathrm{eff}}$ and an effective compressibility $\xi_{\mathrm{eff}}$:

$$ \eta_{\mathrm{eff}} = \eta \frac{\Delta t}{\Delta t + \alpha} \quad \text{and} \quad \xi_{\mathrm{eff}} = \xi \frac{\Delta t}{\Delta t + \beta} \qquad (13) $$

Then the deviatoric stress is given by

$$ \boldsymbol{\tau}^{t+\Delta t} = \eta_{\mathrm{eff}} \left( \hat{D}^{t+\Delta t} + \frac{\boldsymbol{\tau}^t}{\mu \Delta t} + \frac{\boldsymbol{W}^t \boldsymbol{\tau}^t}{\mu} - \frac{\boldsymbol{\tau}^t \boldsymbol{W}^t}{\mu} \right) \qquad (14) $$

and the pressure by

$$ p^{t+\Delta t} = -\xi_{\mathrm{eff}} \left( D_{kk}^{t+\Delta t} - \frac{p^t}{\Delta t K_e} \right) \qquad (15) $$

To model an incompressible material $K_e$ and $\xi$ are made very large such that $D_{kk} \approx 0$.

Our system of equations is thus composed of a quasi-Newtonian viscous part with modified mate-
rial parameters and a right-hand-side term depending on values from the previous timestep. This approach minimizes the modification to the viscous flow code. Instead of using physical parameters for viscosity and bulk modulus, we use effective material properties (13) to take into account elasticity. Then during computations for the force term, we add elastic internal stresses from the previous timestep or from initial conditions.

$$F_{i,t}^{e,t} = \frac{\xi_{\text{eff}}}{K_e \Delta t} p_{i,t}^{t-\Delta t} - \frac{\eta_{\text{eff}}}{\mu \Delta t} \sigma_{ij,j}^{t-\Delta t} \quad (16)$$

We solve (15) and (14) and obtain a solution for $D^{t+\Delta t}$. From this solution we compute the new stress state due to the velocity field and previous stored stresses.

## 2.3 Stability

The approach outlined above is unconditionally stable only if the timestep is larger than the relaxation time for the material, i.e.

$$\Delta t < \frac{\eta}{\mu} \quad (17)$$

in the case of the shear moduli. Alternatively, this means a Deborah number, $De < 1$, indicating that the method is appropriate to the viscous, rather than the elastic, limit.

One difficulty is that the timestep is not necessarily chosen to match the physical problem, but by the Courant condition for the chosen mesh. This means that a convergence demonstration for arbitrarily small elements may not be possible for the general case. We are currently addressing this issue.

In practice, however, for our area of research, viscous flow drives the plate motions, and the lithospheric plates are embedded in a highly viscous material. This may produce a situation where a system-wide relaxation time is more important than the relaxation times of individual materials, since loading and unloading of the elastic materials happens almost exclusively through a low-viscosity medium. Under these circumstances, the relaxation time of an individual layer (such as the lithosphere) may be much larger than the Courant timestep, but stresses are either balanced, or relaxed by driving a flow in one of the viscous materials.

# 3 COMPUTATIONAL METHOD

## 3.1 Choice of Numerical Scheme

In fluid dynamics, where strains are generally very large, but not important in the constitutive relationship of the material, it is common to transform the equations to an Eulerian mesh and deal with convective terms explicitly. Problems arise whenever advection becomes strongly dominant over diffusion since an erroneous numerical diffusion dominates. In our case, the advection of material boundaries and the stress tensor are particularly susceptible to this numerical diffusion problem. Mesh-based Lagrangian formulations alleviate this difficulty, but at the expense of remeshing and the eventual development of a

less-than optimal mesh configuration. This increases complexity and can hinder highly efficient solution methods such as multigrid iteration. The Natural Element Method eliminates remeshing difficulties but is associated with considerable complexity of implementation, particularly in 3D.

A number of alternatives are available which dispense with a mesh entirely: smooth particle hydrodynamics and discrete element methods are common examples from the fluid and solid mechanics fields respectively. These methods are extremely good at simulating the detailed behaviour of highly deforming materials with complicated geometries (e.g. free surfaces, fracture development), and highly dynamic systems. They are, in general, formulated to calculate explicitly the interactions between individual particles which ultimately means that a great many timesteps would be required to study creeping flow where the timescales associated with inertial effects are very many orders of magnitude smaller than typical flow times.

We have therefore developed a hybrid approach – a particle in cell finite element method which uses a standard Eulerian finite element mesh (for fast, implicit solution) and a Lagrangian particle framework for carrying details of interfaces, the stress history etc.

## 3.2 The Particle in Cell Approach

Our particle-in-cell finite element method is based closely on the standard finite element method, and is a direct development of the material point method of Sulsky et al. (1995). The standard mešh is used to discretize the domain into elements, and the shape functions interpolate node points in the mesh in the usual fashion. The problem is formulated in a weak form to give an integral equation, and the shape function expansion produces a discrete (matrix) equation. Equation (1) in weak form, using the notation of (2) becomes

$$\int_{\Omega} N_{(i,j)} \tau_{ij} d\Omega - \int_{\Omega} N_i p d\Omega = \int_{\Omega} N_i f_i d\Omega \quad (18)$$

where the trial functions, $N$, are the shape functions defined by the mesh, and we have assumed no non-zero traction boundary conditions are present. For the discretized problem, these integrals occur over subdomains (elements) and are calculated by summation over a finite number of sample points within each element. For example, in order to integrate a quantity, $\phi$ over the element domain $\Omega^e$ we replace the continuous integral by a summation

$$\int_{\Omega^e} \phi d\Omega \leftarrow \sum_p w_p \phi(x_p) \quad (19)$$

In standard finite elements, the positions of the sample points, $x_p$, and the weighting, $w_p$, are optimized in advance. In our scheme, the $x_p$'s correspond precisely to the Lagrangian points embedded in the fluid, and $w_p$ must be recalculated at the end of a timestep for the new configuration of particles. Constraints on the values of $w_p$ come from the need to integrate polynomials of a minimum degree related to the degree of the shape function interpolation, and the order of the
underlying differential equation (e.g Hughes, 1987). These Lagrangian points carry the history variables which are therefore directly available for the element integrals without the need to interpolate from nodal points to fixed integration points. In our case, the distribution of particles is usually not ideal, and a unique solution for $w_p$ cannot be found, or we may find we have negative weights which are not suitable for integrating physical history variables. We therefore store an initial set of $w_p$'s based on a measure of local volume and adjust the weights slightly to improve the integration scheme.

Moresi et al. (2000) give a full discussion of the implementation of the particle-in-cell finite element scheme used here including full details of the integration scheme and its assumptions. They also discuss the specific modifications to the material point method required to handle a convecting fluid.

# 4 BENCHMARKS

We have benchmarked our numerical scheme against analytic solutions in order to characterize its strengths and weaknesses, and to quantify the likely level of accuracy we can achieve with a given mesh/particle density. We first benchmark the purely viscous flow case to provide a baseline for comparison with viscoelastic cases.

## 4.1 Analytical solution

We study the spreading of a rectangular sample of material under a constant downward velocity $V$ applied on top (see Fig.2).

The specified boundary conditions give:

$$ \sigma_{xx} = \tau_{xx} - p = 0 \qquad (20) $$

which is substituted into the $x$ component of (6) to give an expression in the pressure

$$ \frac{\dot{p}}{2\mu} + \frac{p}{2\eta} = \frac{1}{2}(D_{xx} - D_{zz}) \qquad (21) $$

$D_{xx}$ can be eliminated between (21) and (7), and $D_{zz}$ is given kinematically as

<figure><img src="image_3.png" /><figcaption>Figure 2: Geometry and boundary conditions of for the analytic solution.</figcaption></figure>

$$ D_{zz} = -\frac{V}{h_0 - Vt} \qquad (22) $$

where $h_0$ is the initial height. Hence

$$ a \dot{p}(t) + b p(t) = - \frac{2V}{h_0 - Vt} \qquad (23) $$

with

$$ a = -\frac{1}{K_e} - \frac{1}{\mu} \quad \text{and} \quad b = -\frac{1}{\xi} - \frac{1}{\eta} \qquad (24) $$

The solution to this equation is

$$ p(t) = \frac{2e^{\frac{b}{a}(\frac{h_0}{V}-t)}}{a} \times \left\{ \ln\left(1-\frac{V}{h_0}t\right) + \sum_{i=1}^{\infty} \left(-\frac{b h_0}{a V}\right)^i \frac{(1-\frac{V}{h_0}t)^i - 1}{i \cdot i!} \right\} \qquad (25) $$

We use this relationship to eliminate the pressure derivative from (21) and (7) and express the unknown $D_{xx}$ as a function of $p$ and $D_{zz}$. Consequently we obtain a relation between the velocity $V_x$ of a point and its coordinate $x$.

$$ V_x = -\frac{x}{K_e + \mu} \times \left( p \left( \frac{K_e}{\xi} - \frac{\mu}{\eta} \right) + D_{zz} (K_e - \mu) \right) \qquad (26) $$

## 4.2 Model

The Eulerian mesh does not carry any information from timestep to timestep other than the boundary conditions. Therefore, when convenient, the mesh may be modified, replaced, and refined as necessary. For this problem, compression is applied by a moving boundary condition which causes the mesh to compact in one direction. For simplicity, the mesh is simply scaled to the new aspect ratio without altering the number of elements. In a more complicated situation, however, it would be possible to regrid completely without loss of accuracy. The only detail which needs to be observed is that the updating of the boundary node locations follows the same formulation as that of the particles (here, a second order Runge-Kutta integration procedure) to prevent the boundary conditions from drifting with respect to the stored information on the particles.

In the $x$-direction we have a free surface. In order to investigate the properties of a particle representation of such interfaces, it is important not to simply use a mesh-based boundary condition. Instead we use a mesh (Fig. 2) larger than the specimen and fill the gap with a background material having ($\eta = 10^3$ MPa.s and $\xi = 10^5$ MPa.s).

The square mesh is composed of 4096 elements.
## 4.3 Results

We have tested the code, ELLIPSIS, with three different type of materials: viscous and compressible Fig.(3.c and d), viscous and incompressible Fig.(3.a and b), viscoelastic and compressible Fig.(4). The numerical parameters for the viscous part of (15) and (14) are:

$$ \eta = 10^6 \text{ MPa.s.} \quad \text{and} \quad \xi = 4.10^6 \text{ MPa.s.} \qquad (27) $$

Due to the mesh size and the prescribed velocity on top, the Courant condition requires $\Delta t < 5.10^{-3}$ s.. For our problem we take $\Delta t = 10^{-3}$ s..

As an indicator of accuracy we compare analytical and numerical x-velocity at the point I (Fig. 2). The steps in the numerical solutions are related to the motion of the material interface relative to the element edges. When the interface between the specimen and the background material crosses into a new element there is an immediate discontinuous contribution to the element equations from the sample material. We have verified that jump tends to zero as we increase the element and particle densities.

In the viscoelastic case, in Fig. (4) we can see that, for different relaxation time, the error against the analytical solution remains below 3%. Some steps are still present, but in the viscoelastic case the velocity is considerably more noisy. What is most clear is that the computations with larger relaxation time have greater fluctuation in accuracy. The problem becomes more acute in this case when the Courant time (decreasing due to the compression of the background mesh) becomes comparable to the relaxation time. This can result in a loss of stability which is entirely an artefact of the discretization. The most promising solution to this issue is to compute the time-derivative of the stress tensor over a physically relevant interval, rather than that imposed by the mesh. This is a particular focus of our current research.

<figure><img src="image_4.png" /><figcaption>Figure 3: Viscous deformation results: a. viscous incompressible, analytical, b. viscous incompressible, numerical, c. viscous compressible, analytical and d. viscous compressible, numerical.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Viscoelastic material a. Numerical solution, b. $\mu = 10^4$ MPa., c. $\mu = 10^5$ MPa.</figcaption></figure>

# 5 APPLICATION TO PLATE DYNAMICS

## 5.1 Mathematical Model

We treat the Earth on a large scale as an incompressible, viscoelastic Maxwell fluid with infinite Prandtl number in which motions are driven by internal temperature variations. The force term from equation (1) is a gravitational body force due to density changes. We assume that these arise, for any given material, through temperature effects:

$$ \nabla \cdot \tau - \nabla p = g\rho_0(1 - \alpha T)\hat{z} \qquad (28) $$

where g is the acceleration due to gravity, $\rho_0$ is material density at a reference temperature, $\alpha$ is the coefficient of thermal expansivity, and T is temperature. $\hat{z}$ is a unit vector in the vertical direction. We have also assumed that the variation in density only needs to be considered in the driving term (the Boussinesq approximation).

The equation of motion is then

$$ \nabla(\eta_{\text{eff}}D^{t+\Delta t}) - \nabla p = g\rho_0(1 - \alpha T)\hat{z} - \nabla \left( \eta_{\text{eff}} \left[ \frac{\tau^t}{\mu\Delta t} + \frac{W^t\tau^t}{\mu} - \frac{\tau^t W^t}{\mu} \right] \right) \quad (29) $$

The velocity field u and pressure at $t + \Delta t$ can be solved for a given temperature distribution and the stress history from the previous step.

Motion is driven by the heat escaping from the interior. The energy equation governs the evolution of the temperature in response to diffusion of heat through the fluid. For a given element of fluid,

$$ \frac{DT}{Dt} = -\kappa \nabla^2 T \qquad (30) $$

where $\kappa$ is the thermal diffusivity of the material.

So far, all equations have been written in a purely
Lagrangian framework. The time derivate of temper-
ature and the Jaumann stress rate refer to a frame of
reference which is carried by the fluid. In choosing a
solution method, it is necessary to choose whether to
honour the Lagrangian formulation, or to work with a
fixed reference frame and introduce additional terms
to compensate for the advection of temperature and
stress by the fluid.

5.2 *Brittle failure*

As we discussed above, plate models need to in-
clude a description of the brittle nature of the cold-
est part of the lithosphere. Geologists use this term
quite loosely to distinguish fault-dominated deforma-
tion which may result in seismic activity, from ductile
creep which occurs at higher temperature and pres-
sure. In all recent studies of mantle convection where
the brittle lithospheric rheology has been taken into
account, the brittle behaviour has been parameterized
using a non-linear effective viscosity which is intro-
duced whenever the stress would otherwise exceed
the yield value τyield. This approach ignores details
of individual faults, and treats only the influence of
fault systems on the large-scale convective flow.

To determine the effective viscosity we extend (6)
by introducing a von Mises plastic flow rule:

$$
\frac{\nabla}{2\mu} + \frac{\tau}{2\eta} + \lambda \frac{\tau}{2|\tau|} = \hat{D}_e + \hat{D}_v + \hat{D}_p = \hat{D} \quad (31)
$$

where λ is a parameter to be determined such that
the stress remains on the yield surface, and |τ| ≡
(τᵢⱼτᵢⱼ/₂)⁽¹⁾⁄₂. We again express the Jaumann stress
rate in difference form (10) to give:

$$
\tau^{t+\Delta t} \left[ \frac{1}{2\mu\Delta t} + \frac{1}{2\eta} + \frac{\lambda}{2|\tau|} + \right] = \\
\dot{D}^{t+\Delta t} + \frac{1}{2\mu\Delta t}\tau' + \frac{1}{2\mu}(\mathbf{W}'\tau^t - \tau^t\mathbf{W}^t) \quad (32)
$$

No modification to the isotropic part of the problem is
required when the von Mises yield criterion is used.
At yield we use the fact that |τ| = τyield to write

$$
\tau^{t+\Delta t} = \eta' \left[ 2\hat{D}^{t+\Delta t} + \frac{1}{\mu\Delta t}\tau^t + \frac{1}{\mu}(\mathbf{W}^t\tau^t - \tau^t\mathbf{W}^t) \right] \quad (33)
$$

using an effective viscosity, η' given by

$$
\eta' = \frac{\eta \tau_{\text{yield}} \mu \Delta t}{\eta \tau_{\text{yield}} + \tau_{\text{yield}} \mu \Delta t + \lambda \eta \mu \Delta t} \quad (34)
$$

We determine λ by equating the value of |τ<sup>t+Δt</sup>| with
the yield stress in (33). Alternatively, in this particular
case, we can obtain η' directly as

$$
\eta' = \tau_{\text{yield}} / |\hat{D}_{\text{eff}}| \quad (35)
$$

where

<figure><img src="image_6.png" /><figcaption>Figure 5: Example: compression of a viscoelastic plate with yield stress overlying a low viscosity fluid of equal density</figcaption></figure>

$$
\hat{D}_{\text{eff}} = 2\hat{D}^{t+\Delta t} + \frac{1}{\mu\Delta t}\tau^t + \frac{1}{\mu}(\overbrace{\mathbf{W}^t\tau^t - \tau^t\mathbf{W}^t}^{\text{---}}) \quad (36)
\text{and } |\mathbf{D}| = (2D_{ij}D_{ij})^{1/2}.
$$

The value of λ or η' is iterated to allow stress to
redistribute from particles which become unloaded.
The iteration is repeated until the velocity solution is
unchanged to within the error tolerance required for
the solution as a whole.

The value of the yield stress is, in principle, a func-
tion of strain, yield history, and temperature, and can
be distinct for different materials.

6 PLATE MODELING

As a simple example, we demonstrate the compres-
sion of a viscoelastic-brittle layer which lies on top
of a slightly less dense viscous fluid layer (Figure
5). This system is an analogue of the cool oceanic
lithosphere which rests upon the warm asthenosphere
(though we do not solve the temperature equation
in this case). The viscoelastic layer is initially split
to provide an initiation point for a model subduc-
tion zone. The vertical boundaries are free-slip, and
the right hand edge is given a horizontal velocity to
shorten the system. There is a layer of highly com-
pressible material of very low viscosity above the
elastic layer which accomodates the volume change
associated with shortening of the mesh, and mimics a
free surface boundary condition on the upper surface
of the elastic layer. As compression proceeds, the vis-
coelastic layer flexes and the viscous layer flows to
accomodate the deformation. As stresses build up in
the model lithosphere, a second failure point develops
allowing one half of the material to fold up under the
other half. Further compression forces the two halves
of the lithosphere layer to slide past each other along
a zone of material failure. After this point, the pres-
ence of the bottom boundary begins to interfere with
the evolution of the system.

This particular simulation demonstrates the capa-
bility of the algorithm in the simulation of subduction
zone geometry in the style of Melosh (1978) or Gur-
nis et al, (1996). The fact that the algorithm is imple-
mented within a fluid-dynamics framework suggests
that a viscoelastic analysis of convection with strong
temperature dependence of viscosity and a yield stress
is now possible.

7 DISCUSSION

The algorithm described above is designed to intro-
duce elastic effects into convection simulations where
temperature-dependent viscosity and yielding domi-
nate the mechanical behaviour. The viscosity of the
mantle and the mantle lithosphere is very strongly de-
pendent on temperature (several orders of magnitude
variation over 1000°C) whereas the shear modulus is
not strongly affected (there is only a modest change
in seismic wavespeed due to temperature). Therefore,
elastic effects become unimportant outside the cold
thermal boundary layer where viscosity is extremely
large.

The influence of elastic stresses is likely to be felt
at the subduction zones where the lithosphere is bent
into the interior of the Earth. In these regions stresses
are typically close to the yield stress — a fact which
allows the plates to move in the first place. The ad-
dition of elasticity is likely to complicate the simple
picture presented by Tackley (1998) and Moresi &
Solomatov (1998) for viscous materials with a yield
stress.

Our methodology is limited to a coarse continuum
description of the subduction zone system at a resolu-
tion of a few km. This may be able to give us valu-
able information into the nature of plate tectonics, the
thermal conditions in and around subducting litho-
sphere, and the stress state of the system. However,
the resolution is too coarse to say anything about the
detailed mechanics of the failure of lithospheric fault
zones and the conditions for major failure to occur.
For this we require a coupling of the large-scale code
with an engineering-scale code (e.g. DEM or small-
deformation Lagrangian FEM) using the large-scale
to provide boundary conditions for the small scale.
The issue of scale-bridging is important in many ar-
eas of numerical simulation. Essentially the same dif-
ficulties arise in material science where the atomic
scale is best treated by molecular dynamics codes but
the large scale must be treated as a continuum (e.g.
Bernholc, 1999).

ACKNOWLEDGEMENTS

The authors gratefully acknowledge the suggestions
of an anonymous reviewer which have helped to clar-
ify this paper.

REFERENCES

Bernholc, J., 1999, Computational materials science: The
era of applied quantum mechanics, Physics Today, 52, No
9, 30–35.

Braun, J., Sambridge, M., 1995, A numerical method for solving partial differential equations on highly irregular evolving grids, Nature, 376, 655-660.

Gurnis, M., Eloy, C., Zhong, S., 1996, Free-surface formulation of mantle convection — II. Implications for subduction zone observables, Geophys. J. Int., 127, 719–727.

Hughes, T. J.R., 1984, The Finite Element Method,
Prentice-Hall.

Melosh, H. J., 1978, Dynamic support of the outer rise, Geophys. Res. Lett., 5, 321–324.

Moresi, L., Mühlhaus, H.-B., Dufour, F., Proceedings of the 5th International Workshop on Bifurcation and Localization in Geomechanics, Perth, W. A., Australia, Balkema.

Moresi, L., Solomatov, V. S., 1998, Mantle convection with a brittle lithosphere: thoughts on the global tectonics styles of the Earth and Venus, Geophys. J. Int., 133, 669–682.

McKenzie, D., 1984, The generation and compaction of partially molten rock, J. Petrology, 25, 713–765.

Peltier, W. R., 1974, The impluse response of a Maxwell
Earth, Rev. Geophys. Space Phys. 12, 649–669.

Podlachikov, Yu., Yu., Lenardic, A., Yuen, D. A., Quareni, F., 1993, Dynamical consequences of stress focussing for different rheologies: Earth and Venus perspectives, EOS Trans. AGU, 74 no. 43/Suppl., 566.

Sulsky, D., Zhou, S.-J., Schreyer, H. L., 1995, Application of a particle-in-cell method to solid mechanics, Comput. Phys. Commun. 87, 236-252.

Tackley, P. J., 1998, Self-consistent generation of tec-
tonic plates in three dimensional mantle convection, Earth
Planet. Sci. Lett., 157, 9-22.

Tackley, P. J., 2000, The quest for self-consistent genera-
tion of plate tectonics in mantle convection models, AGU
Monograph on The history and dynamics of global plate
motions, ed. M. Richards, In Press.