SOPALE : a finite-element computer code for the computation of visco-plastic creeping flows with applications to tectonics (short summary) ©P.Fullsack

      Method, limitations, history

1.      What is SOPALE?

SOPALE is a computer code written in Fortran 77 which computes the flow of visco-plastic materials in the ‘creeping flow’ limit, that is, in the limit where viscosities are so high that inertial accelerations may be neglected. It is specialized to the field of geodynamics which studies the thermo-mechanical behavior of the crust, the lithosphere and upper mantle on geological time scales, for which it provides a forward evolution simulation tool.

2.      How does SOPALE work?

The material flow is represented as a sequence of quasi-static equilibria between body forces (due e.g. due the presence of an external uniform gravity field) and the viscous resistance opposed by internal stresses in the deforming material. The material deformation is caused by the combination of applied kinematical boundary conditions and body forces, and limited by the material viscosity. The equations of continuum mechanics are used to cast the equilibrium problem as a Stokes flow problem, which is a saddle-point system of partial differential equations to be solved for the velocity and pressure of the flow. The computed velocity field is then integrated in time to follow the trajectories of the material.

The code embeds a choice among a few parametric constitutive laws relating stresses and deformation (named rheological laws). These laws attempt to mimic experimental observations and depend on the category of material studied.

In SOPALE, they have been somewhat specialized, although in very generic form, to

give a coarse description of the mechanical properties of rocks (see section 5.). It should be emphasized that these rheological laws treat the geological material (rocks) as fluid and not solid. SOPALE and MICROFEM (see section 9) were probably among the first codes written in this perspective (for tectonic applications), which was chosen to allow very large deformation, while the traditional elastic or elasto-viscoplastic models were typically much more limited in the total deformation that could be computed.

In the context of geophysics, temperature is a primary factor affecting the value of viscosities, and the material flows in a non-uniform temperature field (named geotherm, for the earth). The Stokes model must be supplemented with a thermal model describing the conservation of heat. SOPALE uses the mechanical (Stokes) velocities as the transporting speed (one assumption of the physical model used)    and combines transport with heat diffusion, limited by the material thermal conductivity, and heat sources, representing for example the effect of radioactivity, or shear heating. Heat is therefore produced, diffused and transported, and the ratio of heat per thermal mass, or temperature field, by a feed-back effect on the rheologies, contributes to control the mechanical behavior.

This gives rise to a coupled ‘thermo-mechanical’ system.

Constitutive laws coded in SOPALE  for the mechanics and the thermics reflect the primary type of applications that SOPALE has been used for, that is problems in geophysics involving crustal and upper mantle dynamics over time scales typical of orogenesis (see section 3) . Similarly, all external processes (see section 6) applied to the evolving flow deal with specific aspects of tectonics or geomorphology.

The initial geometry of the domain to be deformed is typically a rectangular window, or a perturbation of such, with an aspect ratio selected as a function of the processes and scales considered.

The equations are transformed to a finite matrix form (in the sense of linear algebra), through finite-element discretization (FE for short).

The FE method approximates the material continuum by a grid (nodes and mesh) and uses a variational formulation of the equations of continuum mechanics. The FE method is applied to both the mechanical and thermal equilibrium models (see section 4)

SOPALE performs a finite sequence of incremental motion computations (or ‘time steps’) to move and deform the material. The velocities and temperature are computed with the FE method on a mildly deformed (eulerian) grid, while another (lagrangian) grid (or cloud of points) follows the material motion. The method is at least as old as Harlow and Welsh first particle-in-cell computations in Los Alamos (1965), and is widely used in computational fluid dynamics. One refers to such methods as Arbitrary-Lagrangian-Eulerian or ALE methods

-         an eulerian grid serves as the FE discretization grid, and carries the velocity and pressure computed from the Stokes equations. The same grid is used to solve the thermal equations. This grid conforms to the total domain boundary by adjusting its thickness (to match the current positions of the top and bottom boundaries)

-         A lagrangian grid or set of points (lagrangian  markers or particles) is moved along material trajectories. This grid records the physical history of material particles, in particular the pressure, temperature and accumulated deformation.

This grid is highly deformed and not suitable for a FE computation, which is why a separate eulerian grid is used.

To summarize, the computation initializes a domain geometry and spatial properties (fields), then computes a sequence of steps, each composed of the same subsequence: solution of a mechanical equilibrium (FE-Stokes), a solution of thermal ‘equilibrium’ (FE-Temperatures), and an update integrating in time the lagrangian or material fields, and computing a matching image of these fields on an update eulerian grid (ALE)

