L Moresi

F Dufour

H B M Ühlhaus

A Lagrangian integration point finite element method for large deformation modeling of viscoelastic geomaterials

We review the methods available for largedeformation simulations of geomaterials before presenting a La grangian integration point finite element method de signed specifically to tackle this problem. In our EL LIPSIS code, the problem domain is represented by an Eulerian mesh and an embedded set of Lagrangian in tegration points or particles. Unknown variables are computed at the mesh nodes and the Lagrangian parti cles carry history variables during the deformation pro cess. This method is ideally suited to model fluidlike behavior of continuum solids which are frequently en countered in geological contexts. We present bench mark examples taken from the geomechanics area.

Introduction

A longterm goal of geological modeling is to achieve the degree of simulation capability currently enjoyed by the engineering community. The routine ability to recreate the evolution of geological structures dur ing deformation and simultaneously compute stresses, temperatures, fluid flow vectors and chemical evolu tion would greatly enhance our ability to understand the Earth. Many questions in geology are formulated in terms of inverse analysis because of the very nature of the science but there is little hope of making significant progress until a reliable forward modeling capability is developed (Wijns et al., 2001 [20]). Unlike most engineering simulations, however, the ge ometry of the geological model is a result of the contin uously evolving nonlinear interaction of the structure and the rheology at extremely large strains (M ühlhaus et al, 2001 [START_REF] Moresi | The interplay of material and geometric instabilities in large deformations of viscous rock[END_REF]). Initially the geometry might be rela tively simple, e.g. flat layering, but during the course of the simulation very intricate patterns will develop and need to be accurately resolved by the numerical method. In many cases the complexity reflected in the geometry results from postfailure deformation of the material. To make matters worse, the rheological laws which govern particular geological materials at a given scale are often poorly determined this is partly due to the difficulty in measuring behaviour at realistic strain rates, and partly due to the microstructural complex ity of the materials which make up a given rock suite. The fact that the materials also undergo metamorphic transformations and phase changes during the evolu tion of a geological structure is another complicating factor which makes simulation difficult.

As an example, on a global scale, the continents drift as an integral part of the surface thermal boundary layer of the convecting mantle. They have retained a dis tinct identity within the mantle flow for billions of years while developing a strong physical and chemical fabric along the way. Motions in the mantle are described by the equations of fluid dynamics for very large deforma tion. The rheology needed to describe deformation in the lithosphere is highly nonlinear, and near the sur face where temperatures are less than approximately 600 • C it becomes necessary to consider the role of elas ticity (Watts et al, 1980 [19]). The strong correlation be tween seismicity and plate boundaries (e.g. [START_REF] Barazangi | World seismicity maps compiled from essa, coast and geodetic sur vey, epicenter data 19611967[END_REF] [START_REF] Barazangi | World seismicity maps compiled from essa, coast and geodetic sur vey, epicenter data 19611967[END_REF]) makes it seem likely that plate mo tions are associated with localization of deformation occuring when stresses reach the yield strength of the lithosphere.

From a modeling point of view, it is necessary to consider the fluid convection of the mantle and the historydependent viscoelastic/brittle behaviour of the continental crust as a single coupled system. At the same time, the precise structure and composi tion of the deep continental crust is not well known. The requirements for a geological simulation code are therefore an ability to track boundaries and interfaces through extremely large deformation, including fluid convection, of nonlinear history dependent materi als. The wide range of physical and temporal scales, and the many coupled physical processes also impose a need for computational efficiency. The code should also be very flexible in the rheological laws which it can treat.

Many different numerical methods have been devised for mechanical simulations of this kind. Some derive from standard engineering methods, while others were developed to handle specific problems in the physical sciences. We summarize a number of these methods in order to illustrate the difficulties involved in creating convincing, realistic simulations. We next present the ELLIPSIS Lagrangian Integration Point code: a method for simulating viscoelasticbrittle materials in fluidlike deformation. The method is tested on a number of very simple benchmark cases where analytic solutions are known, or where the accuracy can otherwise be quan tified in large deformation.

Figure 1: Different ways to discretize a problem pro vide a natural representation for systems with differ ent controlling physics. In fully Lagrangian FEM, SPH and DEM all the computational points are also material points. In the Lagrangian integration point method, computational points are not material points but a set of material points is also tracked. In ALE and Eulerian FEM there is no tracking of material points.

modeling techniques

The key to all these problems is to find a way to deal efficiently with the most general case: finite strain vis coelasticity. This is a specialized problem: most en gineering codes are optimized to study the modest strains which accumulate prior to failure. Neverthe less, an enormous amount of research has been done in developing numerical methods suited for large strain problems in different application areas. The principal difficulty in finite strain modeling is the tendency for any computational or logical connecting mesh between material points to become arbitrarily tangled during the simulation. Broadly speaking there have been two approaches to dealing with this problem. One approach is to retain the mesh, but to reorganize it periodically and, per haps, modify the method to cope with a lessthan optimal distribution of discretization points, and the other is to dispense with the mesh entirely for com puting mechanical behavior and derive a formulation where the material points interact with each other di rectly. The approach used in Ellipsis is distinct from these methods in relaxing the requirement that mate rial points and computational points need to be the same (see Figure 1)

Highly deforming mesh methods

