Lithoprobe Supporting Geoscience Proposal, November 1998

• INTRODUCTION
• GENERAL PROGRESS REPORT ON LITHOPROBE SUPPORTED RESEARCH
• THE DALHOUSIE IBM SP-2 COMPUTER
• PROGRESS REPORT ON PROJECT A1 Figure 1
• PROGRESS REPORT ON PROJECT A2
  • A2.I NUMERICAL TOOLS Figure 2
  • A2.II A NEW APPLICATION PROGRAM : OPALE Figure 3, Figure 4
  • A2.III SIMPLE NUMERICAL EXPERIMENTS
• PROPOSED RESEARCH
  • DEVELOPMENT OF MODEL TOOLS
  • APPLICATIONS
• RELEVANCE

Back to Dalhousie Geodynamics home page

INTRODUCTION - During 1998 I was funded by Lithoprobe for research on two related PanLithoprobe proposals:

Most of the Earth’s continental crust consists of deformed metamorphic rocks that formed deep within mountain belts. Lithoprobe has produced seismic images of continental crust within mountain belts ranging in age from over 3 billion years old to the present day. Understanding how, where, and when deformation and metamorphism takes place contributes substantially to how these seismic images are interpreted, as well as to a general understanding of how continents and mountain belts have evolved throughout geologic time.

The Dalhousie Geodynamics Group, in collaboration with other members of the Canadian Institute for Advanced Research - Earth System Evolution Program, has developed numerical modelling techniques that provide computer simulations of the mechanical and thermal behaviour of the Earth’s crust during plate collisions. These ‘mantle-subduction’ models of orogenesis are based on the proposition that the mechanics of small orogens are controlled by sub-orogenic subduction and surface erosion (Willett et al. 1993, Geology 21, 371-74; Beaumont et al. 1994, Tectonophys. 232, 114-32). The model results demonstrate styles of crustal-scale deformation that have proven useful in interpreting data from a number of Lithoprobe transects and other orogens (eg., Beaumont & Quinlan 1994, Geophys. J. Int. 116, 754-83; Ellis et al., 1996, ECSOOT transect meeting report; Schmid et al. 1996, Tectonics 15, 1036-64). The model orogens are doubly vergent, with two crustal-scale thrust zones (termed the ‘pro-" and ‘retro-" shear zones) developed above the point (S) where sub-orogenic lithosphere is detached and subducted.

During 1996-97 Lithoprobe Supporting Geoscience 50 percent funded Dr. Susan Ellis as a postdoctoral fellow at Dalhousie. The main objectives of her research were to use the complete range of numerical model results from the mantle-subduction models to develop a series of templates for processes at convergent boundaries, and to apply these results to various contractional orogens studied by Lithoprobe.

Last year, we reported on: 1) Four Lithoprobe reports from this work; 2) Three Ellis-Beaumont papers (see form 100) and progress on; 3) The assembly of the model templates in an overall comparison with some of the smaller orogens studied by Lithoprobe, and; 4) A model study of mass partitioning of accreted material at convergent margins (some of these results were shown as Figures 2.4-3 and 2.4-4, in the Lithoprobe Phase V Proposal, pages 2-60 and 2-61).

I am now pleased to report that: 1) The overview paper, Ellis, S., Beaumont, C., Models of convergent boundary tectonics: implications for the interpretation of Lithoprobe data, has been submitted to Can. J. Earth Sciences; 2) The paper on mass partitioning at convergent margins, Beaumont, C., Ellis, S., and Pfiffner, A., Dynamics of subduction-accretion at convergent margins: short term modes, long-term deformation and tectonic implications, has been submitted to Jour. Geophysical Research; 3) A second paper concerning crustal scale episodic accretion in subduction zones, and the formation and stacking of fold-nappes has also been completed; Ellis, S., Beaumont, C., Pfiffner, A., Geodynamical models of crustal scale episodic tectonic accretion and underplating in subduction zones, submitted to Jour. Geophysical Research (now accepted subject to revision); 4) A third paper on the subduction-collision transition and prototype models of the evolution of the Alps, Pfiffner, A., Ellis, S. and Beaumont, C. is in preparation; 5) Additionally, the paper, Ellis, S., Beaumont, C., Jamieson, R., and Quinlan, G., Continental collision including a weak zone - the vise model and its application to the Newfoundland Appalachians, is in press in Can. J. Earth Sciences. Lithoprobe supported Susan Ellis in this research.

