Computational approaches to studying non-linear dynamics of the crust and mantle

https://doi.org/10.1016/j.pepi.2007.06.009Get rights and content

Abstract

We outline a mathematical formulation for mantle convection which can deal with the viscoelastic–plastic rheology of the cool parts of the lithosphere. This formulation is then analyzed to expose the numerical challenges inherent in the equations and a suitable solution strategy is outlined. With this strategy in place, we discuss a parallel implementation, explaining how we maintain computational efficiency, in a tiered and modular software environment. We show two examples from geodynamics which demonstrate the capability of the numerical method and its implementation.

Introduction

The dynamic evolution of the Earth at the large scale – mantle convection, plate tectonics and the deformation of continental crust – is a problem which requires both massively parallel computing and advanced numerical techniques. Computational challenges associated with modelling the Earth’s large-scale evolution are dominated by the need to model the non-linearities in the system, the emergent fine-structure which they generate, and to resolve the important range of length and timescales in any model. Specialized numerical approaches can deal with each of these issues very effectively, but efficient implementations of these techniques for parallel architectures are difficult to implement.

The least well-constrained information about the global system of plates is the mechanical behaviour of the lithosphere, particularly the deforming boundary zones. The reason for this is three-fold: (1) the constitutive laws for minerals under lithosphere and mantle conditions (temperatures, pressures and strain rates) are difficult to measure in the laboratory; (2) the precise composition of most of the lithosphere is unknown and heterogeneous; (3) plate boundaries are, in fact, systems with internal structure and evolution in addition to those associated with rock deformation at the cm scale—the processes which are active at the crystal scale interact with those which are active at the lithosphere scale (e.g. the focussing of ground water into shear zones in the crust). In continental regions, deformation is generally more diffuse, but the same difficulties apply when we try to model the details of the structures which accommodate deformation.

Constitutive laws for the lithosphere and mantle are empirical, based on understanding the small-scale physics, and constrained by observations; they depend upon the scale of the model. England and Molnar (1997) showed that a non-linear viscous fluid model gives a good description of the orogen-scale deformation of the Himalayas but not the formation of the individual fault systems. To capture the tendency for the lithosphere to localize deformation, Weinstein and Olson (1992) suggested using very non-linear viscosity models which limit the maximum stress in the lithosphere; and Bercovici, 1993, Bercovici, 1995 advocated a rheology which weakens at large strain rates based on a physical model in which the interaction of voids with interstitial fluids creates a self-generating weakness. Strain softening, in which strength is lost due to accumulated damage during deformation, introduces persistence of localization across discrete phases of deformation. Such models are common in material science, but relatively recent in lithosphere-scale geological modelling Lyakhovsky et al., 1997, Lavier et al., 1999, Lenardic et al., 2000, Huismans and Beaumont, 2003. Regenauer-Lieb et al., 2004, Regenauer-Lieb et al., 2006 and Kaus and Podladchikov (2006) demonstrated how thermal–mechanical feedbacks create short-timescale responses in localized deformation.

When small-scale processes below model resolution are parameterized in a larger-scale constitutive law, the resulting model is usually non-linear, anisotropic, and strongly history dependent. Many techniques have been developed by the computational mechanics community to address similar problems to those faced in the geosciences—fluid/solid systems with viscoelasticity, frictional failure surfaces, emergent anisotropy, strain or damage dependent material properties, and long-lived evolving interfaces. These have been taken up, and extended by computational geodynamics groups to cater for the specific modelling needs of the discipline (e.g. Melosh, 1978, Poliakov et al., 1993, Braun and Sambridge, 1994, Batt and Braun, 1997, Toth and Gurnis, 1998, Moresi et al., 2001, Moresi et al., 2003, Schmalholz et al., 2001, Vasilyev et al., 2001, Frederiksen and Braun, 2001, Muhlhaus and Regenauer-Lieb, 2005). The Lagrangian Integration Point FEM Moresi et al., 2001, Moresi et al., 2003 is one of the conceptually simpler of the specialized methods for geodynamics—it tracks material deformation on moving material points which are then used to create the integrals in the finite element weak form.

Hence a suite of requirements are presented for fabrication into a computational tool that span geophysical workflow, constitutive and numerical concerns. With clear indications that present community codes cannot achieve these requirements Harry et al., 2005, Zhong et al., 2005, Spiegelman and Montesi, 2006 without substantial effort, and realizing there is a need to further evolve the implemented constitutive laws and numerical methods to ultimately achieve these goals, a new approach to computational software development is taken. This problem of obsolete or hard to maintain codes is not unique to computational geophysics, as equivalent problems are evident in programs within the NSF, NASA, ASCI and DARPA (Post and Votta, 2005). Our approach is to enable a spectrum of computational scientists to collaboratively develop the prescribed capabilities. In particular we identify how our over-arching architecture helps nurture the development of advanced numerical schemes, readily integrated into these scale crossing geophysical models.

Section snippets

Problem description

Here we outline a mathematical and numerical formulation for a simplified picture of the coupled thermal/mechanical evolution of the Earth’s crust and mantle. We include the rheological complexity associated with elasticity, non-linear temperature and pressure dependent viscosity and plasticity, and include a description of history dependent, anisotropic plastic damage. We neglect many other important effects which broaden the parameter space considerably but which do not increase what we might

Computational approach