Mesh based Lagrangian methods are not particularly well suited to very large strain fluid flow applications as the mesh is subject to tangling and considerable effort is required to regrid or reconnect the mesh to prevent the computation of derivatives and shape functions from becoming inaccurate. The Arbitrary Lagrangian Eulerian method (Huerta and Liu 1987 [START_REF] Heurta | Viscous flow with large free surface motion[END_REF]) avoids mesh tangling by allowing computational points to move in dependently from the underlying material. This strat egy makes it possible to prevent mesh points from ap proaching one another too closely, while still preserv ing important interface details. However, additional advection terms are required to handle transport of quantities relative to the mesh and these terms are dis persive. The Dynamical Lagrangian Remeshing method (Braun and Sambridge 1994 [START_REF] Braun | Dynamical lagrangian remeshing (dlr): A new algortihm for solving large strain deformation problems and its application to faultpropagation folding[END_REF]) reconnects all mesh points af ter any deformation to ensure the mesh always remains optimally configured (in DLR this means the mesh is always a Delaunay triangulation). Stress histories can be compute at integration points and therefore need to be interpolated when new element connections occur. The extension to 3D is nontrivial and the method is limited to linear triangles/tetrahedra. An alternative approach, the Natural Element Method (NEM, Braun and Sambridge 1995 [START_REF] Braun | A numerical method for solving partial differential equations on highly irregular evolving grids[END_REF]) uses a natural neighbor interpolation scheme over a Delaunay trian gulation to develop finiteelementlike shape functions associated with each computational, material point. The Natural Element Method (NEM) does not suffer from the same difficulties as FEM when elements are highly distorted and, although it uses a triangulated mesh, has the advantage that the interpolation func tions are continuously differentiable away from the computational points. NEM shares with DLR the dif ficulty associated with preserving tensor information at integration points after reconnection of the mesh. Braun and Sambridge ( [START_REF] Braun | A numerical method for solving partial differential equations on highly irregular evolving grids[END_REF]) point out that the use of a standard gaussian integration scheme in NEM is not exact for natural neighbor interpolation functions. An other complication is that the shape functions overlap element boundaries, and therefore descriptions of flux boundary conditions is not trivial and may not be ac curate.

Meshless methods

Meshless methods include the Discrete Element Method (DEM, see for example Sakaguchi and � M ühlhaus 1997 [START_REF] Sakaguchi | Mesh free modelling of failure and localization in brittle ma terials[END_REF]), Smoothed Particles Hydro dynamics (SPH, Monaghan 1992 [12]), Element Free Galerkin (EFG, Belytschko et al 1994 [START_REF] Belytschko | Elementfree galerkin methods[END_REF]), and the Point Integration Method (PIM, Liu & Gu, 2001 [11]). In the DEM, computational points are associated with a finite size and shape. Interaction between points (particles) occurs only where they are in contact or close neighbors according to the specified interac tion rules. This is a genuinely discrete method with no homogenization of the particleparticle interactions through a continuum description. The DEM is very well suited to modeling of fracture since particle inter actions can take the form of breakable bonds. The in teraction of a particle with it's neighbours varies with angle and time, since it is a function of the neighbours' coordinates, and the interaction history. SPH and EFG are continuum methods in which the ma terial points are associated with a basis function which can be used as an interpolant for field variables, and to develop an expansion for PDEs. The basis func tions overlap other particles, making it more compli cated to apply constraints to values at individual ma terial points but allowing particles to pass arbitrarily close to one another during simulations of fluid motion with stagnation points. PIM is also a continuum method. In PIM, global basis functions are defined which can be combined to pro duce nonoverlapping shape functions. The method therefore eliminates one of the principal obstacles in efficient formulation of meshless solvers. However, very little in the way of practical applications for the method has been published to date so the promise of the method is largely unexplored. SPH, EFG and PIM are formulated in such a way that implicit solution methods or explicit methods can be used.

Lagrangian integration points

The method used in Ellipsis is a hybrid mesh / parti cle code -the idea being to retain the generality and robustness of meshbased FEM, and capture the geo metrical flexibility of a fullyLagrangian set of particles for tracking material deformation. The computational points are a set of nodes fixed in space, connected by a mesh of finite elements. An independent set of ma terial points which carry material properties and the solution history is embedded in the mesh. Since the material and computational points have been formally separated, a strategy for coupling the two representa tions is needed. The usual finite element interpolation of nodal point values to element interiors is used to update the locations of particles and history variables. The particle properties are coupled to the computa tional mesh through a nonstandard element quadra ture in which the particles which happen to be in a given element are used as integration points. Bound ary conditions constrain unknowns directly and there fore belong with the mesh points. This formulation is very close to the Material Point Method (MPM) developed by Sulsky, Schreyer and their coworkers (e.g. Sulsky & Schreyer, 1996 [18]) who use an Eulerian mesh with Lagrangian particles. Their for mulation was derived from problems where momen tum dominates: impacts between elastoviscoplastic materials, suspensions in fastmoving fluids, and fast granular flows. The strength of the method comes from the incremental way in which particle quantities (in cluding momentum) are updated from the mesh. In compressible creeping flow can be treated by introduc ing a ficticious inertial term, but this is not ideal for the longduration experiments we wish to simulate. The implicit formulation in our approach requires consid erable modification to the method. However, when we introduce elasticity, there are strong parallels with the MPM.

Finite elements with moving integra tion points

In this section we describe the similarities and prin ciple differences between a standard finite element method and one with moving integration points. Prin cipally these differences are in the updating of inte gration points themselves, and in the formulation of a weak form based on particlederived material prop erty matrices. In later sections we discuss efficient grid based solution methods which can be used with our formulation, and some implementation details.

Updating integration points

On one hand, we aim to build a finite element like for mulation then we need integration points to integrate some quantities. On the other hand material particles are used to track history variables and deformation. The essence of our formulation (like in PIC method) is to use particles as integration points. In a particleincell method the grid generally remains fixed whereas particles move through the mesh during computations. Unlike most fixed gridmethods, how ever, the existence of a Lagrangian reference frame al lows the equations to be formulated in a Lagrangian sense (equating the material derivative with the time derivative, see ( 1)) provided all rates of change are com puted particlebyparticle and updated appropriately during the advection phase. This eliminates the diffi culty with numerical diffusion / dispersion of advected quantities.