3.      What are the applications of SOPALE ?

The technology developed in SOPALE has been designed to solve problems for specific applications, and are therefore narrowly bound to them.  In this section we want to outline a few motivations which were central to the creation of the code and its development over the past 10 years.

SOPALE has been applied to the flow of geological materials (rocks, sediments) at crustal scale (around 30 kilometers) and upper mantle scale (depth of of order a few hundred kilometers). It has been used to study deformation subjected to compressive (e.g. fold and thrust belts) or extensive (rifting, lithospheric extension, necking) boundary conditions.

To briefly outline a few examples, SOPALE has been applied to

-understand the coupling between surface processes (e.g. erosion and deposition, in subaerial or submarine environments) and tectonics. This coupling is a direct result of the pressure sensitive nature of rock rheologies and leads to rock exhumation and metamorphic transformations which highly depend on the narrow interaction between the mechanical, thermal, and geomorphological components of the system (see section 9). The ability to represent such interactions in SOPALE is one feature which distinguishes it from many other public domain or commercial codes.

-understand halokinesis , in studies looking at the response of the gravitationally unstable (heavy plastic layer over a lighter, weak salt layer) to a given parametric progradation sequence.

-understand the mechanics of subduction and roll-back, at upper-mantle scale, for example the balance of  forces such as trench pull, slab push,  the influence of small scale convection, the role of crustal and mantle rheologies on slab break-off.

Applications have been developed primarily, but not exclusively (see section 9), at Dalhousie University in the Geodynamics Group (hereafter mentioned as DGG) led by Professor Chris Beaumont.

They have been published in a series of publicly available research articles. More details on these publications and the authors of various applications may be found on the Dalhousie Geodynamics Group website  http://geodynamics.oceanography.dal.ca/index.html)

Each application has typically enriched the SOPALE code with a module typical of the study at hand. Bonny LEE has implemented many of these features (see section 6 and 9), with the DGG geophysicists who would specify the type of law to be implemented.

4.      Details on the FE computations

·        Mechanical FE solution (Stokes flow)

The resolution of the Stokes equations uses the bilinear velocity/discontinuous pressure or Q1-P0 finite-element. This interpolation supplements the isoparametric 4-noded quadrilateral with a pressure degree-of-freedom representing the average pressure in the element. This ‘internal or bubble’ pressure unknown is eliminated by an approximation of the incompressibility condition named penalization. This approximation replaces the exact requirement of null velocity divergence by a linear relation between pressure and divergence, using a coefficient (the bulk viscosity) much larger than the shear viscosity.

The elimination of the pressure unknown replaces the original, saddle-point, Stokes problem, by an elliptical problem, which leads to a symmetric positive definite ‘stiffness’ or FE matrix. Such systems always admit a square root triangular matrix (the Cholesky factor, L) and can be solved, once L has been computed (Cholesky factorization), by 2 triangular matrix solves (upper and lower back-substitutions). This is the major benefit of the penalized approach, at the expense of approximate incompressibility, compared to the full indefinite solver in the velocity-pressure variables.

As viscosities are functions of the velocity gradient field (strain rate), through the rheological laws, the velocity is found by re-solving several times (non linear iteration loop) the equilibrium between external forces and internal (visco-plastic)

forces, with an update of viscosities from the last predicted strain rate to account for yield stresses and power-law rheologies (see section 5)

In each non-linear iteration:

       -the body forces are computed. They combine the effect of temperature on rock densities (coefficient of thermal expansion at constant pressure, reference density,reference temperature), and ,for some applications, the weight of a water column loading the top ‘mechanical’ surface.

         -parametrized  ‘phase change’ may occur during iterations, depending on local thermodynamic (pressure/temperature) conditions.

         -discontinuous pressures and nodal pressures are computed from the last predicted velocities. An input flag allows the filtering of a spurious 0-energy mode asscciated with the Q1-P0 element (a mode named checkerboard mode or ‘element-by-element oscillating mode’). The mode is removed by L2 projection which amounts to computing the bilinear field closest in square energy norm to the discontinuous ‘elemental’ pressure field. Note that other quasi-singular modes , with a energy dissipation rate tending to zero in the limit of 0 grid size, exist for the Q1-P0 element, but are not filtered or stabilized in this code.

          -viscosities are updated, using discontinuous pressures and other element-averaged quantities (see section 5) controlling the rheologies

          -the FE stiffness matrix of each element is computed and assembled in the form required by the sparse linear solver BLKFCT (see section 10) (the right-hand of the linear system, a column vector, is also assembled). The code uses selective integration, and integrates the spherical part of the matrix with the one-Gauss-point integration rule, and the deviatoric part of the matrix with a four-Gauss-point rule.

