Time Dependent Behaviour of Congested Subduction

Louis Moresi, School of Earth Sciences, University of Melbourne, Victoria, Australia
David Willis, School of Earth, Atmosphere and Environmental Sciences, Monash University, Victoria, Australia

Introduction

The movies in this supplement are based on work described in [1] and [2]. We model a wide, three dimensional subduction zone which we collide with a zone of thickened crust representing, for example, an oceanic plateau.

These movies highlight three aspects of the time-dependent dynamics of subducting slabs. First are the fluctuations in the rate of plate motion caused by geometrical instabilities in the slab as it hits the transition zone and bends or buckles or by changes in the relative negative-buoyancy of the slab. These have been thoroughly documented in the past [3] but need to be accounted for as a reference model. Second, there are changes in the subduction rate which occur when exotic material in ingested by the subduction zone due to colliding material influencing the plate coupling and the buoyancy force balance. Third, congestion generally influences the local retreat rate which leads to a time dependence in the shape of the subduction zone to accomodate differences in the rate of slab subduction along a uniformly converging margin.

Model

Our models are built upon Model 80w of [1], substituting a buoyant, weak, circular plateau for the continental ribbon material (See Figure 1). We present four models in which we vary the radius of the colliding plateau (see Table 1) from 0 to 500 km. The plateau is represented by a 30km deep cylinder of weaker (constant viscosity of 200 times the reference value), material with the density of basalt.

We use the particle-in-cell finite element code Underworld [4] which has been thoroughly tested for this class of problems [5]. The oceanic lithosphere is assumed to have cooled and to be negatively buoyant relative to the mantle according to a half-space cooling relationship. Temperature is mapped to density variations as a persistent property of a given material point during the simulation. The overriding plate parameterisation and geometry are identical to that study: 1) a region of weaker back-arc crust and mantle lithosphere, 2) a strong, cratonic backstop, and 3) a transitional region of crust and mantle lithosphere with a strength between these two. The asthenosphere viscosity is set to 3 times the reference value at the base of the lithosphere to account for the modest increase with pressure over the depth of the upper mantle ([2]).

Figure 1: (a) The initial distribution of distinct material domains that represent variations in density and strength. The rheological layering is as described in . The density and strength structure of the subducting plate is computed assuming a half-space cooling model to map lithospheric age to temperature and then averaging properties within 25 km thick layers: (b) yield strength in each layer of 80 Myr old ocean based on the lithostatic pressure gradient, (c) density of “normal” oceanic lithosphere with a 7km basaltic crust, (d) relative viscosity of “normal” oceanic lithosphere assuming the wet-olivine diffusion creep mechanism of  for all materials, (e) density profile through the plateau. In the viscosity plot, the geotherm is also shown (green line) and viscosity is truncated between \(10^{-2}\\) and \(10^5\\)
Figure 1: (a) The initial distribution of distinct material domains that represent variations in density and strength. The rheological layering is as described in [1]. The density and strength structure of the subducting plate is computed assuming a half-space cooling model to map lithospheric age to temperature and then averaging properties within 25 km thick layers: (b) yield strength in each layer of 80 Myr old ocean based on the lithostatic pressure gradient, (c) density of “normal” oceanic lithosphere with a 7km basaltic crust, (d) relative viscosity of “normal” oceanic lithosphere assuming the wet-olivine diffusion creep mechanism of [6] for all materials, (e) density profile through the plateau. In the viscosity plot, the geotherm is also shown (green line) and viscosity is truncated between \(10^{-2}\) and \(10^5\)

Notes on the Visualisations

The movies were rendered by the gLucifer system [7] which is part of the Underworld software suite. gLucifer was designed to allow efficient transport of particle visualization data from remote, parallel supercomputers. To show the location of the two plates, we use “sheets” of passively moving particles to represent, separately, the surface layer of the downgoing and over-riding plates. Each sheet carries a grid of strain markers which we rendered in a contrasting shade (See Figure 2c).

Tracers constrained to lie at the surface at the initial trench location (red) allow us to track the zone of convergence prior to the collision. The location of these tracers is used to compute the advance and retreat rates at different points along the subduction zone. Additional tracers were used to indicate the edges of the slab and these can be used to follow the re-initiation of subduction after the collision.

The movies themselves are made from snapshots taken at every 10 timesteps in the simulation and rendered at 15 frames per second. Elapsed time is approximately 1 Myr per frame but the timesteps of the simulation may vary due to numerical stability requirements.

Results

Model Number Movie Name Description
1 Plateau000.mov Reference model with no plateau material. Camera angle along the line of the trench
2 Plateau125.mov As Model 1 with Plateau radius 125km
3 Plateau250.mov As Model 1 with Plateau radius 250km
3a Plateau250i.mov Model 3 rendered oblique to the convergence direction
4 Plateau500.mov As Model 1 with Plateau radius 500km

Table 1: Model numbers, Movie file names and a description of the model parameters and rendering parameters

Model 1 is the reference model which has no plateau material. There is an inherent time dependence associated with the interaction of the subducting slab and the viscous material of the lower mantle (the lowermost 150km of the box). The evolution of the trench is largely two dimensional with approximately one third of the subduction due to trench retreat and two thirds due to the forward plate motion. We assume that the velocity of the plate and the trench retreat that we obtain from this simulation can be subtracted from the equivalent measurements from the other models to determine the influence of the collision.

Models 2, 3 and 4 have increasingly larger plateaus (doubling radius each time, i.e. 4 and 16 times the volume of buoyant material). In these models, we compare the behaviour of the colliding part of the margin with the “undisturbed” part, bearing in mind that the reference model described above has a time dependence of its own.