Lithoprobe Supporting Geoscience has also supported the development and application of our research on the coupling of mechanical and thermal processes during orogenesis. One aspect of this research addresses the source of heat for the ubiquitous Barrovian style of metamorphism in orogens. Barrovian metamorphism, with its characteristic progression from chlorite to sillimanite zones in pelites and greenschist to amphibolite in metabasites, is the most common style of regional metamorphism in time and space. Virtually every orogen covered by Lithoprobe transects includes at least some metamorphic rocks of this type (e.g., Southern Cordillera, parts of the eastern Grenville Province, parts of the Newfoundland Appalachians). Although it is widespread, this style of metamorphism is not adequately explained by existing thermal models of orogenic processes - in particular, most models do not account for observed temperatures of 600-700° C at depths of 20-30 km. The required crustal temperatures can be replicated by models that assume high crustal heat production (e.g., Chamberlain & Sonder 1990, Science 250, 763-9) or that consider tectonic accretion and redistribution of radioactive material (e.g. Huerta et al. 1996, Science 273, 637-39; Jamieson et al. 1996b, Huerta et al., J. Geophys.Res., 103, 15287-302, 1998). This approach is supported by data from some orogens indicating high heat production in metasedimentary sequences (e.g. Omenica Belt; Lewis et al. 1992, CJES 29, 1197-1214). In model orogens that incorporate efficient erosion, loss of radioactive upper crust partly offsets the effect of high heat production. However, radioactive material that is tectonically accreted beneath the crust, e.g. by partial subduction of crustal material ("tarm"), is protected from erosion and can act as an effective heat source for regional metamorphism in the lower to mid-crust. This effect is partly offset by continued subduction beneath the orogen as long as convergence continues at normal plate tectonic rates. Results from models involving a progression from subduction to accretion to collision, where some radioactive material is accreted to the base of the orogen, produce P-T-t paths that resemble those characteristic of Barrovian sequences, with early high P - low T metamorphism overprinted by higher T and lower P assemblages. As the lower crust heats up, the basal detachment propagates toward the foreland, and a tectonically stagnant zone develops immediately above the "tarm"; temperatures high enough to form granulite are eventually reached in this zone. Partial melting of "tarm";within 20 My of the end of convergence could yield the late- to post-orogenic granites observed in many Barrovian terranes.

The conceptual framework described above is based on quantitative experiments that used a coupled thermal-mechanical finite element model (developed by Philippe Fullsack). The results of the model experiments and the conceptual interpretation are described in, Jamieson, R.A., Beaumont, C., Fullsack, P., and Lee, B., 1998. Barrovian regional metamorphism: Where's the heat? (see form 100).

Our ability to develop and apply numerical techniques (eg. for Arbitrary Lagrangian-Eulerian (ALE) finite element calculations) requires a moderate level of computational capability. At a minimum, numerical experiments at medium levels of resolution require current state of the art single processor workstations (specfp 95>25, RAM > 256Mb), or preferably multiprocessor parallel computers (for example the IBM SP-2). The Dalhousie Geodynamics group was a partner with Dalhousie Unviersity (University Computing and Information Services (UCIS)) and other researchers in the acquisition of a small SP-2 computer (4 processors each with 66 MHz wide nodes, 128Mb RAM/processor). In the 1997/98 NSERC Major Equipment competition I was successful in obtaining partial funding (see form 100) for 4 more processors (166 mhz P2SC thin Nodes, 256 Mb RAM/processor). These new nodes (specfp95=25) are approximately 2.5 x faster than the old ones.

To maximize the advantage of the 4 new nodes required the complementary development of a parallel version (OPALE) of the computational fluid dynamics programming environment (MOZART) already completed by Philippe Fullsack. The development of OPALE has been Philippe's major focus during 1998. Results using 4 parallel P2SC nodes are discussed later in this proposal.

OPALE will be even more computationally efficient on larger parallel computers and we have joined a Canada Foundation for Innovation consortium proposal, with Jacob Slonim, Dal's dean of the new Faculty of Computer Science as P.I., to obtain a 16 node IBM SP-2 (or an equivalent or better computer from another vendor). This proposal has advanced to the second round of the competition.

In our 1997 proposal we described and illustrated results from ALE finite element calculations of crustal scale deformation. We proposed two subprojects which use and further develop these capabilities:

1) INVESTIGATION OF STRAIN SOFTENING/HARDENING IN MODEL OROGENS;

2) TO INVESTIGATE THE DEVELOPMENT OF LARGER-SCALE MODEL OROGENS AND COMPARE THE RESULTS WITH THE LARGER OROGENS STUDIED BY LITHOPROBE

In 1, we proposed that the next logical step is to investigate the consequences of deformation-dependent changes of material properties on the evolution of the models. Such changes include strain-dependent hardening and softening of both frictional (plastic) and ductile (power law viscous) rheologies in the models, melting or other mineralogical changes, and the feedback effect of shear heating on temperature-dependent rheological properties. The initial approach will be to use parametric models that relate rheology to strain and/or heating in a simple way. Although such models will not correctly represent rock behaviour, they will allow a systematic sensitivity analysis of model response to limited material property changes.

An important tectonic application is to determine the necessary changes in material properties that will lead to crustal-scale changes in tectonic regime. For example, it is generally agreed that ‘crustal weakening’ can lead to syn-convergent extension in orogens (orogenic collapse). However, the relationships between crustal strength measures and gravitational/buoyancy forces that lead to this behaviour are poorly understood. For example, is the late-tectonic normal faulting in the orogenic core of the Southern Cordillera caused by crustal weakening (?pegmatites/migmatites), or was extension caused by a change in boundary conditions? An equivalent problem exists for the Ontario part of the Grenvillian orogen which exhibits extension between contractional episodes.

In 2, we proposed that these crustal-scale models could be used to investigate more complex orogens with relatively large size (sizes comparable to the Grenvillian and Cordilleran orogens) involving several phases of subduction, accretion and collision, and large convergence. The models also have sufficient resolution to address the internal dynamics of the model crusts that are heterogeneous, either as a consequence of their initial properties or their evolution.

Specifically, we proposed a systematic investigation in which complexity is progressively introduced to the models. Complexities to be considered include: layered crustal properties (both Coulomb and thermally-activated ductility);random heterogeneity (introduced statistically); intrusives (eg. existing large-scale plutons and syn-tectonic intrusives and melting); surface denudation and deposition.

We illustrate progress with results from one particular line of research - an investigation of whether tectonic accretion of radioactive material (the "tarm" process described earlier) increases crustal temperatures sufficiently to cause melting and the possible consequences of this melting. These results were presented at the recent GSA Annual meeting in Toronto. The results (Figure 1) demonstrate the conclusion that lithospheric temperatures during orogenesis are controlled by the thermal Peclet number of the system (which measures the relative rates of heat transport by advection and diffusion, P_e=vL/k ) and, more importantly, the Damköhler III number (which measures the relative rates of heat production, in this case by tarm, and heat advection). When lithosphere subducts beneath an orogen, such that P_e>1, orogenic crust will normally be cooler than the equivalent crust that conducts the mantle heat flux to the surface. This is the well-known cooling effect of the over-riding lithosphere at subduction zones and leads to low temperature (blueschist) and high pressure (eclogite) metamorphism, not Barrovian or granulite metamorphism, or crustal melting. However, if the tarm process leads to sufficient increases in radioactive internal heating within the orogenic crust, and DIII>1, internal heating will more than compensate for the loss of conductive mantle heat flux owing to subduction beneath the orogen. A complete range of orogenic metamorphic and structural styles can occur depending on the values of P_e and DIII, and these styles will record the evolution of these dimensionless quantities.

The example shown here considers accretion of upper continental crust with A=2m W/m³ and demonstrates crustal temperatures >900° C (Figure 1). These conditions lead to melting and lower crustal convection in the model. Such conditions have been proposed for the crust of the Tibetan plateau and conditions modified by melting certainly occurred in the Grenville orogen and the Canadian Cordillera, although they were not as extreme as those shown in the model. The present model is overly simplistic. The melting model has to be made more physically realistic to consider melt production and extraction from the residue and the corresponding changes to the bulk properties of the system.

Melting and the ensuing changes to crustal rheology, heat transport processes, and density are generally considered to be intimately related to the onset of extension within orogens (often termed 'extensional' or 'gravitational collapse'). Extension, if characterized by normal faulting, can also be shown to occur during continuous plate convergence (eg. Himalayan-Tibetan orogen). Models like the one described above exhibit outward lower crustal flow leading to thinning of the crust beneath the model plateau (Olivier Vanderhaeghe, pers. comm.). Although this is a preliminary result, it may, if confirmed, provide insight into the ways gravitational potential energy is released/redistributed in large orogens.