Boundary conditions are enforced by penalization: at prescribed degrees of freedom, the diagonal of the matrix is augmented by a term much larger than the diagonal maximum, and a reaction using the same lagrangian multiplier and the prescribed value is added to the right-hand side vector. Boundary conditions are not discussed in detail here.

·        Thermal FE solution

The thermal FE calculation uses the standard 4 noded quadrilateral element with an streamline-upwind Galerkin (SUPG) stabilization. In the limit where heat transport is faster than heat diffusion, that is in the high Peclet number range, the Galerkin variational method becomes inadequate for the solution of the diffusion-advection equation. The presence of a dominant first-order advecting term destabilizes the Galerkin method. We use the SUPG method proposed by Hughes (and now widespread in the FE community) and compute upwind positions of the 4 standard Gauss integration points. The points are moved by a function of the local Peclet number (1/tanh(P)-1/P) which solves exactly the 1d advection-diffusion problem with uniform advecting velocity. This is equivalent to a Petrov-Galerkin method and stabilizes computation at high Peclet number.

A Cranck-Nicholson scheme is used for time discretaization and the advection matrix term responsible for the lack of symmetry of the FE matrix is reported in the right-hand side. The system matrix is therefore a combination of mass and diffusion, is symmetric positive definite, and uses the same Cholesky solver as the mechanical system. The right-hand side combines the complementary Cranck-Nicholson term (mass and diffusion), the updated effect of advection and the heat sources (e.g. radio-active production or frictional heating). Even for linear thermal materials (constant conductivity, thermal inertia) iterations are required to solve the (linear) non-symmetric problem.

·        Linear solver : BLKFCT (see section 10)

5.      Details on constitutive laws

Code users may prescribe several material zones, which will be advected (i.e moved with the material motion). New material may also we added as a mass flux, e.g. by filing depressions left by the tectonic deformation in the topography (top surface shape).

Each material zone is assigned a ‘physical material’, which is characterized by a particular thermal or mechanical constitutive law.

Mechanical constitutive laws are all visco-plastic. They are characterized by a bulk viscosity (a numerical parameters enforcing the penalization of incompressibility)  and an effective viscosity which depends on various controls (e.g. temperature, effective pressure. total accumulated strain). The deviatoric stress tensor and deviatoric strain rate tensor are proportional, have the same principal axes, and the effective viscosity measures their ratio.  The visco-plastic material are assumed to behave viscously until the norm of the deviatoric stress reaches a yield envelope, at which point deformation becomes visco-plastic and adjusts to the stress yield level.

Plane strain assumptions are built-in the model (the model plane , typically vertical, is assumed to contain the minimum and maximum principal axis, no deformation occurs in the direction  perpendicular to the plane).

The viscous behavior sums (as in parallel rheological circuits) the forces from 3 creeping mechanisms :

non-linear or power-law creep, grain-size sensitive creep, and Peierls creep. The forces are driven by the ambient strain rate and the various controlling fields (see below). The plastic yield stress obeys a Drucker-Prager law.

These laws give crude approximations to the complex rheological behavior of rocks

Controls on the local rheology include:

-         various parametric (non dynamic) forms of the fluid pressure. Fluid pressure may be prescribed e.g. by a user-defined function of strain (interpreted as having an effect on porosity) multiplied by the fluid hydrostatic pressure.

Terzaghi’s effective stress principle is implemented in various forms (e.g. the Drucker-Prager yield stress is an affine function of effective pressure)

