Subduction initiation with vertical lithospheric
heterogeneities and new fault formation
Xiaolin Mao
1
, Michael Gurnis
1
, Dave A. May
2
Xiaolin Mao, xlmao@gps.caltech.edu
1
Seismological Laboratory, California
Institute of Technology, Pasadena, CA,
91125, USA
2
Department of Earth Sciences, University
of Oxford, Oxford, United Kingdom
This article has been accepted for publication and undergone full peer review but has not been through
the copyediting, typesetting, pagination and proofreading process, which may lead to differences be-
tween this version and the Version of Record. Please cite this article as doi: 10.1002/2017GL075389
c
©
2017 American Geophysical Union. All Rights Reserved.
Abstract.
How subduction initiates with mechanically unfavorable lithospheric het-
erogeneities is important and rarely studied. We investigate this with a geo-
dynamic model for the Puysegur Incipient Subduction Zone (PISZ) south
of New Zealand. The model incorporates a true free surface, elasto-visco-plastic
rheology and phase changes. Our predictions fit the morphology of the Puy-
segur Trench and Ridge and the deformation history on the overriding plate.
We show how a new thrust fault forms and evolves into a smooth subduc-
tion interface, and how a preexisting weak zone can become a vertical fault
inboard of the thrust fault during subduction initiation, consistent with two-
fault system at PISZ. The model suggests that the PISZ may not yet be self-
sustaining. We propose the Snares Zone (or Snares Trough) is caused by plate
coupling differences between shallower and deeper parts, the tectonic sliver
between two faults experiences strong rotation, and low density material ac-
cumulates beneath the Snares Zone.
Keypoints:
•
New fault forms when subduction initiates with unfavorable dipping litho-
spheric heterogeneities
•
Uplift and subsidence on the Puysegur Ridge are related to new fault
formation and maturation
•
Snares Zone is caused by differences in plate coupling at shallower and
deeper fault depths
c
©
2017 American Geophysical Union. All Rights Reserved.
1. Introduction
Subduction initiation is a vital phase of the plate tectonic cycle since it fundamen-
tally alters the global force balance on tectonic plates. Numerical studies have advanced
our understanding of under what circumstances and with what physical processes a new
subduction zone can develop [e.g.,
Toth and Gurnis
, 1998;
Regenauer-Lieb et al.
, 2001;
Gurnis et al.
, 2004;
Nikolaeva et al.
, 2010;
Thielmann and Kaus
, 2012], but uncertainty
still exists and key parameters related to subduction initiation remain poorly quantified
mainly due to the lack of good constraints on numerical models. Subduction initiation
can be either induced or spontaneous: induced subduction initiation begins with strong
compression and uplift [
Gurnis et al.
, 2004], whereas spontaneous initiation begins with
rifting and subsidence [
Stern
, 2004]. Therefore, topographic changes that result from sub-
duction initiation can be used to distinguish different initiation modes and can potentially
be used to quantify parameters that control the initiation process. While the importance
of topographic change in subduction initiation has been noticed previously [
Gurnis et al.
,
2004], applying topographic changes as a constraint on subduction initiation process at
a specific subduction zone has not been addressed, as most of the early record needed
to constrain the dynamics is over-printed by later deformation and volcanism for mature
subduction zones, like the well known examples of the Eocene initiation of Izu-Bonin-
Marianas and Tonga-Kermadec subduction zones [
Sutherland et al.
, 2006]. On the other
hand, some possible incipient subduction zones are so young (Gorringe Bank, the Owen
Ridge, the Hjort Trench, Mussau Trench) that the slab may not have yet started to bend
into the mantle [
Gurnis et al.
, 2004].
c
©
2017 American Geophysical Union. All Rights Reserved.
Luckily, one subduction zone overcomes these limitations: the Puysegur Incipient Sub-
duction Zone (PISZ). The Puysegur Trench and Ridge form the northern end of the
Macquarie Ridge Complex (MRC) defining the Australian-Pacific plate margin south of
New Zealand (Fig. 1a). Since about 20 Ma, highly oblique convergence beneath the
Puysegur Ridge results in a maximum total convergence of 150-200 km at Puysegur as
suggested by a Benioff zone with seismicity down to 150 km depth (Fig. 1c) [
Sutherland
et al.
, 2006]. Subduction related igneous rocks, especially adakite, which is formed by
the partial melting of young oceanic crust under eclogitic facies conditions, are sparsely
distributed on the overriding plate at Solander Island [
Reay and Parkinson
, 1997]. This
confirms that the slab enters the mantle. The morphology of the Puysegur Ridge (Fig.
1b) shows a characteristic change from uplift in the southern part, where the total con-
vergence is less, and subsidence in the northern part, the Snares Zone (or Snares Trough),
where the total convergence is largest, and is roughly consistent with geodynamic mod-
els of induced subduction initiation [
Gurnis et al.
, 2004]. Discrete flat-topped segments,
which are interpreted as the results of subaerial exposure and erosion, are also evident at
both northern and southern parts of the Puysegur Ridge (Fig. 1d). The southernmost
segment is close to sea level (-120 m), while a peak subsidence of
∼
1800 m is found in
the Snares Zone in the north, suggesting that the southern part has only experienced up-
lift, while there was uplift followed by subsidence in the Snares Zone [
Collot et al.
, 1995;
Lebrun et al.
, 1998;
Gurnis et al.
, 2004]. The width of the Puysegur ridge also widens
northwards from less than 50 km at 49.5
◦
S to
∼
80 km at 47.5
◦
S . A confined strike
slip fault zone is found near the peak of the ridge in the south, while a splayed fault
zone structure is suggested in the trough of the Snares Zone (Fig. 1d). Together, they
c
©
2017 American Geophysical Union. All Rights Reserved.
show that the overriding plate close to the trench is under compression in the south and
potentially in extension in the north [
Collot et al.
, 1995;
Lamarche and Lebrun
, 2000].
The corresponding trench depth in the south is about 1 km shallower than in the north
[
Collot et al.
, 1995]. These spatial variations in structure along the PISZ are thought to
represent different time periods in the evolution, and a space-for-time substitution can be
made to compare the time evolution from 2D models with the spatial variation along the
PISZ.
The two-fault system at PISZ, with a thrust fault at the trench and a vertical fault
inboard of the thrust fault, is recognized through the dual rupture mode for large earth-
quakes and interpretations of multibeam bathymetric, sonar imagery, seismic reflection
and geopotential data [
Ruff et al.
, 1989;
Collot et al.
, 1995]. However, how this two-fault
structure formed is still open to debate.
Ruff et al.
[1989] propose that the thrust interface
formed through the propagation and connection of disconnected small thrust faults behind
the vertical fault, while this hypothesis preceded the mapping of these two faults.
Collot
et al.
[1995] propose that it developed through progressive adjustments of two adjacent
vertical weak zones under compression, which requires the development of the ridge and
differential uplift of the crustal block at one weak zone, and the rotation of the other.
Here, we propose that the two-fault system formed through developing a new thrust fault
near the preexisting vertical fault during subduction initiation. Previous studies suggest
that there is a transition in the force balances from being forced externally to a state of
self-sustaining subduction under its own negative buoyancy for induced subduction initi-
ation [
Gurnis et al.
, 2004;
Leng and Gurnis
, 2011]. However, it is not clear whether PISZ
is self-sustaining, or if the along strike variation in the uplift and subsidence of the Puy-
c
©
2017 American Geophysical Union. All Rights Reserved.
segur Ridge represents this transition. Here we use 2D geodynamic models, which have a
true free surface to track topographic changes, and model setup and boundary conditions
tailored for PISZ, to test our hypothesis for the formation of the two-fault system at PISZ,
and to explore the factors that control the transition in the force balance, while focusing
on the evolution of topography and state of stress.
2. Method
Although a number of studies have tracked the topography in subduction zone models
[e.g.,
Zhong and Gurnis
, 1994;
Billen and Gurnis
, 2001;
Kaus et al.
, 2008;
Gerya and
Meilick
, 2011], predicting reliable topographic evolution in subduction zone remains a
challenge. By reliable topography, we mean that not only is the predicted topography
consistent with that observed, but the model should also be able to include most of the
important geophysical, petrological and geochemical processes that affect the force bal-
ance, and the model setup, boundary conditions and the evolution of material properties
need to be self-consistent and properly constrained. During subduction initiation, the
driving forces must overcome the resisting forces from the friction on the sliding interface,
the bending of lithosphere and the buoyancy of oceanic crust (before the basalt-to-eclogite
phase change) [e.g.,
McKenzie
, 1977;
Toth and Gurnis
, 1998]. With the subduction of the
down-going slab, phase changes in the crust lead to an increasing crustal density, and this
part of the resisting force evolves to drive subduction. Fluids released from the subduct-
ing crust contribute to the decoupling between subducting and overriding plates, which
may reduce the resisting force. Also influencing the force balance is the development of
topography and surface process, which generate loads that affect lithospheric and mantle
dynamics [
Kaus et al.
, 2008, 2010].
c
©
2017 American Geophysical Union. All Rights Reserved.
A true free surface is tracked in pTatin3D [
May et al.
, 2014, 2015], based on the Arbitrary
Lagrangian Eulerian (ALE) finite element method, and is used to follow the dynamic
mantle-surface interactions and the topographic evolution. Initial topography is calculated
from isostasy (Fig. 2), and topography is updated between time steps with surface velocity
under the constraint that the vertical topographic change is smaller than 20 meters to
avoid topographic oscillation [
Kaus et al.
, 2010]. A simplified surface process model, based
on linear topographic diffusion, is implemented [e.g.,
Avouac and Burov
, 1996]. A 5 km
thick altered basaltic crust is placed on the top of dry pyrolite, and sediments are generated
at the surface with our surface process. Density and free water content for different phase
assemblages are gained by referring to precalculated 4D (temperature, pressure, rock type
and total water content) phase maps using Perplex [
Connolly
, 2005]. Darcy’s law is used
to migrate free water, and a linear water weakening is applied to the mantle material
[
Hirth and Kohlstedt
, 1996]. The Drucker-Prager yield criteria with a maximum yield
stress is employed for material plasticity, and the accumulated plastic strain is recorded
on tracers and used to reduce the material friction coefficient and cohesion. Elasticity
alters the stress within the slab, and we include a new visco-elastic formulation called
the Elastic Viscous Stress Splitting (EVSS) method [e.g.,
Keunings
, 2000] in our model.
Energy change with shear heating is treated as heat source terms and stored on tracers.
Thermal and rheological parameters are given in Table S1. More details on the numerical
implementations of surface process, phase change, and EVSS are found in supplementary
material.
c
©
2017 American Geophysical Union. All Rights Reserved.
3. Results
We first show the detailed evolution of a subduction initiation case (SI1) that starts
with a vertical weak zone with imposed far field plate velocities (Fig. 3), which are key
aspects of PISZ. Initially (by 2.26 Myr), both the overriding and subducting plates are in
a strong sense of compression, and the first lithospheric scale fault forms at the top of the
initial weak zone with a high dip angle (fine structures shown in Fig. 4a). Displacement
on this new fault is small, and wide spread shear bands together with buckling of the
lithosphere develop within the overriding plate. The buckling and shear bands absorb a
considerable amount of the convergence. Some strain begins to localize at a shear band on
the subducting plate side, and this shear band connects to the deeper part of the initial
weak zone (Fig. 4a). Responding to strong compression, a topographic pair of uplift
adjacent to subsidence appear near the plate boundary at the overriding and subducting
sides respectively, and amplitudes of both are around 1.5 km (Fig. 3d). Topographic
fluctuation on the order of a hundred meters is predicted due to the lithospheric buckling
in the overriding plate (Fig. 3d). By 3.87 Myr, the shear band on the subducting plate
side evolves into a new fault which later becomes the subduction interface. A
∼
30
km wide triangular lithospheric block which was on the subducting plate attaches to
the overriding plate with the formation of the new plate boundary. A two-fault system,
with the main thrust fault at the subduction interface and a vertical fault inboard of
the thrust fault, starts to appear and dominate the deformation and structure of the
overriding plate close to the trench. A sharp topographic signal is manifest with the
uplift and rotation of the triangular block. The rough subduction interface is slowly
smoothed by continuous tectonic erosion and plate bending from 3.87 to 11.65 Myr. With
c
©
2017 American Geophysical Union. All Rights Reserved.
smoothing, the compressive stress in the overriding plate decreases (Fig. 3c and Fig. 4c),
and shear bands in the overriding plate become inactive (Fig. 3a). Topography of the
overriding plate close to the trench changes from uplift to subsidence, and the trench depth
deepens (Fig. 3d). The bending of the subducting plate causes near surface extension
and compression deeper within the subducting lithosphere near the trench, which leads
to development of the forebulge topography and normal faults in the upper part of the
lithosphere. Interestingly, these normal faults become temporarily inactive after they pass
the bending region and enter the subduction zone (Fig. 3a). Basalt transitions to eclogite
at 70-80 km depth resulting in a density jump in the oceanic crust, and decoupling occurs
between the overriding and subducting plates in the mantle wedge caused by released
fluids. Both tend to bend the subducted slab to a higher dipping angle. A trough
forms at the overriding plate close to the trench region when the slab reaches
∼
200 km
depth (Fig. 3d), and weak material consisting of sediments and crustal and lithospheric
components (scrapped off along the interface) accumulate beneath the trough (Fig. 4a).
We compare this case with one with a 30
◦
dip angle weak zone (Fig. S1). Initially,
the subduction interface develops easily within the weak zone, and strain localizes to this
thrust fault efficiently. The shape of the fault is favorable for subduction, and compression
within the overriding plate decreases much faster compared to the case with a vertical
weak zone. Only one fault develops during the entire subduction initiation process, and
few shear bands appear during the early stage of compression. Lithospheric buckling of the
overriding plate is absent in this case. The overriding plate close to the trench experiences
uplift and subsidence during the initiation, while no trough-like structure develops. Unlike
the vertical weak zone case, where the whole overriding plate is under weak compression
c
©
2017 American Geophysical Union. All Rights Reserved.
or at a nearly neutral stress state after the smoothing of the subduction interface, the
stress state in the overriding plate transitions into the extension phase quickly (Fig. 4c).
To fill the gap between the two end members of shallow and steep dipping weak zones,
we computed three additional cases with dip angles of 45
◦
, 60
◦
and 75
◦
(Fig. S2-S4).
These models can be divided into two categories: one starts with high dip angles and
evolves to a two-fault system (Fig. 4a), while the other starts with a shallow dip and
converges to a single-fault system (Fig. 4b). During initiation, resisting forces in the
two-fault system are at higher stress level, and the compression in the overriding plate
can last much longer compared to the single-fault system (Fig. 4c).
To test the other geometrical parameter that may influence the fault structure during
subduction initiation, we conduct two more cases varying the width of the weak zone with
a 90
◦
dip angle (Fig. S5-S6). With a wider weak zone, the evolution pathway is generally
similar to case SI1. The main differences are the triangular block between the thrust and
vertical faults becomes smaller and the compressive stress within the overriding plate is
at a slightly lower level. With half the width of the initial weak zone, the model fails to
localize the deformation near the plate boundary. These two models confirm that once
the weak zone is wide enough to localize the initial deformation, the dip angle plays a
more important role to determine if a new fault is needed to initiate subduction.
In the previous cases, we focus on the influence of weak zone geometry on fault structure
and stress state within the lithosphere during subduction initiation under compression,
and we ignore a harzburgite layer in our density structure for simplicity (Fig. S7). How-
ever, harzburgite has a slightly lower density compared to the undepleted mantle, and
could increase resisting forces for slab subduction [e.g.,
Oxburgh and Parmentier
, 1977;
c
©
2017 American Geophysical Union. All Rights Reserved.