D ∂ Dt ≡ ∂t (1)
Particle velocities are interpolated from nodal veloci ties and then particle position is updated using a suit able integration scheme such as:

t+Δt x = x t + Δt v n N n (x p ) (2)
p p node � � � �
where v is the nodal velocity and N are the shape func tions associated with the nodes of the element in which the particle currently resides. In practice, a higher or der scheme such as the 4 th order RungeKutta scheme gives a more accurate result. Particle updates can be done in a predictorcorrector fashion, initially updat ing particle locations to obtain velocity solutions and then repeatedly correcting the final locations to obtain a converged velocity. We have found that the improve ments in accuracy from iterating on the particle loca tions are generally too small to justify the increase in complexity, and, for nonlinear rheological laws, the it erative scheme may become unstable. Our preference is therefore to use a 2 nd order RungeKutta scheme also known as the midpoint method and reduce the timestep whenever higher accuracy in the time integra tion is required.

Variational form

The governing equations in domain Ω bounded by Γ are standard for creeping fluids, a conservation equa tion for momentum:

τ ij,j -δ ij p ,j + f i = 0, (3)
subject to a continuity requirement

v i,i = 0, (4)
In the above, v is the velocity, τ is the deviatoric Cauchy stress tensor, p its volumetric component (i.e. pres sure) and f is the specific body force. The notation p ,j denotes differentiation of p with respect to x j . The boundary conditions are given as follows:

τ • n -p = t ¯ on the natural boundary Γ t (5) u = ū on the essential boundary Γ u (6)
in which the superposed bar denotes prescribed boundary values and n is the unit outward normal to domain Ω. Inertial terms are not considered in (3); since the vis cosity associated with creeping flow in rocks is enor mous we can assume the Prandtl number is infinite. This means that the strain rate instantaneously equili brates with the applied boundary conditions and body forces -a factor which strongly influences the way we approach the solution of the equations. Body forces are assumed to arise through density dif ferences in the material. This may be due to composi tional differences between adjacent parcels of material, or temperature variations within the fluid. The latter is not addressed in this paper.

f = (0, -gρ) ( 7
)
where g is the acceleration due to gravity (downwards) and ρ is the material density. We assume the Boussi nesq approximation in only considering the density variations when they contribute to the body force term.

The weak form is obtained in the usual way multiply ing (3) with the test functions N i and (4) with Q. Inte gration over the volume and application of Gauss' the orem yield:

τ ij N i,j dΩ-pQ i,i dΩ = f i N i dΩ+ t ¯ i N i dΓ (8) Ω Ω Ω Γt and � Qv i,i dΩ = 0 (9)
Ω Equation ( 8) must hold for all admissible N i vanishing on boundary Γ u where displacements are prescribed. The notation N i,j indicates the symmetrical form of the derivative: ∂N i /∂x j + ∂N j /∂x i .

Solution strategy

The particle in cell formulation produces a set of ma trix equations in nodal unknowns which have precisely the same form as the standard FEM formulation for the same mesh (though the coefficients in the matrices will be different in most cases). Our solution scheme is very similar to that detailed in Moresi & Solomatov (1995) [START_REF] Moresi | Numerical in vestigation of 2d convection with extremely large viscosity variations[END_REF].

Mixed formulation

We use a mixed formulation in which nodal unknowns are velocity and pressure, and the incompressibility constraint is considered as an independent equation. This gives

� � � � � � K G T G 0 v p = F 0 ( 10
)
where K is the global stiffness matrix, G is the discrete gradient operator. K is built up from particle derived material properties as well as shape function deriva tives, whereas G is purely geometrical and is associated with the mesh. In a mixed method, an arbitrary choice of shape func tions for pressure and velocity can produce unreli able results. The particleincell formulation is no dif ferent in this respect, and we have found it neces sary to combine bilinear velocity shape functions with constant, discontinuous pressures, for example. In other words, despite generally being a lower order ap proximation than gaussian quadrature, the integration scheme based on particles does not underintegrate the pressure field in a way which avoids locking of the elements.

Whereas in many fluid applications or elasticity formu lations it is common to eliminate the pressure equa tions through a penalty method, we have found that this leads to very poor performance of iterative solvers and cannot be used in this case. Instead we follow the � Uzawa scheme in which we apparently eliminate the velocity by writing

G T K -1 Kv + G T K -1 Gp = G T K -1 F (11)
and substituting the expression for G T p to obtain an expression entirely in the pressure:

G T K -1 Gp = G T K -1 F (12)
The procedure outlined in general by Cahouet & Chabard (1988 [8]) and for this specific system by Moresi & Solomatov (1995) [START_REF] Moresi | Numerical in vestigation of 2d convection with extremely large viscosity variations[END_REF] solves this pressure equation through a conjugate gradient series of solu tions for velocities which can be accumulated to give the full solution. Hence the notion that we solve the pressure equation alone is an illusion.

Multigrid solver

The static, regular grid makes it possible to use (un modified) a multigrid solver developed for purely Eu lerian flow problems (Moresi & Solomatov, 1995 [13])

The basis of the method is to build an approximate solution to the problem on a coarse grid by comput ing the stiffness matrix, boundary conditions and force terms at this resolution. The approximate solution is then interpolated to a finer mesh and used as a start ing guess in the iteration. By repeated transfer between coarse and fine meshes, and using an iterative scheme which reduces errors on a length scale defined by the mesh (such as GaussSeidel relaxation), reduction of errors on all lengthscales occurs at the same rate.