We have expressions for the stress balance, mass conservation, temperature, and strain history. We first turn our attention to the challenges to accurate and efficient numerical solution inherent in this system of equations, and then address our particular numerical approach to these problems.

Software approach

The mathematical and numerical approaches outlined above deliberately make very little concession to possible software implementation issues. This helps to ensure that the modelling is dictated by the applicable mathematical representation of the physics and/or chemistry and not by the availability of an efficient software tool needing an application. In this section, we discuss how such general purpose software can be constructed which is still efficient on a range of parallel architectures.

Performance

Fig. 2(a) is a parallel speedup graph of the titled plume conduit example in Section 6, performed on the SGI Altix 3700 Bx2 cluster the Australian Partnership for Advanced Computing (APAC). These are 3D problems of mesh sizes 32×32×32 to 128×128×128 elements (approximately 10 million unknowns). It demonstrates good scaling behaviour, with scalability to 64 processes improving with problem size.

Problem sizes greater than 128×128×128 elements are too large to run serially on the controlled

Example applications

In this paper we describe two sample applications of our software framework which are each computationally demanding for very different reasons. The first example is the rise of a thermally or chemically buoyant near-horizontal cylinder of viscous fluid in a viscous background representing a tilted mantle plume conduit. The second example is the stretching of a layer of viscous fluid underlying a second layer of viscoelastic–plastic material which represents the extensional deformation of the

Conclusion

The problem-size requirements on our computational tools span at least four orders of magnitude, as shown in the example of the basin extension models, ranging from 9000 elements, to a desired 90 million for high fidelity geological relevance. Simultaneously, the complexity and computational demands of the constitutive models have evolved from simple geometrical descriptions of rheological domains to fully coupled models which result in emergent structure on a range of scales and produce steep

Acknowledgements

Funding for StGermain–Underworld and gLucifer has been predominantly provided by VPAC, Monash University, the ACceSS MNRF Snark and data assimilation projects, and the APAC CTT program. We would also like to thank the contributions from the Computational Infrastructure in Geodynamics, Caltech’s GeoFramework NSF ITR, and Deakin University’s Xanthus project. The VPAC development team also includes Patrick Sunter, Alan Lo, Luke Hodkinson, Raquibul Hassan, Kathleen Humble and John Spencer. The

References (64)

  • G.E. Batt et al.

    On the thermomechanical evolution of compressional orogens

    Geophys. J. Int.

    (1997)
  • D. Bercovici

    A simple-model of plate generation from mantle flow

    Geophys. J. Int.

    (1993)
  • D. Bercovici

    A source-sink model of the generation of plate-tectonics from non-Newtonian mantle flow

    J. Geophys. Res.: Solid Earth

    (1995)
  • B. Boroomand et al.

    Recovery by equilibrium in patches (REP)

    Int. J. Numer. Meth. Eng.

    (1997)
  • A. Brandt

    Multigrid techniques: 1984 guide, with applications to fluid dynamics

    GMD-Studien, vol. 85

    (1984)
  • A. Brandt

    Multi-level adaptive solutions to boundary value problems

    Math. Comput.

    (1977)
  • S.J.H. Buiter et al.

    The numerical sandbox: comparison of model results for a shortening and an extension experiment

  • H.P. Bunge et al.

    A sensitivity study of three-dimensional spherical mantle convection at 10(8) Rayleigh number: effects of depth-dependent viscosity, heating mode, and an endothermic phase change

    J. Geophys. Res.: Solid Earth

    (1997)
  • U.R. Christensen et al.

    Layered convection induced by phase transitions

    J. Geophys. Res.

    (1985)
  • P. England et al.

    Active deformation of Asia from kinematics to dynamics

    Science

    (1997)
  • Harlow, F.H. A machine calculation method for hydrodynamic problems. Technical report LAMS-1956. Los Alamos Scientific...
  • Harry, D., Lavier, L., Willet, S. Report to CIG from the NSF workshop on tectonic modelling. Technical report....
  • R. Huismans et al.

    Symmetric and asymmetric lithospheric extension: relative effects of frictional-plastic and viscous strain softening

    J. Geophys. Res.

    (2003)
  • R.S. Huismans et al.

    The effect of plastic-viscous layering and strain-softening on mode selection during lithospheric extension

    J. Geophys. Res.

    (2005)
  • S. Karato et al.

    Rheology of the upper mantle—a synthesis

    Science

    (1993)
  • B.J.P. Kaus et al.

    Initiation of localized shear zones in viscoelastoplastic rocks

    J. Geophys. Res.: Solid Earth

    (2006)
  • R.C. Kerr et al.

    Structure and dynamics of sheared mantle plumes

    Geochem. Geophys. Geosyst.

    (2004)
  • L.L. Lavier et al.

    Self-consistent rolling-hinge model for the evolution of large-offset low-angle normal faults

    Geology

    (1999)
  • A. Lenardic et al.

    The role of mobile belts for the longevity of deep cratonic lithosphere: the crumple zone model

    Geophys. Res. Lett.

    (2000)
  • J.R. Lister et al.

    The effect of geometry on the gravitational instability of a buoyant region of viscous fluid

    J. Fluid Mech.

    (1989)
  • V. Lyakhovsky et al.

    Distributed damage, faulting, and friction

    J. Geophys. Res.: Solid Earth

    (1997)
  • H.J. Melosh

    Dynamic support of the outer rise

    Geophys. Res. Lett.

    (1978)
  • Cited by (0)

    View full text