In last year's proposal we showed the results of two different types of medium resolution ALE finite element calculations of deformation using the serial (one processor) MOZART programming environment. These model experiments were chosen to illustrate the computational capabilities of the techniques. The first type of experiment was the multilayer Coulomb model forced by basal subduction. The second experiment concerned a Lithospheric-Upper Mantle scale model, designed to provide a better physical model for subduction/collision processes in both oceanic and continental zones of convergence. One motivation for this type of model is to investigate the circumstances under which a subduction-type boundary condition (the type that we have used in our earlier modelling) actually occurs. The model also has a broader range of capabilities.

Project A2 focusses on the development of improved computational techniques designed to enable numerical experiments like the Coulomb and Lithosphere-Mantle models. We seek the maximum resolution, timestepping and large deformation that can be achieved on the Dalhousie SP-2.

Progress on this project includes 3 parts: I) The development of numerical tools for parallel computers (ie. 'parallel' MOZART capabilities) in addition to the ‘serial’ ones reported last year; II) Assessment of the performance/efficiency of the application program OPALE, and III) Illustrative simple model experiments that demonstrate some of these capabilities using OPALE.

A2.I   NUMERICAL TOOLS

The MOZART programming environment has been substantially enriched this year. MOZART is the software server we have developed to provide efficient numerical functions for application programs such as those used for the simulations of various thermomechanical models. It has evolved towards a platform interfacing high level libraries conceived and distributed by other research groups specialized in the field of scientific computing. In short, MOZART allows its users (program developers) to tune the algorithmic solution to the specific numerical problem at hand. The following diagram (Figure 2) shows the general organization of MOZART. Complex applications such as OPALE require the participation of a large fraction of the whole library totalling more than 100,000 lines of coding.

Data parallelism has been the major focus of MOZART'S development this year and we have reached the goal of writing a fully scalable finite element, arbitrary Lagrangian-Eulerian CFD code, OPALE, described in more detail below. To achieve scalability, parallel solutions had to be found for matrix assembly, linear system solutions and Lagrangian advection/regridding to Eulerian fields.

Linear System Solution for Incompressible Fluid Velocities - previous tests have convinced us that direct linear solvers are superior, in 2D, to iterative solvers and we have selected the PSPASES library (copyright IBM Corporation and University of Minnesota) as the best choice available today for the direct solution of large sparse symmetric positive definite linear systems. PSPASES was written by Mahesh Joshi, U of MN, Anshul Gupta, IBM Corp., and George Karypis, U of MN. All phases of the computation, namely the symbolic factorization, numerical factorization and triangular solvers are parallel and scalable. Excellent performances, up to 52 GFlops, have been recorded. This solver is portable on both distributed and shared memory computers. Reordering of the equations is performed in parallel by the PAROMETIS library using sophisticated graph tools and multilevel nested dissection in order to both reduce the amount of fill-in resulting from Cholesky factorization and produce a balanced elimination tree that is mapped onto the array of processing elements. While the minimum fill problem is NP-complete, nested dissection orderings have been proven to achieve fill within O (log(number of equations)) of minimum, and multilevel nested dissection offers a trade off between overall computation complexity/storage requirement and operation count balance across processors. The ordering phase therefore both reduces computational complexity and prepares data for parallelism. Matrix graph structures may show high compressibility. For example, the banded Linpack factorization of the velocity system discretized on a 600x200 grid (see the results below) would have required around 800 Mb of RAM, to say nothing about factorization time.

PSPASES uses column stripes partitioning and a multifrontal approach for the elimination process. This technique, intermediate between delayed-update (or fan-in) and immediate update (or fan-out), allows a better use of the memory hierarchy. Supernodal elimination then benefits from dense, vendor highly optimized, BLAS-Level 3 kernels and the algorithm becomes both parallel, partly vectorized, and tunable to the machine architecture.

To assess sparse direct solvers, such as PSPASES, we have written a MOZART interface which generates the matrix structure graph from any, structured or unstructured finite element grid. Serial functionality has been kept for BLKFCT, (named genspd in our previous reports) with the addition of an iterative method to allocate memory for unpredictable matrix compression rates.

A parallel matrix assembler has been developed, based on geometric domain splitting and either use of bissection or redundancy to relocate the matrix element within the matrix structure.