If carefully tuned, multigrid methods are capable of solving problems with n unknowns in a time propor tional to n, in contrast to n log n for conjugate gradient method and n 2 or even n 3 for direct solver.

For nonlinear problems, the solution rate is also greatly enhanced by using the coarse levels for numer ous fast iterations to approach the approximate solu tion before interpolating to finer meshes. For example, see Moresi & Solomatov (1998) [START_REF] Moresi | Mantle con vection with a brittle lithosphere: thoughts on the global tectonic styles of the earth and venus[END_REF]. In Ellipsis we use most of the time a full multigrid scheme with a direct solver scheme at the coarsest level.

Numerical implementation

Numerical integration

At the heart of the finite element method is the con cept that the weak form of the equations can be split up into a number of subintegrals over individual ele ments, and that these integrals can be computed ap proximately using a suitable quadrature scheme.

� P e Ω e ψdΩ ≈ w p ψ(x p ) (13)
p=1
where Ω e is an element volume, ψ is a field quantity which is evaluated at a set of sample points with co ordinates x p ; w p are weights for each integration point.

In standard FEM, the locations of the quadrature points, and their weights, are chosen to optimize the integration accuracy for a given set of interpolation functions. The criterion for choosing the quadrature scheme is usually computational efficiency: the mini mum number of locations required to achieve exact in tegration of a specific degree polynomial.

In our formulation, we assume the sample points coin cide with the particles notionally attached to the fluid which therefore move with respect to the mesh. The lo cations of the quadrature points are, as a consequence, given for each element, and it is necessary to vary the weights in order to obtain the correct integral for a given element. The procedure is similar to the deter mination of weights for any quadrature rule: deriving a set of constraints based on the requirement that poly nomials of a certain order have to be integrated exactly, and then equating coefficients to obtain the set of w p values.

We shall consider the one dimensional case in which ψ in ( 13) is a polynomial, and the integral is over the range 1 to 1 (a typical "master element").

α n x n ψ(x) = α 0 + α 1 x + α 2 x (14) 2 + • • •
We integrate [START_REF] Moresi | Mantle con vection with a brittle lithosphere: thoughts on the global tectonic styles of the earth and venus[END_REF] algebraically over the domain, and equate coefficients with the quadrature expansion of the integral from ( 13) to obtain n + 1 constraints on the set of w p values:

nep � w p = 2 (constant terms) p=1 nep � w p x p = 0 (linear terms) p=1 nep � w p x 2 p = 2 3 (quadratic terms) ( 15
)
p=1 nep � w p x 3 p = 0 (cubic terms) p=1
and so on.

In principle, we can imagine choosing a suitable num ber of constraints, and performing an inversion for the w p values, but in practice this would be very time consuming, and may produce very poor results if the particles are unevenly distributed within an element, including negative w p values. If we wish to associate the value of w p attributed to a particle with the rep resentative mass or volume of fluid which it occupies, then w p must be positive. The Material Point Method uses constant weights for the particles to ensure that mass is conserved exactly in the system. Although constant particle weights are not guaranteed to satisfy even the lowest order constraint once the particle's positions evolve to a general con figuration (the numbers of particles per element may vary), the accuracy of the integration scheme is rela tively good for simulations where material strain is not extreme (see benchmark section). With this in mind, we store the values of w p determined from the local volume occupied by the particle in the initial config uration. These reference values are adjusted to best fit the constraints up to a certain degree, and subject to the further constraint that no tracer should have a neg ative weight. For viscous flow using bilinear elements we can obtain optimal convergence rates, in the limit of an infinitely fine mesh, if only the constant and the lin ear constraint terms are used (e.g. , Hughes et al, 1987 [START_REF] Hughes | The Finite Element Method[END_REF]) . To integrate a given order polynomial using particle lo cations, we generally need two to three times the num ber of particles per element (in 2D) than there would usually be Gaussian quadrature points. This is because the locations of the integration points cannot be op timized as they are in Gaussian or other integration schemes but are determined by where the particles happen to be within an element.

� (a) (b) (d) (c) -1.4 Vx Vz 1 0.2 -0.2

Integration scheme

In order to verify the integration scheme for dealing with large deformation problems, we model driven convection in an unit square box (Fig. 2). The default boundary condition is freeslip everywhere except on the top between x = 0.4 and x = 0.6 where a horizon tal velocity is applied towards the right. We track with a darker line the motion of the material along the vertical midline through time (Fig. 2). The velocity solution is uniquely determined by the boundary conditions and is independent of time. The flow carries particles around the box through time, however, which will stir any material interfaces, alter the distribution of particles within elements, and con sequently may disrupt the integration scheme. Intense deformation in the corners results in a need to intro duce new particles resulting in an increasing number of particles as a function of time (Fig. 3). We treat the issue of particle splitting (a local remeshing step) in a subsequent section. We compute velocity errors by comparison with a very fine mesh solution. This solution is obtained at the first timestep for a regular square mesh of 192×192 el ements using a Gaussian integration scheme of 4 inte gration points per element. Our study model is a 48×48 elements with either Gaussian integration points or evenly distributed and weighted points. The error due to the use of a coarser mesh is about 2.8% (Fig. 4.(a) and (c)), the remaining error is due to the integration scheme approximation. The error is computed as follows:

V z -v z 2 V x -v x Err = | � | 2 + | | (16) V 2 + V 2 x z
where V is the reference velocity field obtained with the fine mesh and v is a coarsegrid velocity field. Firstly we compare results between a Gaussian integra tion scheme and a scheme with evenly distributed and weighted particles. The Gauss scheme does not allow for the particles to move, so is only applicable on the first timestep. The particlebased integration scheme was used over several tens of timesteps with two dif ferent numbers of integration points. At each step the particle weights are recomputed to satisfy constant and linear terms in the constraint equations [START_REF] Moresi | The interplay of material and geometric instabilities in large deformations of viscous rock[END_REF]. Although the particle based integration scheme is never better than a fourpoint Gauss scheme, the accuracy is com parable through time, particularly when larger num bers of particles are used (Fig. 4) Recomputing the particle weights is a computationally costly step not required in the Material Point Method which is a close ancestor of the scheme presented here.