-         parametric user-defined dependence of cohesion and angle of friction on total strain (implying strain hardening or strain softening

-         temperature (through pre-exponential activation energy factors e.g. in thermally activated dislocation creep) or pressure (activation volume)

-         the total accumulated strain (which integrates in time the lagrangian strain rate invariant, and should be considered as a measure of dissipated work, or of accumulated damage, rather than a geometrical deformation measure.

-         grain size

6.      Details on external controls

Users of SOPALE can provide information relative to the non-modeled part of the geological/geophysical system, by prescribing various forces, mass, or geometric loads in the form of imposed functions of time and space. The code response to these inputs needs to be checked and interpreted, as no consistency is required for these parameters.

Users may ,for example :

·        Subject the model to surface processes, which will either redistribute mass in a diffusive (curvature-dependent), or slope-dependent manner, or add mass to the system by partial or total filling of accommodation space below a reference height. (e.g. to mimic basin fill). The added material (sediment) may then be reworked by tectonics and lead to structures similar to e.g. inverted basins.

Sediment addition may take the form of uniform agradtion or a gaussian-shape progradation (progressing seaward)

·        Subject the model to a waterload, e.g. at an oceanic margin. This will later the surface processes (subaerial versus submarine), the fluid hydrostatic pressure and the total weight applied to the model supporting base

·        Subject the model to various boundary kinematical boundary conditions (velocities, with of without mass influx) and thermal boundary conditions (given temperature or heat flux)

·        Provide an external support to the model in the form of a broken beam, made of 2 elastic members meeting at a common point. This allows users to rest the model on an underlying elastic support (e.g. a piece of mantle) floating in an inviscid fluid. Models combining Airy and flexural isostasy are therefore available, and various rotation/deflection/force/torque boundary conditions may be applied to the beam ends. No opening is allowed however between the beams, which are assumed to join continuously. Flexural rigidities may be assigned per column.

7.      Limitations of SOPALE

What SOPALE can do is best documented in the published research papers using it (see section 3). We want here to explain what the code cannot do or is not guaranteed to do, and why. This will help preventing a misuse of the code, warning code users, and circumscribing the limit of the technology, which might be useful in issues involving distinctions between code improvements and new technologies.

We see 3 major sources of limitations to the present technology: mathematical model limitations, physical model limitations and algorithmic model limitations:

·        Mathematical model limitations

The rheologies used in SOPALE belong to the family of non-associated plastic models, for which the plastic potential and plastic yield functions differ. The standard tools of mathematical analysis (variational inequations and convex analysis) fail to provide a method for proving the existence and unicity of solutions in this case. The numerical model used in SOPALE can only be

seen as a pragmatic extension of a technique fully proven for other (e.g. linear, or associative plastic) rheologies, without much (if any) theoretical support/justification for the computed solution for the general non-linear case.

 The mathematical problem at hand is very challenging for it combines incompressibility, wide range of viscosities, and a yield stress which becomes null (singular) at the surface (for Coulomb plastic materials which have no cohesion). Most (all?) existing codes adopting the same method, face the same problem. This problem is aggravated by a lack of coherence in the formulation itself, i.e. the continuum mechanics description is unable to provide a consistent definition for faults (stress and velocity characteristics differ).

SOPALE can be used to approximately model  rigid-plastic systems. When the system is built with either one physical zone operating in the viscous range, or with a mix of viscous and plastic layers, the viscous part of the system ‘gives a time scale ‘ for evolution and  constrains the system to have a well defined solution. Hence we discuss here the more difficult question of rigid-plastic systems.

The traditional approach of solution uses a system of hyperbolic equations (Hilda Geiringer’s famous equations) which is not a boundary value problem, but a propagation problem. The theoretical question is: what validity does the plasticity penalization, boundary value method, used in SOPALE (at the continuum level) have for such problems? We note that rigid-plasticity  may be seen as the limit of a Bingham fluid rheology when the fluid viscosity is null. Bingham fluids use the Von Mises potential as a minimum stress below which no flow occurs. Duvaut and Lions have proven the existence and unicity of the interior Bingham fluid evolution problem, using the method of variational inequalities. However in SOPALE rigidity is not modeled precisely : it is approximated (regularized) by a user-chosen large viscosity below yield.

Known solutions for nonlinear plastic problems are rare and hard to find.  This prevents the calibration of numerical errors for problems involving these rheologies, that is for quasi-all problems in the range of applications of SOPALE.

It is therefore impossible to provide an estimation of numerical errors, and the code provides no indication of accuracy (even if it has of course been tested in a series of simple benchmark problems: Rayleigh-Taylor, thermal convection, critical wedge average slope , punching problem, advection-diffusion, simple corner flows…)

·        Physical model limitations

We have selected in this section a few examples showing some limitations of the physical model used in the code. These examples are far from exhaustive, and given only for illustration.

The code uses a plane strain assumption and models a two-dimensional vertical plane, without deformation in the direction orthogonal to the plane. It ignores the earth curvature even on the scale of 500-1000 kilometers. It isolates a ‘slice’ of material which is part of a larger system, therefore requires external assumptions

regarding the fluxes/forces applied to the modeled ‘physical window’.

Consequently, many external parameters are introduced (see section 6), without any restriction imposed by a system of equations (these ad-hoc parameters will give opportunities to users to bias/modify the results in various directions).

This implies that good users should be ‘field-specialists’

They should be aware of the meaningful ‘forcing’ time scales and spatial scales (e.g. in surface processes), and capable of interpreting the pertinence of numerical results.

Rheologies involve the use of the effective pressure supported by the rock matrix However, fluid pressures are parametrized without an explicit solution of the equations of Darcy, (involving explicit porosity, permeability model etc.)

Rheologies are fluid and do not capture the elastic behavior of rocks.

Material phase changes are allowed but there is no associated production or absorption of heat..

Grain size sensitive creep rheologies are available but no model of grain size evolution is available.

All models are isotropic, and remain so even for material which have been highly sheared along a constant plane direction (it is questionable that this material would not develop an orthotropic symmetry, even at large scale).

The role of hydration / dehydration is not taken into account in any rheology.

However this role is important:

-in the dynamics of subduction

 The water cycle is perceived to play a key role in the lubrification of plate tectonics and dehydration of the water-loaded subducting plate, or dehydration related to partial melting is believed to have a major role on viscosities and on the time of transport of diapiric intrusions

-in halokinesis: the rheology of evaporates and salt rocks is sensitive to water and the cycle of water may play a role in the timing of salt diapirism. We also note that no sediment lithification modeled is included in the code. Other research work (e.g. Tuncay and Ortoleva) have attempted to use ‘finer scale’ rheological models which couple the physics and chemistry of sediment diagenesis with the evolving (dynamic) rheology of the sediments. Quantitative timing of halokinesis may be a very tricky business, because the Rayleigh-Taylor instability will amplify any model deficiency!

-erosion, deposition, progradation are all geometric postprocessors of the surface shape. They typically do not conserve mass and there is no attempt to understand where material is eroded from or where it is transported. Chronostratigraphy is imposed as a user-defined sequence of sedimentation rates and material labelling corresponding to a predetermined cycle, irrespective of the evolution of the model. These postprocessors are therefore different from numerical models (such as Cascade) used in the Geophysics community to compute the evolution of the surface geomorphology.

Hence, again, computations are rather generic in the code, and should only be seen as qualitative.

·        Algorithmic model limitations

-limitations of the discretization used :

The code uses a very primitive mesher, which only allows for a tensor product (given number of rows and columns) rectangular mesh. The top and bottom mesh boundary may be initially perturbed. They follow at each time the domain boundary, as prescribed by their net motion (e.g. tectonic minus erosional). This is done by ‘thickening’ or thinning elements. This may lead to a loss of accuracy (element distorsion).

The eulerian grid only moves its nodes vertically (by thickening) to follow the outer shape.

There is no addition of nodes, no adaptation of the grid where gradients would be higher, or near physical material boundaries.

Thickening involves a loss of resolution due to discretisation of a greater height column by a constant number of elements. Fortunately, in crustal-scale tectonic applications, angles are quiet small the thickening bounded by a factor around 2. However in other models the loss of accuracy does become very important.

Another deficient aspect of the discretization is the lack of  any form of discretization for material zones boundaries, i.e. the limit between a physical material and another physical material (e.g. a layer of salt and a layer of sediments) is simply not represented. Each material zone is represented by an array of particles and  patching is used to avoid cases where an eulerian grid would not contain any lagrangian particle. Much better algorithms (e.g. using a lagrangian Delaunay grid) clearly exist and have been tried, but not in SOPALE.

Some eulerian finite-elements may contain particles belonging to different material zones. In this case a crude majority rule is used. This rule does not take into account any of the actual flow stress that the participant materials would actually compute !

-treatment of the coupling between temperature and rheologies

The code never solves a coupled equilibrium involving velocities, pressure and temperature. It simply alternates a mechanical and a thermal solution. A fully coupled treatment, or an iterative solution of this coupling has been implemented by other experts in the field (e.g. Peter Van Keken) and would compute a more accurate equilibrium

-limitation of the linear solver used

Better, faster sparse solver exist. If a parallel hardware platform can be used, the model spatial resolution may be increased (FE grids of 1000x1000, or equivalent) should be feasible on  32 processors machines (e.g. IBM  p690 or Linux clusters)

            This increased resolution may be used e.g. to include other physical scales in the model and/or to increase or test the model accuracy

One of the ‘co-inventors’ (PF) thinks that new technologies and continued research work will be needed to address most of the limitations mentioned in this section. 

8.      Where is the code documentation?

In addition to the present document, users of SOPALE may read the following manuals :

1. The sofware manual explains the code flowchart, with an emphasis on the concepts used in SOPALE.
2. The source manual explains how the source code is organized and discusses the static memory allocation of variables
3. The theory manual explains the numerical techniques used.
4. The user manual explains how to use the software (inputs,run,postprocessing)