Chris Beaumont and Philippe Fullsack
In the summary of our proposal we outlined the planned development of finite element numerical computational techniques designed to allow various aspects of the geodynamics of the lithosphere to be investigated through numerical model experiments. We regard the problems to be addressed as important general problems that have been identified partly by Lithoprobe research but, as yet, have no clear answers. As such, our aim is to contribute to the Pan-Lithoprobe synthesis through synthetic investigations that focus on mantle-lithosphere interactions, particularly during plate convergence, subduction, and ensuing continent-continent collision.
Our earlier numerical model research has focussed on mechanical deformation and thermal consequences of orogeny at the crustal scale. In this research it was assumed that mantle lithosphere beneath the crust subducted in a manner that is similar to the subduction of oceanic plates. Although increasingly improved geological and geophysical observations provide a general picture of the crust, the fate of lithospheric mantle during continental convergence remains enigmatic.
One purpose of the proposed numerical developments is to allow efficient computation of mantle-lithosphere flows at the scale of the upper mantle with sufficient dynamics that the model is no longer under strong kinematic control. In addition to the large-scale and associated large convergence and large deformation, the calculations must have sufficient resolution to compute accurately effects such as strain hardening/softening, which appear to be critical processes in subduction. One method to address these requirements is to develop adaptive finite element techniques such that the finite element mesh can be automatically refined where high spatial resolution is required. These techniques are computationally intensive and require parallel processing to be efficient.
In this report we first discuss progress on upper mantle scale model experiments using a single processor finite element code (SOPALE). The results provide some preliminary insights into the fate of continental mantle lithosphere under a range of circumstances, such as different boundary conditions, material properties and distribution, and a range of temperature regimes. We then describe the limitations of the current computations. Lastly, progress in the development of improved numerical techniques, as discussed in the proposal, is described. Although these developments are described in the context of lithosphere-mantle scale flows, the techniques are much more general and are designed to address the thermal mechanical deformation of quite general elastic-viscous-plastic materials during creeping flows.
Given that the high resolution parallel computation techniques are currently under development, we will illustrate the types of problems that will be addressed by showing results obtained using a relatively coarse resolution calculation and with no adaptive/ refinement of the finite element mesh. The calculations were made using the program SOPALE. The models are examples of interactions between simple plastic/viscous lithospheres and an underlying uniform linear low viscosity upper mantle. We use 'plasticity' in the mechanical/engineering sense to imply a material with a finite instantaneous yield strength. This usage is different from the microrheological 'crystal plasticity', which may be an effective viscous creep bulk rheology. These models were provided by Chris Beaumont, Russ Pysklywec and Ritske Huismans.
The model configuration (figure 1) shows an early stage in a typical calculation; in this case model DP13. Actual results are shown for the models DP14 (figure 2) and DP15 (figure 3). The purpose is to investigate lithospheric deformation as a consequence of far-field convergence specified by the velocity boundary condition at the right hand boundary. A small weak region embedded in the strong upper mantle lithosphere acts to seed local deformation.
In both DP14 and DP15 the lithosphere has a 2-layer crust comprising a frictional plastic upper layer and a lower layer with the same properties (DP15) or that is weak and viscous (DP14). The uppermost mantle is also frictional plastic with the same internal angle of friction as the crust. The lower mantle lithosphere is uniform linear viscous and is twice as thick in DP14 as DP15. The density field changes at the surface, the base of the crust, and at the base of the lithosphere (figures 2 and 3).
The deformation, shown at 100, 200 and 300km convergence, is a consequence of the superimposed effects of the imposed horizontal contraction and density driven flow. In particular, because the lower lithosphere is more dense than the asthenosphere, it is unstable and exhibits a Rayleigh-Taylor instability causing it to 'drip'. This is the type of behaviour that Greg Houseman, Peter Molnar and colleagues have advocated for the convective removal of thickened lower mantle lithosphere beneath continental orogens. They have, however, only investigated this instability for purely viscous models of the lithosphere. In the current models we see a mixed mode of mantle lithosphere deformation. The strong (plastic) upper mantle lithosphere 'subducts' (similarly to oceanic lithosphere), while the lower mantle lithosphere drips. There is a strong positive feedback from the lithospheric thickening, caused by the horizontal contraction, to the finite growth of the Rayleigh-Taylor instability. Of course, this may not represent Earth-like behaviour. However, there is some indication of symmetric thickening of the lower mantle lithosphere (like the 100km convergence results) beneath the Southern Alps orogen New Zealand, and a sheet-like body beneath the Transverse Ranges, California that may be a drip (intermediate to the 100 and 200km stages of DP14).
The DP14 result also shows that in this case the crust is subject to a basal boundary condition that is close to 'subduction'. That is, the subduction assumption used in our earlier modelling can be correct provided that the uppermost mantle is strong and subducts. This boundary condition will act on the crust even though the lower mantle lithosphere may be weak, viscous and gravitationally unstable.
Comparison of DP14 and DP15 also shows that limited subduction/underthrusting of the entire crust is possible when the lower crust is strong and couples the entire crust to the uppermost mantle (DP15). In contrast, weak lower crust favours decoupling and there is no crustal subduction in DP14. Again this may not be Earth-like, but it indicates that the rheological layering of the lithosphere may be as important a control on what subducts as the density structure.
Figure 4 illustrates the evolution of a similar model except in this case the lithosphere has three layers, Von Mises plastic crust and uppermost mantle, and a uniform linear viscous lower mantle lithosphere. Convergence is again specified at the right side boundary. The density structure is the same as the models described above. The initial style shows a regional thickening but this evolves to give subduction of the plastic mantle lithosphere beneath the decoupled crust. In this regard the results are similar to crustal scale models subject to subduction-style basal boundary conditions. Later, the viscous lower mantle lithosphere becomes unstable and drips, followed by slab breakoff of the subducted plastic upper mantle lithosphere.
This evolution is one of a range of behaviours discussed by Pysklywec et al. (Geology, in press). The importance of these models is that they exhibit a richer range of deformation styles than purely viscous models, for example subduction in addition to symmetric thickening. Two key factors that control these differences are the finite strength plastic uppermost mantle and the strain softening that occurs in these plastic materials (shown graphically in figure 4, for example). The results demonstrate the critical need for improved knowledge concerning lithospheric bulk rheology. Under what circumstances do finite strength layers exist in the lithosphere and how does mantle lithosphere strain soften or harden?
The last example (figure 5) illutrates behaviour of a back-arc region in the upper plate of a subducting oceanic lithosphere. This model is similar in design to some of Chemenda's laboratory models and has the same scaled material properties for the upper part of the lithosphere with an additional unstable viscous lower mantle lithosphere. The model evolution illustrates the coupling of two Rayleigh-Taylor instabilities; the downward dripping of dense lower mantle lithosphere, as in the earlier models, plus the buoyant intrusion of asthenosphere into the thinned upper plate. The coupling between these two flows produces convective cells that cause the upper plate adjacent to the subduction zone to rift and migrate toward the subduction zone under some combination of 'ridge push', 'subduction suction', and 'basal drag' forces. This motion at 6.5Ma occurs even though active convergence and subduction of the oceanic lithosphere ceased at 1.6Ma. The time lag is a measure of the finite growth time of the two gravitational instabilities. The model may represent a prototype for the recent evolution of the Tyrrhenian Sea (figure 6). These results are preliminary and the model properties are not necessarily close approximations to the Earth. Nevertheless, they do represent dynamical experiments of lithosphere-mantle scale interactions in which the models have sufficient resolution that crustal and mantle lithosphere layers can be distinguished and in which plastic, as well as viscous, rheologies exist. Additional properties of the SOPALE models are that the upper surface is a true free surface and the full density and pressure fields are calculated. The latter is important for frictional plasticity.
It is first necessary to state that we are very pleased to be able to report on the success of the current (SOPALE) lithosphere-mantle scale models... the 'Deep model' that we have focussed on developing for some time. It is also important to note that SOPALE is a fully-coupled thermal-mechanical code that includes a wide range of rheological models (grain size sensitive creep, thermally activated power law creep, Peierl's law creep), in addition to the more basic rheologies (linear viscous, Von Mises plasticity and frictional plasticity (Drucker-Prager)) used in the examples outlined in section 2.
There are many interesting and potentially important problems that can be addressed through model experiments using SOPALE and only very simple ones have been completed to date.
However, it must be recognized that SOPALE is limited in that the ALE calculations are carried out on structured grids with a fixed topology. If the size of the FE grid is increased, the time required for each solve of the corresponding matrix equation, even with efficient direct solvers, increases rapidly such that serial solution on one processor is no longer efficient. The computational requirements limit SOPALE calculations to Eulerian meshes of approximately 100 by 200 for models with 1000's of timesteps.
A limitation is therefore one of spatial resolution such that boundary layers, regions of high gradients in model fields and material properties cannot be accurately resolved and tracked during the dynamical evolution of model experiments. This, together with other limitations, necessitates the development of computational techniques using adaptive/self refining meshes on parallel systems in order to achieve accurate calculations. These developments are the core of our proposal.
Most of the capacity of today's computers for processing larger amounts of information lies in their parallel architecture. One has therefore to rewrite or rethink algorithms such that data and operations are shared more or less equally (so that waiting is minimized) across an array of processing units concurrently treating the information. Communication between these units, needed to coordinate tasks, also has to be minimized. Data structures need to be mapped (i.e. cut into pieces) on the array of processors so that memory requirements are shared and no processor is saturated with excessively large data arrays. All aspects of information processing, namely flow computation and pre/post-processing, need to be parallelized so that bottlenecks are avoided in the processing chain.
The following sections present in more detail the algorithmic components that we have chosen in our FE environment Mozart to achieve these objectives to some degree. These sections should be related to the corresponding parts of the proposal.
Scalable algorithms offer prospects of higher spatial resolutions (e.g. 2D 1000x400 and 800x800 experiments have been run, approximately 1 km resolution at the upper mantle scale). In particular, deep models, equivalent to those described earlier, can achieve a much higher crustal resolution than in SOPALE. Local high gradients (whether originating in loss of stiffness, material contrast or boundary layers) may be adaptively followed. Similarly, localization in strain or temperature dependent viscoplasticity (shear heating, strain softening/hardening, grain size reduction/growth) may be followed. Refinement may allow convergence of the ALE method to a quasi-Lagrangian method. In particular, large deformation may be followed for elasto-viscoplastic materials.
Practically all aspects of Mozart , ranging from code management to parallel Finite-Element applications have been modified by new developments. Our major goals were the clarification of existing data structures and the creation of new data structures/operations allowing for both structured and unstructured (serial or parallel) FE computations. The relationship of these tools to the overall structure of Mozart is shown in Figure 7.
We quickly review these developments here.
We have adopted a strategy of mesh refinement to enable mesh adaptivity. In order to circumvent the drawbacks of both Delaunay meshers (e.g. TRIANGLE, discussed in the proposal) and our previous longest-edge-bisection routine, we have created a new mesh refinement tool named "HALFNER" after Mozart's Haffner serenade. HALFNER uses topological 'origami-like' rules to refine grids. All components of our algorithm generalize to meshes of arbitrary dimension and are therefore applicable to 3D. HALFNER has been used in many tests involving Eulerian and/or Lagarangian mesh adaptivity. To give a few examples, computations have been performed that develop dense mesh regions near a free surface boundary layer, in zones of high deformation (or rate of deformation), in the vicinity of material interfaces, etc.
Tools dealing with the representation of graphs mapped from structured or unstructured meshes have been further developed (GRAPH) in association with structures representing linear systems (LSYSTEM).
A randomized algorithm based on conformal mapping onto a sphere, along the lines proposed in  has been implemented to optimize load balancing and minimize cross-edge cuts of the mapping of arbitrary meshes onto an array of processors (GEOMPART). As described in , the algorithm randomly samples the mesh, uses radon reduction to find approximate counterpoints, from which a reference stereographic projection is computed. Randomly chosen equators of the sphere created by this process are expected to provide optimal vertex separators for meshes satisfying mildly restrictive regularity conditions.
Our prototype parallel structured FE code OPALE (see proposal) has been improved into a new code: WOLFGANG.
By contrast with OPALE, WOLFGANG is based on new Mozart primitives. It allows Lagrangian refinement, uses MPI-gathering routines for Lagrangian parallel advection (routine master_search), allows user-definable primitives for boundary conditions and constitutive equations (USERSPACE: rheologies, densitologies, thermologies are Mozart acronyms corresponding to functions defining the mechanical, mass and thermal physical models to be used in evolution equations), and uses a new interface (LSYSTEM) with PWSSMP. This last library was provided by A. Gupta  (IBM Watson Research Center) and is optimized for use on Dalhousie's new SP supercomputer whose architecture combines distributed and shared memory (the programming model is a hybrid of message-passing and multi-threading). Master- and peer- mode parallel plotting and IO routines are also available in WOLFGANG (but were not available in OPALE).
Preliminary runs of thermal convection  and sandbox experiments have been carried out with WOLFGANG on up to 24 processors of the Dalhousie SP computer.
CUBE is our structured 3D incompressible flow solver rewritten in April 2000 in order to create a new set of Mozart primitives. This tool, still in progress, is our testing ground for multi-scale, highly structured computations and iterative Krylov methods. New matrix assemblers are being designed in CUBE to form completely disconnected (atomic) or partially disconnected (molecular) systems. CUBE combines vector and parallel features. A stabilized context [4,5] UZAWA algorithm has been written to accelerate convergence. This algorithm is evolving towards an abstract tool in which different matrix types will result in a variety of acceleration methods.
Low level finite element primitives have been written and used to reconstruct (from the same functions) the Q4/1 , Q9/3, and CR elements. The Pian-Sumihara 4-noded element has been implemented , because it is regarded as the most accurate FE quadrilateral and both works in the incompressible limit and accurately represents bending modes. Anisotropic elements have been introduced and used.
Significant efforts have been made to clarify developments.
Experiments are listed in arbitrary order and are just indicators of Mozart's potential.
Experiment 1: (Bern, July 1999) Successful comparison of 2-layer viscous folding numerical and analytical solutions.
Experiment 2: Adaptive advection of 1D-Lagrangian closed curves in a steady convection flow shows the capability of following fine-scale mixing.
Experiment 3: Adaptive CR computation of visco-plastic flows in the presence of gravity and singular (Heavyside function) kinematical boundary conditions.
Experiment 4: Successful comparisons of solutions obtained from various elements.
Experiment 5: Successful comparisons of thermal Boussinesq steady convection flows with published numerical experiments (cases 1 and 2 in  and SOPALE). Similar tests were performed in parallel with WOLFGANG.
Experiment 6: Lagrangian grids are equipped with a 'microstructure', which takes the form of a material vector polarizing the rheology. The resulting anisotropic flow is then followed and it advects the anisotropy direction and updates the anisotropy contrast. This serves as an example illustrative of the behaviour of foliated material. In this experiment, rheologies characterized by foliation-parallel and foliation-perpendicular viscosities were used. The ratio of these viscosities may evolve as a function of other variables such as strain, thus showing combination of anisotropy development and strain softening/hardening. Nonlinear plastic generalizations could be most interesting.
Experiment 7: Tests of adaptive Lagrangian /Eulerian refinement. A mode of deformation is chosen in a series of models (analytical stream function, singular modes corresponding to 0 or infinite rates of deformation). Lagrangian grids are advected and adapted to deformation, Eulerian grids are created to conform to the Lagrangian node distribution and these grids are geometrically partitioned onto an array of processors.
Experiment 8: Pure shear force driven compression Lagrangian tests were successfully compared to analytic solutions for elasto-viscous Maxwell rheologies.
A word of caution is in order here, to give a proper measure of the difficulties to be expected sooner or later. Higher spatial resolutions imply, via the Courant-Friederichs-Lax convergence condition, the requirement for higher temporal resolution, and therefore a higher number of operations to get to a desired deformation level. Adaptive computations are subjected to many simultaneous refinement criteria and the solution must be proven to be robust, e.g. in the sense of being mesh-independent. Nonlinear convergence may be slowing down as the grid is refined (forever in the case of singular problems); this would indicate that alternatives to global matrix reconstruction should be found. Refinement algorithms coded in Mozart only allow pseudo-coarsening in the sense of using a mesh at a lower level in the nested grid hierarchy. This of course is not general enough. CR FE solvers are more memory and computation intensive than simpler (e.g. Q4/1) computations, a drawback that we hope will be more than compensated by the use of mesh-adaptivity. Sparse parallel direct solvers are likely to be unable to provide, alone, acceptable computation times for 3D computations (we confirmed this view by some symbolic tests on the SP showing that grids larger than 50x50x50 would represent a big computational challenge on 32 processors). All these difficulties will be addressed as developments continue.
Results of the SOPALE model calculations have been presented at several meetings,including:
and will be presented at:
The paper Pysklywec, R.N., C. Beaumont and P. Fullsack, Modelling the behavior of the continental mantle lithosphere during plate convergence is in press in Geology.
 Gilbert, J.R., G.L. Miller and S.H. Teng. Geometric mesh partitioning: Implementation and experiments, SIAM J. Sci. Comput., vol.19, no.6, 2091-2110, November 1998.
 Gupta, A., J. Mahesh and V. Kumar. WSSMP: Watson Symmpetric Sparse Matrix Package, The User Interface: Version 3.7.4, IBM Research Report RC 20923 (92669), July 17, 1997.
 Blankenbach, B., F. Busse, U. Christensen, L. Cserepes, U. Hansen, H. Harder, G. Jarvis, M. Koch, G. Marquardt, D. Moore, P.L. Olson, H. Schmeling and T. Schnaubelt. A benchmark comparison for mantle convection codes, Geophys. J., vol. 98, 23-28, 1989.
 Silvester, D. and A. Wathen. Fast iterative solution of stabilized Stokes systems, Part II: Using general block preconditioners, SIAM J. Numer. Anal., vol. 31, no.5, 1352-1367, October 1994.
 Vincent, C. The influence of the stabilization parameter on the convergence factor of iterative methods for the solution of the discretized Stokes problem, International Journal for Numerical Methods in Fluids, vol. 20, 1237-1252, 1995.
 Zienkiewicz, O.C. and R.L. Taylor. The Finite Element Method, 4th ed.,Vol 1: Basic Formulation and Linear Problems, McGraw-Hill, 1989.