Several parallel advection/search/regridding algorithms have been written for both structured and unstructured grids.

A2.II    A NEW APPLICATION PROGRAM : OPALE

OPALE (Optimized Parallel Arbitrary Lagrangian Eulerian) is a fully scalable finite element program using the functions described above. It has both serial and parallel pre- and post-processors and is specialized to tensor product meshes of quadrilaterals. Various versions are in development and we show some model results in Figures 3 and 4. We have exploited the special character of the grid to write a quasi-1D version of the routines tracking the Lagrangian markers, resulting in a big decrease in computation time. We believe OPALE has a real potential for high resolution 2D thermomechanical models provided good access to appropriate computing facilities is available.

The following table gives an idea of some OPALE performances, on Dalhousie's IBM SP-2 (4 processors, each P2SC-166MHz-256Mb RAM). The global Lagrangian grid size is 1221x401 for all cases.

NEQ: Number of equations N1~26000
LS: Linear solve (=1 numerical factorization and solve) (s = seconds)
A: Assembly phase (pressure, viscosity and matrix computation included)

This table is likely to give an underestimate of OPALE's performance because the solver has not been finely tuned and assembly/search/regridding algorithms are still being improved. Furthermore, we expect PWSSMP (for which we have written an interface) to roughly double the linear solver performance on SP-2 machines. The 70% increase in solver efficiency is due to more efficient vectorization in the supernodal elimination.

We have successfully tested OPALE for 2 problems where analytical solutions are available: punching of a Von Mises material and initial Rayleigh-Taylor instability. Tests to come will include Boussinesq convection, Stokes flow and Coulomb critical wedges.

OPALE, being parallel and scalable, would be best used on larger machines with grid sizes further increased. The minimum resolvable length scale may be reduced or the size of the global geometry increased, or both. In the first case physical models relevant to previously unreachably small length scales (eg. rheological laws related to localization) can become part of the continuum model. The second alternative is illustrated below with a series of lithospheric scale subduction models.

A2.III     SIMPLE NUMERICAL EXPERIMENTS

We now discuss a series of numerical experiments (Figures 3 and 4) performed with OPALE which investigate how subduction may proceed, after some hypothesized initiation, in different modes, according to an assumed mechanical and density state of the lithosphere. The grossly simplified physics used here should not mask the main feature of this type of model: the approach is nonpertubative and can be used to solve problems ranging from purely kinematically driven (e.g. our previous models with prescribed basal velocities assuming a sharp ‘singular (s)’ transition) to purely thermally driven (eg. thermochemical convection), with a consistent treatment of the free surface.

The models presented here use a multiphase viscoplastic flow driven by some lateral entry flux (that can eg. be taken as a proxy for ridge push) and slab-pull induced by a negatively buoyant root. The pro- and retro-mantle are assumed to be separated by a transition zone of assignable strength. The view taken here is that motion is somehow being imported (outside of the domain under consideration) to the upper lithosphere which acts as the active part of the system while the upper mantle resists passively and viscously to slab penetration.

Viscoplastic bending, buckling and stable or unstable Rayleigh-Taylor effects may occur in various combinations depending on parameter choices. Introducing lateral slab heterogeneity (density and/or strength nomaly) may give rise to necking, or localized thickening, or slab breaking. [Some models not shown here also look at the mechanics of blocked subduction which is one proposed mechanism for obduction zones]. One can also look at the peeling off of positively buoyant subducted material as an alternative mechanism as proposed by Chemenda. High resolution models clearly show that crustal deformation accommodated by brittle faulting is resolvable, without the arbitrarily imposition of basal velocities, say at the Moho level. Trench migration with slab retreat or advance is an output of the model when active convection is not important.

Proposed research can be divided into two categories.

DEVELOPMENT OF MODEL TOOLS

APPLICATIONS

Here we briefly list some of the tests required and applications that are directly related to the interpretation of Lithoprobe data.

Although these are certainly general problems, they address many of the processes that are important to the overall interpretation of Lithoprobe results.

RELEVANCE - Geodynamical modelling is recognized as an important component of Lithoprobe (Phase IV/and Phase V proposals). As explained in the Phase V proposal, an important component of the PanLithoprobe synthesis will be to gain an improved understanding of the assembly of North America based on process models. This proposal is submitted in the belief that models derived from this research will be valuable to this process-based understanding.