Computational approaches to studying non-linear dynamics of the crust and mantle
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 to elements (approximately 10 million unknowns). It demonstrates good scaling behaviour, with scalability to 64 processes improving with problem size.
Problem sizes greater than 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)
- et al.
Dynamical Lagrangian remeshing (DLR): a new algorithm for solving large strain deformation problems and its application to fault-propagation folding
Earth Planet. Sci. Lett.
(1994) - et al.
Numerical modelling of strain localisation during extension of the continental lithosphere
Earth Planet. Sci. Lett.
(2001) - et al.
The accuracy of finite element solutions of stokes’ flow with strongly varying viscosity
Phys. Earth Planet. Inter.
(1996) - et al.
A Lagrangian integration point finite element method for large deformation modeling of viscoelastic geomaterials
J. Comput. Phys.
(2003) - et al.
A numerical method for suspension flow
J. Comput. Phys.
(1991) - et al.
Antisymmetric form of the material point method with applications to upsetting and Taylor impact problems
Comput. Meth. Appl. Mech. Eng.
(1996) - et al.
A study of some numerical viscoelastic schemes
J. Non-Newton. Fluid Mech.
(1991) - et al.
Non-linear effects from variable thermal conductivity and mantle internal heating: implications for massive melting and secular cooling of the mantle
Phys. Earth Planet. Inter.
(2002) - Appelbe, W.F., Moresi, L., Quenette, S., Sunter, P.D., 2006. Scientific software frameworks and grid computing:...
- et al.
Efficient management of parallelism in object oriented numerical software libraries