In Figure 5, we show the accuracy of integration for the driven convection problem when the weights are fixed, when they satisfy the constant terms of the constraints, and when the constant and linear terms are satisfied. We ran each model for hundreds of timesteps with ini tially 16 particles per element and evaluated the accu racy of the results as well as the fluctuation of the so lution. Fluctuation was computed using a moving win dow of 15 point width to give the average error (Ave(t)) and its variance (V ar(t)). Then we plotted for each case Ave(t) + V ar(t) and Ave(t) -V ar(t). When we keep the particle weight constant so that the particle mass is conserved, an element mass may change when particles cross an element boundary. This simple scheme is associated with large errors with large fluctuations which becomes worse for very large deformation (Fig. 5.(a)) despite an ever increas ing number of particles.

Modifying the weights stabilizes the error through time with a mean value close to that of the Gaussian inte gration scheme. The variation in the error is consider ably smaller when the constraint equation ( 15) for both constant and linear terms is applied (Fig. 5(b,c)). Another reason to include the constraints from (15) up to the linear terms comes from the typical problems we tackle in geological modeling. Invariably there are a strong vertical gravitational forces which are almost exactly balanced by vertical pressure gradients. These must be integrated very carefully to ensure that they do not produce a component of flow since even a small error contribution from these hydrostatic terms would easily dominate the solution.

The Gaussian integration is the most efficient of the schemes we have implemented, but is not suited to large deformation particleincell methods. We obtain reasonable results using four particles (initially) per el ement and adjusting their weights as they are moved with the flow. Sixteen particles are needed before we have enough degrees of freedom to achieve the same integration as a four point Gauss scheme in 2D, allow ing for the fact that some particles will be in the wrong position within an element to contribute to the inte grals.

Element matrices and particle properties

In a standard finite element formulation, the stiffness matrix for an element, k e is built up in a segregated form, often written as:

e k e = B T DB b dΩ (17) ab a Ω e
where B a is a matrix comprising shape function derivatives (associated with node a) obtained from the constitutive relationship, and D is a matrix of material properties. For our formulation, the material property matrix is considered to be an attribute of an individual particle together with its current state and history and not a property of the mesh. Therefore, in the context of the numerical integration scheme discussed above,

nep k e = w p B T a (x p )D p B b (x p ) (18)
ab p=1

where n ep is the number of particles which happen to be in an element.

Particle splitting

Fluid flow close to a stagnation point produces an elon gation in one direction and a corresponding shorten ing in the perpendicular direction. This has the effect of distorting the original volume local to each tracer into a narrow filament (Figure 6 (a) and (b)). This may mean that the fluid initially associated with a particle lies in several different elements, while the integration scheme simply lumps the entire amount into the ele ment containing the tracer itself. In the worst case this may leave some elements entirely empty of tracers and others considerably over represented.

� � � � � � � � �
The remedy is to ensure that the volume of fluid asso ciated with a particle never becomes too distorted. We keep track of a local measure of strain associated with each particle and use this to generate new particles nearby which occupy the extremities of the distorted local volume. This is illustrated in Figure 6(b) where the heavily shaded particles sitting at the centroid of the salamishaped local regions represent the original occupants of the volumes from Figure 6(a). The ma terial within the element is poorly represented by the particles which actually contribute to the element in tegrals. The lightly shaded particles are later additions which aim to correct this problem. Local volumes cor responding to the new particle distributions are indi cated in Figure 6(c).

When splitting particles, we attribute the same history variables to the copies as were held on the original par ticle. To ensure that this approximation is a reasonable one the particle splitting should occur when the distor tion is relatively small.

Element inverse mapping

It is usual to change variables in the element integrals such as [START_REF] Sulsky | A particle method for historydependent materials[END_REF] to a regular master element. This greatly simplifies the computation, but in our case there is an additional step required before we can map to a master element. The particle positions are known in the global coordinate system and must first be mapped into the element local coordinate system. Zhao et al (1999) [START_REF] Zhao | A consistent pointsearching algorithm for solu tion interpolation in unstructured meshes con sisting of 4node bilinear quadrilateral elements[END_REF] found an algebraic mapping for bilinear quadrilateral elements, but a more general approach is required for higher order elements and for 3D.

The notation is given in Figure 7 for the 2D case -the extension to three dimensions is straightforward. e ξ and e η are unit vectors in the 'natural directions' of the distorted element which map to the ξ and η axes re spectively in the master element. h ξ and h η are char acteristic dimensions of the element in the appropriate directions. We wish to map the coordinates, x p of the particle p to the coordinates in the master element, ξ p . The procedure takes the form of a predictorcorrector iteration: we first guess an initial value of ξ p and use 0 this to predict the global coordinates x p ,

ξ p = (0, 0) nen nen � � ( 19
) 0 x = N n (ξ p )x n , N n (ξ p )y n p n=1 n=1
where n en is the number of nodes in the element, and x n are their coordinates. We compute ξ p through a number of corrector steps:

i ξ p ← ξ p + β e ξx x p + e ξy x p i /h ξ i η p ← η p + β e ηx x p + e ηy x p i /h η � � ( 20
) nen nen i+1 = x p N n (ξ p )x n , N n (ξ p )y n n=1 n=1
i is the iteration index, β is a relaxation parameter which we chose to be 1.0 for the first iteration and 0.9 for subsequent iterations. When e ξ and e η are orthog onal, the iteration completes in one step, otherwise it is necessary to repeat the correction step until the pre i dicted x p is within a satisfactory tolerance of the known value.