Figure 2a shows the velocity of the downgoing plate through time for each of the models in the set of curves indicated by the number 1. The envelope of all the curves is shaded to indicate the range of variability which occurs as a result of the plateau collision. The curves marked by the number 2 show the trench retreat rate away from the collision: there is less than 1cm/yr variation in this observable between any of the models, a significantly smaller difference than in the overall convergence rate. The curves indicated by the number 3 and plotted with solid circles indicate the velocity of the congested part of the trench. At the onset of the collision, retreating of the trench rapidly switches to advancing at close to the overall convergence rate.

Figure 2b focusses on Model 2 and shows the total plate subduction rate (trench retreat plus trenchward plate motion) for the “undisturbed” end of the trench in comparison with the collision. At 39Myr, when the collision starts, deformation in the over-riding plate accomodates the motion and the overall convergence between the downgoing plate and the suture between the plateau and the over-riding plate is less than 1cm/yr (about 1/3 the convergence rate for the retreating segment). In the animations, a zone of distributed deformation in the oceanic lithosphere behind the plateau begins can clearly be seen developing before a sudden re-initiation of the trench. This corresponds to the increase in convergence rate at 60Myr in Figure 2b.

Figure 2: (a) The velocity of the subducting plate, roll back of the trench (at \( z=0 \\) ), and the collisional suture (at \( z=z_\textrm{max}\\)) for each of the four models from the initiation of the model to a scaled run time of 85 Myr. The shading indicates the range of variation between the different models for the undisturbed parts of the system. (b) Model 2: a comparison of the plate velocity in the reference frame of the computational domain to the net subduction rate at either end of the trench. Roll back of the slab contributes to higher flux of subducted material than results from the plate movement alone (light shading). During collision, the trench advances at roughly the same rate as the plate velocity producing a period from 45-60Myr with no subduction in that region (heavily shaded). (c) Diagram showing a snapshot of model 3a with the plateau plus the slab and over-riding plate tracer-sheets visible. The numbers correspond to sample points that were used to measure the velocities for each model in sub-panel (a)
Figure 2: (a) The velocity of the subducting plate, roll back of the trench (at \( z=0 \) ), and the collisional suture (at \( z=z_\textrm{max}\)) for each of the four models from the initiation of the model to a scaled run time of 85 Myr. The shading indicates the range of variation between the different models for the “undisturbed” parts of the system. (b) Model 2: a comparison of the plate velocity in the reference frame of the computational domain to the net subduction rate at either end of the trench. Roll back of the slab contributes to higher flux of subducted material than results from the plate movement alone (light shading). During collision, the trench advances at roughly the same rate as the plate velocity producing a period from 45–60Myr with no subduction in that region (heavily shaded). (c) Diagram showing a snapshot of model 3a with the plateau plus the slab and over-riding plate tracer-sheets visible. The numbers correspond to sample points that were used to measure the velocities for each model in sub-panel (a)

Summary

Van Hunen and Miller, in this volume, have given an overview of the time-dependent behaviour of subduction during ocean closure and collision. Here we have presented a short series of simulations which show some of the natural time dependence of "two dimensional" subduction and the three-dimensional signatures of small-scale collisions in the deformation record of the over-riding plate.

These movies explore the influence of plateau size on the dynamic evolution of the slab. The larger collisions produce the more dramatic signatures of shortening and re-configuration of the over-riding plate and reduce the speed of the downgoing plate (Figure 2a). We observe, however, that even in this case, the overall rollback rate for the uncongested section of the trench is only weakly affected.

We observe local changes in the trench geometry by which the subduction zone can accommodate a combination of collision and trench retreat along its length. The global subduction rate is not much affected by this change (Figure 2b) and this provides a continuing source of negative buoyancy to drive compressional deformation in the over-riding plate. Continued subduction in the retreating section of the margin also provides a compressive stress regime in the downgoing plate behind the plateau which leads to the eventual re-initiation of the subduction zone there.

References









  1. Moresi, L., Betts, P. G., Miller, M. S., & Cayley, R. A. (2014). Dynamics of continental accretion. Nature. doi:10.1038/nature13033

  2. Betts, P. G., Moresi, L., Willis, D., Miller, M. S. (2014). Geodynamics of oceanic plateau and plume head accretion and their role in Phanerozoic orogenic systems of China (Geoscience Frontiers, In Press, July 24, 2014)

  3. Stegman, D. R., Farrington, R., Capitanio, F. A., & Schellart, W. P. (2010). A regime diagram for subduction styles from 3-D numerical models of free subduction. Tectonophysics.

  4. Moresi, L., Quenette, S., Lemiale, V., Mériaux, C., Appelbe, W., Mühlhaus. (2007). Computational approaches to studying non-linear dynamics of the crust and mantle. Phys. Earth Planet. Inter., 163, 69–82. doi:10.1016/j.pepi.2007.06.009

  5. Stegman, D. R., Moresi, L., Freeman, J., Schellart, W. P., & May, D. A. (2006). Influence of trench width on subduction hinge retreat rates in 3D models of slab rollback. Geochem. Geophys. Geosys., 7, Q03012.

  6. Karato, S.-I., & Wu, P. (1993). Rheology of the Upper Mantle: A Synthesis. Science, 260(5109), 771–778. doi:10.1126/science.260.5109.771

  7. Stegman, D. R., Moresi, L., Turnbull, R., Giordani, J., Sunter, P., Lo, A., & Quenette, S. (2008). gLucifer: next generation visualization framework for high-performance computational geodynamics. Visual Geosciences, 13(1), 71–84. doi:10.1007/s10069–008–0010–2