of 8
Geophysical Research Letters
Subduction Initiation With Vertical Lithospheric
Heterogeneities and New Fault Formation
Xiaolin Mao
1
, Michael Gurnis
1
, and Dave A. May
2
1
Seismological Laboratory, California Institute of Technology, Pasadena, CA, USA,
2
Department of Earth Sciences,
University of Oxford, Oxford, UK
Abstract
How subduction initiates with mechanically unfavorable lithospheric heterogeneities is
important and rarely studied. We investigate this with a geodynamic 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 Puysegur
Trench and Ridge and the deformation history on the overriding plate. We show how a new thrust fault
forms and evolves into a smooth subduction 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 that 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 accumulates beneath the
Snares Zone.
1. Introduction
Subduction initiation is a vital phase of the plate tectonic cycle since it fundamentally alters the global
force balance on tectonic plates. Numerical studies have advanced our understanding of under what cir-
cumstances and with what physical processes a new subduction zone can develop (e.g., Gurnis et al., 2004;
Nikolaeva et al., 2010; Regenauer-Lieb et al., 2001; Thielmann & Kaus, 2012; Toth & Gurnis, 1998), but uncer-
tainty 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 sponta-
neous: 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 subduction 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 overprinted 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 sub-
duction zones are so young (Gorringe Bank, the Owen Ridge, the Hjort Trench, and Mussau Trench) that the
slab may not have yet started to bend into the mantle (Gurnis et al., 2004).
Luckily, one subduction zone overcomes these limitations: the Puysegur Incipient Subduction 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 (Figure 1a). Since about 20 Ma, highly oblique con-
vergence 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 (Figure 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
& Parkinson, 1997). This confirms that the slab enters the mantle. The morphology of the Puysegur Ridge
(Figure 1b) shows a characteristic change from uplift in the southern part, where the total convergence is
less, and subsidence in the northern part, the Snares Zone (or Snares Trough), where the total convergence
RESEARCH LETTER
10.1002/2017GL075389
Key Points:
• New fault forms when subduction
initiates with unfavorable dipping
lithospheric 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
Supporting Information:
• Supporting Information S1
Correspondence to:
X. Mao,
xlmao@gps.caltech.edu
Citation:
Mao, X., Gurnis, M., & May, D. A. (2017).
Subduction initiation with vertical
lithospheric heterogeneities and
new fault formation.
Geophysical
Research Letters
,
44
, 11,349–11,356.
https://doi.org/10.1002/2017GL075389
Received 22 AUG 2017
Accepted 24 OCT 2017
Accepted article online 30 OCT 2017
Published online 20 NOV 2017
©2017. American Geophysical Union.
All Rights Reserved.
MAO ET AL.
MAO ET AL.: SUBDUCTION INITIATION
11
,
349
Geophysical Research Letters
10.1002/2017GL075389
Figure 1.
(a) Tectonic outline of the Puysegur region. AUS: Australian Plate; PAC: Pacific Plate; MRC: Macquarie Ridge
Complex. (b) Bathymetry of the Puysegur Incipient Subduction Zone. The red vector is the relative velocity of AUS to
PAC (DeMets et al., 1994). Stars show the location of young volcanic features (Sutherland et al., 2006), and black lines
show position of cross sections in Figure 1d. Dashed lines show possible position of the Puysegur Fault. (c) Filled circles
are earthquakes with magnitude between 2 and 5 from the ISC Bulletin (International Seismological Centre, 2014), and
the focal mechanism solutions are from CMT catalog (Dziewonski et al., 1981). (d) Bathymetry cross sections and inferred
fault structures (modified from Lebrun et al., 1998).
is largest and is roughly consistent with geodynamic models 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 (Figure 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 uplift, while there was uplift followed by
subsidence in the Snares Zone (Collot et al., 1995; Gurnis et al., 2004; Lebrun et al., 1998). The width of the
Puysegur ridge also widens northward from less than 50 km at 49.5
Sto
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 (Figure 1d). Together, they 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 & 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 2-D 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 earthquakes and interpretations of multibeam bathy-
metric, sonar imagery, seismic reflection, and geopotential data (Collot et al., 1995; Ruff et al., 1989). 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 pre-
existing 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 initiation (Gurnis et al., 2004; Leng & Gurnis, 2011). However, it is not clear
MAO ET AL.
MAO ET AL.: SUBDUCTION INITIATION
11
,
350
Geophysical Research Letters
10.1002/2017GL075389
Figure 2.
Model setup. The gray area shows the shape of the weak zone, with dip angle,
, and width,
d
wz
.Theweak
zone is initialized with random plastic strain within 0–0.4, while other material properties remain unchanged. The finest
resolution is 1 km
×
1 km near the subduction zone, and the lowest is 3 km
×
3 km in the asthenosphere.
whether PISZ is self-sustaining, or if the along strike variation in the uplift and subsidence of the Puysegur
Ridge represents this transition. Here we use 2-D 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., Billen & Gurnis,
2001; Kaus et al., 2008; Gerya & Meilick, 2011; Zhong & Gurnis, 1994), predicting reliable topographic evolu-
tion 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 impor-
tant geophysical, petrological, and geochemical processes that affect the force balance, and the model
setup, boundary conditions, and the evolution of material properties need to be self-consistent and prop-
erly 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 & Gurnis, 1998). With the subduction of the
downgoing slab, phase changes in the crust lead to an increasing crustal density, and this part of the resist-
ing force evolves to drive subduction. Fluids released from the subducting 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).
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 topo-
graphic evolution. Initial topography is calculated from isostasy (Figure 2), and topography is updated
between time steps with surface velocity under the constraint that the vertical topographic change is smaller
than 20 m to avoid topographic oscillation (Kaus et al., 2010). A simplified surface process model, based on
linear topographic diffusion, is implemented (e.g., Avouac, 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 4-D (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 & Kohlstedt, 1996).
The Drucker-Prager yield criteria with a maximum yield stress are employed for material plasticity, and the
accumulated plastic strain is recorded on tracers and used to
reduce the material friction coefficient and cohe-
sion. 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
in the supporting information. More details on the numerical implementations of surface process, phase
change, and EVSS are found in the supporting information.
MAO ET AL.
MAO ET AL.: SUBDUCTION INITIATION
11
,
351
Geophysical Research Letters
10.1002/2017GL075389
Figure 3.
Model results (case SI1). (a) Effective viscosity evolution. (b) Accumulation of plastic strain. The black box
shows the corresponding region for Figure 4a. (c) Density evolution. Black lines show direction and magnitude of
maximum principal stress. Rock types and free water contents are shown in the insets with different colors and white
contours. (d) Topography changes. Blue and red lines show initial and final topography for each time interval.
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 (Figure 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 Figure 4a).
Displacement on this new fault is small, and wide spread shear bands together with buckling of the litho-
sphere 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 (Figure 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 (Figure 3d). Topographic fluctua-
tion on the order of a hundred meters is predicted due to the lithospheric buckling in the overriding plate
(Figure 3d). By 3.87 Myr, the shear band on the subducting plate side evolves into a new fault which later
becomes the subduction interface. An
30 km wide triangular lithospheric block which was on the subduct-
ing plate attaches to the overriding plate with the formation of the new plate boundary. A two-fault system,
MAO ET AL.
MAO ET AL.: SUBDUCTION INITIATION
11
,
352
Geophysical Research Letters
10.1002/2017GL075389
Figure 4.
Model summary. (a, b) The zoomed-in plastic strains for cases SI1 and SI2, respectively. WM: weak material.
(c) Averaged horizontal stresses within overriding plate for different cases. The averaging is performed within the region
that spans from 100 km to 300 km distance and from surface to 25 km depth.
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 manifested 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
smoothing,thecompressivestressintheoverridingplatedecreases(Figures3cand4c),andshearbandsinthe
overriding plate become inactive (Figure 3a). Topography of the overriding plate close to the trench changes
from uplift to subsidence, and the trench depth deepens (Figure 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 litho-
sphere. Interestingly, these normal faults become temporarily inactive after they pass the bending region and
enter the subduction zone (Figure 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 man-
tle 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 (Figure 3d),
and weak material consisting of sediments and crustal and lithospheric components (scrapped off along the
interface) accumulates beneath the trough (Figure 4a).
MAO ET AL.
MAO ET AL.: SUBDUCTION INITIATION
11
,
353