We use ξ within the element search algorithm since the master element coordinates give an immediate de termination of whether the particle lies inside or out side the element even for highly distorted elements (if the element is too distorted for this procedure to work, then it is also too distorted to use to create the element matrices)

Moving boundary conditions

One major difficulty with the Eulerian mesh is in the application of velocity boundary conditions normal to the surface as the mesh does not automatically track the boundary position. If the sense of the boundary ve locity is into the domain, then the method requires the generation of new particles to represent the incoming material. If the sense is out of the domain then parti cles, and their stored history, will be lost.

To avoid such difficulties, we recommend updating the mesh in between solution steps to follow any velocity boundary conditions normal to boundary surfaces. As pointed out by [START_REF] Sulsky | A particle method for historydependent materials[END_REF] [START_REF] Sulsky | A particle method for historydependent materials[END_REF], the formula tion is unaffected by changes to the mesh occurring af ter the particle positions have been updated since the mesh carries no information other than domain geom etry and boundary conditions: material configuration and history information is entirely determined by the particles.

In general the mesh can be regenerated completely to have the optimal number of nodes and configuration for the updated geometry. In the compression and ex tension examples given here, however, the mesh topol ogy remains unchanged -the mesh is simply scaled in the appropriate direction to minimize the recompu tation of the connectivity matrices. It is important to use the same updating procedure for nodal points as particles to avoid the boundary nodes overtaking the particles, for example.

Constitutive Relationships

Viscoelastic formulation

We use a Maxwell viscoelastic constitutive relationship which assumes that the deformation rate is the sum of viscous and elastic parts

τ 2µ + τ 2η = Dv + De = D (21)
� where τ is the Jaumann corotational stress rate for an element of the continuum, µ is the shear modulus and

� � � � � � � � �
η is shear viscosity. D is the deviatoric part of D. which we solve by substituting for τ t+Δt e to give a set � τ of equations for velocity unknowns.

= τ ˙+ τ W -Wτ (22)
In choosing a material timescale (Δt e ) independent of the numerical advection timestep (Δt), it is necessary where W is the material spin tensor, for a given particle to carry with it a stress history of

-W t τ t + τ t W t (24)
Yielding t 1 ∂V i ∂V j several advection timesteps (corresponding to several (23) W ij = complete solutions for the stress field). This is achieved 2 ∂x j ∂x i through an averaging procedure in which the stress tensor stored on a given particle is averaged with the newly calculated stress tensor (τ t ) at the same location:

τ t+Δt Δt Δt τ t + τ t+Δt e (30)
e e

ˆ= 1 -Δt Δt

This stress is then advected and rotated with the parti cle to give the updated stored stresses τ t+Δt .

The W terms account for material spin during advec tion which reorients the elastic storedstress tensor.

As we are primarily interested in solutions where very large deformations may occur -such as buoyancy driven fluid convection, we prefer to work with a fluid like system of equations from the outset.

Hence we obtain a stress / strainrate relation from ( 21) by expressing the Jaumann stressrate in a difference form: where α = η/µ is the shear relaxation time.

τ t+Δt e -
We can simplify the above equations by defining an ef fective viscosity η eff :

e Δt η eff = η (26) Δt e + α
Then the deviatoric stress is given by

τ t+Δt e = η eff D + e + 2 ˆt+Δt e τ t W t τ t -τ t W t (27) µΔt µ
Our system of equations is thus composed of a quasi viscous part with modified material parameters and a righthandside term depending on values from the previous timestep. This approach minimizes the modi fication to the viscous flow code. Instead of using phys ical parameters for viscosity we use an effective value (26) to take into account elasticity, then, during com putations for the force term, we add elastic internal stresses from the previous timestep or from initial con ditions.

F e,t = η eff τ t (28) i e ij,j
-µΔt Therefore (3) becomes

τ t+Δt e -δ ij p ,i + f i + F e,t = 0, ( 29
) ij,j i
In a geological context we frequently deal with situa tions where part of the system is subjected to very high stresses. Under such conditions the material fails, but, unlike many practical engineering simulations, we are interested in simulating in the postfailure behavior up to very large strains. On the basis that the post yield deformation trends in creasingly to dominant, simple structures with increas ing strain (BenZion and Sammis, 2001 [START_REF] Benzion | Characterization of fault zones[END_REF]), our geo logical modeling at scales of tens to thousands of kilo meters uses very simple descriptions of yielding. Brit tle behavior is parameterized using a nonlinear effec tive viscosity which is introduced whenever the stress would otherwise exceed the yield value τ yield . This ap proach ignores details of individual faults, and treats only the influence of fault systems on the largescale convective flow.

To determine the value of the effective viscosity at any point we extend (21) by introducing a PrandtlReuss flow rule for the plastic part of the stretching:

τ τ τ 2µ + 2η + λ 2 |τ | = D e + D ˆv + D ˆp = D ( 31
)
where λ is a parameter to be determined such that the stress remains on the yield surface, and 2) . The plastic flow rule introduces a non linearity into the constitutive law which, in general, re quires iteration to determine the equilibrium state. The implementation is as follows, starting from equa tion (31), we again express the Jaumann stress rate in first order difference form (using the Lagrangian parti cle reference frame):

|τ | ≡ (τ ij τ ij /2) (1/
τ
We determine λ by equating the value of |τ t+Δt e | with the yield stress in (33). Alternatively, in this particular case, we can obtain η � directly as ing and separation of the two parts of the sample oc curred with barely any movement of the end boundary.

Even with a material which has no strain softening, there is a tendency for deformation to localize in a particleincell representation of the sample. This oc curs because the sample boundary is never perfectly flat (as in real life) due to numerical fluctuations in the particle locations, and to a mild interference (moiré) effect between the array of particles and the underly ing grid. These effects produce small fluctuations in the stress field which can result in early failure at cer

η = τ yield / Deff (35)
where tain points. Once nucleated, shear bands can propa gate from these points -ultimately resulting in neck ing and complete separation of the two halves of the sample.

In Figure 10a, we plot the stress at each of the sam ple points in the material as a function of time for a fixed end velocity. The evolution of stress within the beam was close to linear apart from the influence of the changing of the beam thickness during deforma tion. The yield stress of the material was 3 × 10 5 . The stress increased within the sample at the same rate for all the sample points until the yield stress was reached. At this stage, the stress was not able to increase any fur ther, and the material deformed uniformly at the yield stress. Once localization had occurred, however, points outside the necking area begin to unload, and the stress dropped dramatically. The rate at which stress drops from yield back to zero is governed by the viscous part of the rheology, and the presence of a low viscosity background material.

The unloading is more clearly seen in the plots of the displacement of the sample points through time in Fig ure 10b. Before yielding, the displacement of each sample point increased monotonically. Once yielding occured, and the deformation localized, the sample points on the left of the break (a,b), under the action of stored elastic stresses, rapidly retreated towards their original locations. The sample point on the right of the break (c) moved rapidly to the right as the elastic defor mation relaxed.

It is worth discussing at this point a consequence of the fact that the yield criterion only applies to the devia toric stress. During the separation of the layer, the pres sure becomes enormous at the constriction, which ob viously could not occur in a real material. To model this situation in a more realistic manner we would need to complement the yield criterion on the deviatoric stress with a suitable tension cutoff condition.

ˆD µΔt e τ t + 1 (W t τ t -τ t W t ) (36) D eff = 2 ˆt+Δt e + 1 µ and |D = (2D ij D ij ) 1/2 .
| The value of λ or η � is iterated to allow stress to re distribute from points which become unloaded. The iteration is repeated until the velocity solution is un changed to within the error tolerance required for the solution as a whole.

Applications Extension of a layer

The yielding algorithm is benchmarked by measuring the second invariant of the stress and displacement at points within a viscoelastic beam which was extended or compressed at a fixed rate, v = 5, by an imposed velocity boundary condition at one end. Figure 8 in dicates the geometry of the numerical experiment: the mesh was initially 3 units long by 1 unit high. The sam ple was 0.5 units thick, occupying the central half of the mesh, and was surrounded by a low viscosity, com pressible material. Three sampling points (a,b,c) for recording the stress invariant and displacement were chosen within the sample initially placed along the midline at x = 0.2, 0.5, 0.8. The material parameters of the sample (η = 10 8 ,µ = 10 6 ) were chosen such that the relaxation time was long (α = η/µ = 100 ) compared to the duration of the ex periment ( 0.25 ) so that the material behaved almost as an elastic solid (high Deborah number, De, defined as relaxation time / observation time). Figure 9 shows the progress of the experiment. Ini tially, deformation was uniform, resulting in gradual stretching of the sample (t ≤ 0.180). The entire sample reached the yield point at the same time (t = 0.212) and initially deformed uniformly with all points yielding. However, the deformation soon localized to a number of shear bands (t = 0.220), then to two places along the

Folding of rock strata

Geologists are intrigued by the formation of folds in sedimentary and metamorphic rocks on a large range of scales. Different conditions are necessary to create such buckling but it is often associated with tectonic compression. Other parameters are rock mineralogy, temperature, pressure, water content etc. We are pre dominantly interested in studying buckling as a func tion of elastic and viscous properties of the different layers. As a first step to show the capability of this new formulation to handle large deformation, we present in this analysis a folding problem of a competent layer embedded into two weaker purely viscous semi infi nite layers. The initial geometry (Fig. 11a) is perfectly linear which means that there is no trigger for the instability to start except numerical noise. An unit inward velocity boundary condition is set up on the right side of the mesh. All the others boundaries are freeslip. We model the strong layer as viscoelastic with a yield criteria to follow qualitatively the evolution of the folding and the failure. Figure 11b shows the evolution during the ho mogeneous shortening phase where numerical noise triggers the instability and internal elastic stresses are built up. Once the instability is triggered the folding of the competent layer starts (Fig. 11c) and show a wave length which is in accordance to the Biot theory (1965) [START_REF] Biot | Mechanics of incremental deformations[END_REF]. If we keep compressing then bending stresses are built up in the beam with the largest tensile stress in the outer part where the curvature is the highest. Then due to tensile failure the beam breaks into a number of pieces (Fig. 11d) and the internal elastic stresses are released which straighten the layer. The final numeri cal plot is to be compared with what geologists can ob serve in situ (Fig. 11e). Even if the reason of such a mi crostructure is not perfectly understood (brittle failure, perpendicular shearing, anisotropy ... etc) the model presented in this paper looks reasonably accurate.

Thermal convection with suspended ellipti cal crystals

In addition to the set of equations described so far, an energy equation can be solved explicitly in conjunction with the timemarching scheme:

∂T ∂T ∂ 2 T + v i ∂x i = κ ∂x 2 j , (37)
∂t where x i are the coordinates, v i is the velocity, T the temperature and κ is the thermal diffusivity. This thermomechanical formulation is applied to the convection solution used as a benchmark by Blanken bach and coworkers in [ [START_REF] Blankenbach | A benchmark comparison for mantle convection codes[END_REF]] The convection solution in the absence of suspended particles is uniquely described by a single dimension less number: the Rayleigh number, Ra = gραΔT /ηκ, where g is the gravitational acceleration, ρ is the refer ence density, ΔT is the temperature contrast across the layer, and η is viscosity. The density contrast in the fluid due to temperature variations is (αρΔT ), and the ratio of density between fluid and solid at equivalent tem perature is denoted by β.

We suspended a set of solid crystals with β = 10 in a steadystate convecting viscous fluid with Ra = 10 5 in an initially regular pattern (see Fig. 12 (a)), and let them settle (Fig. 12 (b)). When the crystals had sunk to the bottom of the layer (Fig. 12 (c)), the convective circula tion was mainly confined to the upper part of the layer where there were no crystals. Within the lower part of the layer, fluid percolated between the densely packed solid crystals.

Once particles had settled (see Fig. 12 (c)) β was ar bitrarily changed to 0.5 which immediately reversed the density contrast between fluid and crystals and re sulted in an fullscale overturn of the whole layer (see Fig. 12 (d) and (e)). In the end particles reached a new equilibrium ( Fig. 12 (f )) on top of the convecting fluid. This simple demonstration shows the promise of the method in simulating the dynamics of crystalrich magmas, potentially accounting for evolving compo sition which would alter the density contrast between liquid and solid during the course of the simulation

Conclusions

We have presented a Lagrangian integration point fi nite element formulation designed to handle large de formation for viscoelastic materials. The scheme de rives from the Material Point Method but differs in a number of important aspects, including the fact that it is based on a fastimplicit solution method, and that it includes various particle reweighting steps which im prove accuracy in the fluiddeformation limit.

We have demonstrated that the particleincell finite element scheme is comparable in accuracy with tra ditional Eulerian finite element methods for fluid dy namical problems.

The principal advantages over other methods are: The ability to continue to extremely large deformation without significant change in accuracy; The ability to track material history and interfaces through time with accuracy comparable to more specialized methods; the fact that a regular grid is preserved throughout the run which allows fast numerical solvers to be employed.

Disadvantages include the fact that the hybrid mesh/particle approach does not allow particles to communicate directly, and therefore lacks some of the flexibility of pureparticle codes. Secondly, the resolution is related to the grid point spacing, not the finer particle spacing. Finally, while it is conceptu ally simple to retain the mesh, and computationally convenient, the storage requirement is very high.

Figure 2 :

 2 Figure 2: Geometry of the convection model at initial timestep (a) with velocity profiles, timestep = 150 (b), 500 (c) and 800 (d).

Figure 3 :

 3 Figure 3: Total number of particles versus timestep for the driven cavity benchmark.

Figure 4 :

 4 Figure 4: Error versus timestep for 4 Gaussian integra tion points (a), and for 4 (b) and 16 (c) particles ini tially evenly distributed within the element with equal weights.

Figure 5 :

 5 Figure 5: Error versus time for different particle weight ing schemes: (a) particle weights are conserved exactly through time (except for particle splitting), (b) parti cle weights are adjusted to fit the constant terms of the constraints, (c) particle weights are adjusted to fit the constant and the linear terms of the constraints.

Figure 6 :

 6 Figure 6: Integration schemes become more compli cated when large material strains produce elongated "local volumes" for particles. (a) Initial configuration, (b) flow near a stagnation point elongates domains: solid circles are original particle points, open circles in dicate new locations within the original domain which have similar spacing to the original interparticle spac ing, (c) remapping the domains to the solid and open circles again allows the particles to contribute to the in tegrals of the correct element. The fact that this proce dure produces overlapping domains is an indication of the approximate nature of the integration scheme.

Figure 7 :

 7 Figure 7: Coordinate systems in global mesh, distorted element, and master element reference frames

Figure 8 :

 8 Figure 8: Geometry for simulation of the extension of a viscoelastic bar with a yield stress.

Figure 9 :

 9 Figure 9: Simulation of the extension of a viscoelastic bar with yield stress. Black shading indicates regions deforming at yield. Embedded marker points which follow the material deformation are indicated by a,b,c.

Figure 12 :

 12 Figure 12: Numerical simulation of the settling of dense crystals and the uplift of light crystals in a convecting Figure 10: Stress, and displacement at sample points viscous fluid. (a) Initial setup, (b) settled crystals, (c) a,b,c as a function of time for the extension experiment time at which density is decreased, (d) and (e) light of figure 9 crystals rising, (f ) steadystate

Figure 11 :

 11 Figure 11: (a) Initial configuration of the layers, (b) homogeneous shortening, (c) viscouslike folding, (d) failure of the layer and (e) Quartzofeldspathic lay ers (light colors) defining asymmetric folds in Archean migmatitic gneiss, Simo, Northern Finland.

  No modification to the isotropic part of the problem is sample (t = 0.227), and finally to a single location (t = required when the von Mises yield criterion is used. At 0.2359) which focussed all subsequent deformation un yield we use the fact that |τ | = τ yield to write til the sample failed entirely (t ≥ 0.2404). The frames are not uniformly spaced in time since the postfailure � given by ητ yield µΔt e η � = ητ yield + τ yield µΔt e + ληµΔt e

	τ t+Δt e	� = η � 2 D	t+Δt e	+	1 µΔt e	� τ t + (W t τ t -τ t W t ) 1 µ (33)	behaviour occurred on a much shorter timescale than the gradual loading. Note, for example, that the neck
	using an effective viscosity, η
			�			� � �	� � �
								t+Δt e	1 2µΔt e	+	1 2η 2 |τ | λ + + =
								ˆt+Δt e D	+	1 2µΔt e τ t +	1 2µ	(W t τ t -τ t W t ) (32)