Earth and Planetary Science Letters 353354 (2012) 2937

Contents lists available at SciVerse ScienceDirect

Earth and Planetary Science Letters

journal homepage: www.elsevier.com/locate/epsl

Complex mantle ow around heterogeneous subducting oceanic plates

Fabio A. Capitanio n, Manuele Faccenda 1
School of Geosciences, Monash University, Clayton, 3800 VIC, Australia

a r t i c l e i n f o

abstract

Article history:
Received 12 May 2012
Received in revised form
27 July 2012
Accepted 28 July 2012
Editor: Y. Ricard
Available online 5 September 2012

The foundering of oceanic lithospheres controls the circulation patterns of the mantle around
subducting slabs. Here, we investigate the sensitivity of the mantle ow to slab buoyancy variations
along convergent margins using three-dimensional numerical models of subduction in a viscous
mantle. The models illustrate that in a buoyancy-driven system varying subduction velocity arising
from negative buoyancy variations effectively drives pressure gradients and confers a general ow
component sub-parallel to the margins strike, allowing for material transport over large distances
around the slabs. The along-slab velocity component introduces widespread horizontal simple shear in
the mantle ow which is maximized beneath the slab between  100 and  350 km. Mantle ow
complexities develop rapidly, although not instantaneously, upon subduction of heterogeneous plates.
The resulting slab pull gradients are mostly accommodated by internal slab deformation, decoupling
the mantle ow from motions at surface. Moderate slab pull variations have a minor impact on plate
velocity, and might not result in trench motions, although effectively rearrange the ow. Slab buoyancy
heterogeneities are rstly associated with age-dependent thickness variations, but also with slab break
offs and windows, varying depth of subduction and entrainment of buoyant blocks. Because these are
observed at all subduction zones, the process shown here should have a global relevance for the ow
around slabs.
& 2012 Elsevier B.V. All rights reserved.

Keywords:
subduction modeling
mante ow
trench-parallel ow
South American subduction

1. Introduction
The negative buoyancy of sinking slabs fundamentally controls
the circulation of the mantle around subducting lithosphere. Yet,
the complexities of mantle ow and mass transport on Earth
remain difcult to constrain from the available indirect observation. Simplied geodynamic models of subduction show that the
mantle is entrained with subducting slabs, owing perpendicular
to the trench in both the subslab and mantle wedge, and mass
exchange is conned to the slab edges where the ow is mostly
toroidal (Funiciello et al., 2004, 2006; Kincaid and Grifths, 2003;
Morra and Regenauer-Lieb, 2006; Schellart, 2004; Stegman et al.,
2006). However, the exposed geological analog of the Talkeetna
arc, Alaska, shows the occurrence of trench-parallel ow in
exhumed mantle wedge rocks, indicated by seismic polarization
on horizontal planes (Mehl et al., 2003). Furthermore, isotope
geochemistry, seismic velocity anisotropy and attenuation suggest that the ow direction is arc-parallel in the mantle wedge
beneath Costa Rica and Nicaragua (Hoernle et al., 2008).

Corresponding author. Tel.: 61 3 99053160.

E-mail address: fabio.capitanio@monash.edu (F.A. Capitanio).
1
Present address: Department of Geosciences, University of Padua, 35131
Padua, Italy.
0012-821X/$ - see front matter & 2012 Elsevier B.V. All rights reserved.
http://dx.doi.org/10.1016/j.epsl.2012.07.042

Analytical solutions indicate that lateral gradients in slab

properties such as subduction velocity and dip induce similar
gradients in pressure (Turcotte and Schubert, 1982), under which
the mantle ow may acquire a trench-parallel component (Hall
et al., 2000; Kneller and van Keken, 2007). Models of subduction
that include realistic slab morphologies, velocity and slab dip
variations along the trench (Honda and Yoshida, 2005; Kneller
and van Keken, 2007; Morishige and Honda, 2011; Natarov and
Conrad, 2012) show some complex patterns of ow, although no
trench-parallel and toroidal ow that may reproduce the observations (Hoernle et al., 2008; Rychert et al., 2008) has been
found yet.
Observables of subduction velocities and dips at convergent
margins show that along-trench variations are common features
of subducting slabs (Lallemand et al., 2005), thus suggesting that
this mechanism might be of general relevance. In models of
buoyancy-driven subduction, variations of subduction rates and
dips are rstly related to the balance between the slabs own
buoyancy force and plate stiffness and mantle wedge suction
(Capitanio and Morra, 2012; Capitanio et al., 2007; RodrguesGonzalez et al., 2012). Similar buoyancy variations in the subducting lithosphere are primarily due to the age-dependent
thickening of plates, so that the observed variations of oceanic

plate age along trenches (Muller

et al., 2008) must correspond to
similar variations in thickness. If so, these should have a role on

30

F.A. Capitanio, M. Faccenda / Earth and Planetary Science Letters 353354 (2012) 2937

the distribution of driving buoyancy force along trenches on

Earth, and likely contribute to the organization of circulation
patterns at depth.
Here, we use buoyancy-driven subduction modeling to test
this mechanism in a subduction system. This approach allows
modeling subduction velocities, trench motions, slab dips and
mantle ow as emerging self-consistent features of the force
balance between slab buoyancy, mantle drag and upper plates
resistance (Capitanio et al., 2010b, 2011). Rather than focusing on
the role of single parameters, we draw implications for the largescale dynamics of the coupled lithosphere-mantle system. By
modeling heterogeneous buoyancy along the width of subducting
slabs, we introduce trench-parallel gradients in subduction rates,
trench migrations and dips along the trench and test their impact
on ow patterns. In these models, slab pull variations induce
sinking rates and pressure lateral gradients introducing a lateral
component of the mantle velocity eld, perturbing the convective
pattern driven by down-going slab. The lateral ow component
induces large trench-parallel simple shear strain in a layer at 100
350 km depth, beneath the slabs. Where the trench-ward ow
component vanishes, trench-parallel ow occurs, and mass can be
transported over large distances around wide slabs. Slabs in the
mantle and non-Newtonian mantle ow rearrange readily to
minor buoyancy perturbation, although this might not impact
plate and trench motions at surface.

2. Modeling strategy, method and setup

Aim of modeling in this work is to test the effect of oceanic
plates buoyancy variation on subduction. The buoyancy of the
plate is determined by the density contrast with the mantle
integrated over its volume (times gravity). The temperaturedependent density of mantle forming minerals has a primary
control on the buoyancy of oceanic lithospheres, modulating the
vertical density distribution as well as the thickness of the
thermal boundary layer. A detailed analysis of lithospheric mantle
density integrated over h  105 km yields average density contrast of Dr 54:3 kg m3 for a 50 Myr oceanic plate, varying to
28.7 kg m  3 for 25 Myr and 63.4 kg m  3 for oceanic lithospheres
of 90 Myr (Afonso et al., 2007). However, the average density
contrast is more likely constant on Earth,  40 kg m  3 (Afonso
et al., 2007), whereas the buoyancy varies with the age-dependent thickness of the plate, which scales as h 2:32kt1=2 , where
k is the thermal diffusivity (Turcotte and Schubert, 1982). For
plates of 25, 50, 90 Myr this yields thickness of 65, 92 and 123 km,
respectively. Variations of subducting oceanic lithosphere age

along the trench are commonly observed (Muller

et al., 2008)

varying over distances of O(103 km) (Sdrolias and Muller,

2006),
resulting in similar variation of thickness in the subducting slabs.
Additionally to this rst order buoyancy heterogeneity, oceanic
lithospheres host extinct ridges or oceanic plateaux, where lithosphere is locally depleted. Although the integrated buoyancy
could vary only by a few percent (Afonso et al., 2007) these can
be wide enough to affect the sudbuction, as the Ogasawara
Plateau (Miller et al., 2006). Another source of buoyancy heterogeneity during subduction is provided by the entrainment of
continental lithosphere segments or blocks. The density and
thickness of continental lithosphere vary largely from those of
mature oceanic plates, impacting subduction dynamics (Capitanio
et al., 2010a). In the subduction models presented here, we use
the integrated buoyancy of the plate as a free parameter and
discuss later the implications for Earth.
Subduction is modeled as the viscous ow of an innite
Prandtl number uid at very low Reynolds number in a threedimensional Cartesian geometry. The ow caused by internal

buoyancy is described by the mass and momentum conservation

equations:

r v 0

rpr s Drg

where p is the pressure, r the density, g the gravity vector, v the

velocity and t is the deviatoric stress tensor:


@vi @vj
tij Z

3
2Ze_ ij ,
@xj
@xi
where Z is the dynamic viscosity and e_ is the strain rate tensor.
These equations are solved in their non-dimensionalized form
with Underworld, a Lagrangian integration point nite element
method, in a 3D Eulerian grid (Moresi et al., 2003). The details of
the numerical method, software implementation and relevant
numerical benchmarks are described in Moresi et al. (2003) and
Stegman et al. (2006). The numerical resolution of the models is
96  96  64, in the x, y, z directions respectively.
The model space is a 3D Cartesian box, 4000 km long, 8000 km
wide and 1000 km deep, and is derived from Capitanio et al.
(2011). The boundary conditions are periodic on the front and
rear walls (x 0, 4000 km), free-slip on top and sidewalls, no-slip
on the bottom and a symmetry plane in y0 (Supplementary
Fig. 1). For the symmetry adopted only half model space is
modeled, i.e. y [0, 4000 km]. We neglect the role of the Earths
sphericity since it has a negligible role on the subduction
dynamics (Morra et al., 2006), and does not affect the pressure
gradients addressed. The subducting lithosphere is 6000 km wide
and is a Newtonian body, with a strong, 1023 Pa s, 70 km thick top
layer, and a 1021 Pa s, 30 km thick bottom layer, with an additional Byerlees plasticity law in the top 30 km:

tY t0 mp,

where t0 35 MPa is the cohesion, m 0.01 the friction coefcient,

in agreement with commonly used values (Krien and Fleitout,
2008; Sobolev and Babeyko, 2005). The subducting plate-mantle
density contrast chosen is Dr 53 kg m, 3 compatible with realistic buoyancy estimates (Afonso et al., 2007). The models upper
plate has no density contrast with the mantle, a thickness of
40 km and Newtonian viscosity of 1024 Pa s, reducing the trench
curvature to a minimum (Supplementary Table 1). An incipient
slab, extending to 150 km depth initiates subduction, after which
it is self-sustaining.
The upper mantle viscosity is
 1n=n
e_
ZM Z0 _ II
,
5
e II0
p
where e_ II e_ ij e_ ij =2 is the second invariant of the strain rate
tensor, and e_ II0 the reference strain rate. Every model has been run
with values of n 1 and Z0 2  1020 Pa s, i.e. Newtonian mantle,
and n 3.5 and Z0 2  1018 Pa s, i.e. non-Newtonian mantle
(Suppl. Table 1). A linear viscosity of 1021 Pa s is imposed below
660 km depth, and a density increase of 50 kg m  3 (Suppl. Fig. 2).
The models run to steady state subduction of a reference
model with constant buoyancy slab and Newtonian mantle (that
is more computationally efcient), then the model is perturbed.
The models enter a steady state after the slab has reached the
transition zone, i.e. 660 km have been subducted, and lasts
throughout the nal stages of the subduction, until the unsubducted plate has reduced to  800 km. Slab and mantle velocities
and slab morphology are constant throughout the steady state.
The distribution of plate buoyancy along the trench direction
is the free parameter in our models. We have modeled it testing
three different buoyancy perturbations, that is varying the thickness throughout the slab length l (Suppl. Fig. 1, Suppl. Table 1,

F.A. Capitanio, M. Faccenda / Earth and Planetary Science Letters 353354 (2012) 2937

Models 111), varying the density contrast (Mods. 1219) and

varying the thickness of a smaller length m (Mods. 2021).
The thickness is varied by changing the rheology of the plates
lowermost 30 km to that of the mantle. The rheology is changed
in the outer portion of the slab width, and is viscous (i.e.,
Newtonian), 100 km thick in the center of the slab of width d
(Suppl. Fig. 1), modeled values of d are 100, 500, 200 km. In a
second setup we tested the role of a smaller length of plate of
varying thickness m to model self consistently the entrainment
of buoyancy heterogeneities. Finally, we tested the variation of
buoyancy keeping the thickness constant along the slab and
changing the density contrast (Suppl. Table 1). All models are
run with both a Newtonian and a non-Newtonian mantle.
After the perturbation the models run until the lithosphere is
entirely subducted. The measurement used for comparisons are
conducted 23 Myr after the slab is perturbed, this is the time
required to rearrange the ow.
Measured quantities are the pressure, where in the remainder
of this paper we refer to the dynamic (i.e. non-hydrostatic)
pressure, velocity and shear stress. To characterize the ow we
measure the vorticity, i.e. the curl of the velocity:

orv

Finally, we measure the dissipation partitioning and provide

an estimate of the dissipation in the lithosphere over the total:
R

FL RL

tij e_ ij dV
,
tij e_ ij dV

where L is the lithosphere volume. The lithospheric dissipation is

integrated in the volume of the plate where the viscosity has not

31

been decreased by plasticity, that is the shear zone at the trench is

not considered part of the subducting plate.

3. Results
3.1. The effect of buoyancy variations
We rst compare the model where the buoyancy is constant
along the slab with the varying buoyancy model, with
d2000 km (Fig. 1, Suppl. Table 1, Mods. 8 and 11). Large
buoyancy in the subducting plate center drives faster subduction
locally (Fig. 1A and D). The slab accommodates the differential
motions at depth through internal deformation. Yet, the differences in trench motions are negligible in these models, and
trenches remain straight due to the high viscosity of the upper
plate. The trench-ward velocity component, vx, measured at the
models surface, z 0, and y1000 km, differs due to the average
buoyancy increase, 4 and 6 cm yr  1 in constant and perturbed
buoyancy models, respectively. However the ow pattern is in
general very similar (Fig. 1B and E), forming a convective cell
sustained by the downgoing slabs and a vigorous return ow. The
component of the velocity vy that is perpendicular to the crosssection plane xz is negligible when the slab model has a constant
buoyancy (Fig. 1C), instead it is larger in the models where slab
buoyancy varies laterally (Fig. 1F). In this model a broad area
beneath the slab extending for  400 km from the slab where the
velocity vy is between 1 and 3 cm yr  1.
The viscosity prole beneath the subducting plate (Suppl.
Fig. S2A) is very similar for both models, reaching  7  1019 Pa s
in the asthenosphere, beneath the plate and slightly increasing to

Fig. 1. Homogeneous vs. heterogeneous models with non-Newtonian mantle. Constant thickness Mod. 8 (A, B, C) and variable thickness Mod. 11 (D, E, F) models at the
steady-state, contour of the density eld with velocity magnitude in color. Black line on top, trace of the upper plate. Thick purple line in D the location of the thickness
boundary. In yellow (A, D) the trace of the cross sections. Sections are taken at y 980 km. (B, E) Horizontal velocity component magnitude, vx, and velocity vectors. (C, F)
Velocity component perpendicular to the section plane, vy. Density contour in black line. Star symbol in E, F location of symbol in Fig. 2. (For interpretation of the
references to color in this gure legend, the reader is referred to the web version of this article.)

32

F.A. Capitanio, M. Faccenda / Earth and Planetary Science Letters 353354 (2012) 2937

mid-mantle depth, and decreasing at the bottom of the upper

mantle at  660 km depth. The component vy in the cross-section
plane is negligible when slab buoyancy is constant and it is as
high as 23 cm yr  1 when slab buoyancy varies along this
direction. In a vertical prole of the normalized velocity (Fig. 2)
we show the component vy is around a third of the vx component
in heterogeneous buoyancy models throughout the upper mantle
thickness, and is negligible in the constant buoyancy models.
Largest vy are found between 300 and 550 km depth, although the
whole upper mantle ow is characterized by a lateral velocity
component.
A depth at which the component vx is zero is found at
200 km. In this 30 km thick layer, although the ow velocity
is low, the trench-parallel ow is dominant (Fig. 2). The location

Fig. 2. Normalized velocity vertical proles from the two models in Fig. 1.
Location of the section is indicated in Fig. 1B by the star. Dashed lines, model in
Fig. 1A (Constant thickness slab), solid lines model in Fig. 1D (Variable
thickness slab).

of this point (Fig. 1E and F, star symbol) coincides with the

location where ow sheared by the plate reverts to the return
ow found deeper, and velocity in this plane is zero, i.e. no
horizontal or vertical motions. Yet, the lateral component is nonzero, so that the lateral transport in this layer is maximized. In the
plane at this depth pressures of  2 MPa occur beneath the
homogeneous slab, decreasing only at 1000 km from the slab
edge (Fig. 3A). Above the slab, in the mantle wedge, a broad low
pressure extends all along the slab, where pressure is  9 MPa.
The ow beneath the slab follows the pressure eld, thus
circulation is mostly perpendicular to the trench, except near
the slab edge, where toroidal ow develops. The along-trench
velocity, vy, is almost negligible around the slabs center, and is
only large close to slab edges, when it has constant thickness. The
path of particles immersed in the velocity eld in this plane
(Fig. 3B) shows that the mass exchange between the mantle
wedge and below the slab is limited to  1000 km from the slab
edge, whereas particles beneath the slab follow the trenchperpendicular path. Instead, in the model that includes plates
buoyancy heterogeneity, the ow is more complex and has a
stronger trench-parallel component (Fig. 3C). The thicker, more
negatively buoyant slab segment induces larger subduction velocity locally and a strong pressure gradient along the slab. When
slab thickness varies along the slab, the trench-parallel velocity vy
component is largest around the thickness variation (Fig. 3D). The
dynamic pressure gradient above the slab is larger, and results in
faster ow towards the center of the slab. The path of the particles
immersed in this ow shows that mass from beneath the center
of the slab can be transported along the trench-parallel direction
to the slab edge, and from there is entrained in the mantle wedge
and reach as far as the center of the 6000 km wide slab. In models
where the slab has constant thickness, the vorticity is large only
around the slab edge (Suppl. Fig. S3A). Instead, when the slab

Fig. 3. Detail of the horizontal (half) section through the models in Fig. 1. Depth of the sections is  200 km, indicated by star symbol in Fig. 2. Density contour in black line.
Dynamic pressure and velocity vectors of constant (homogeneous, A) and variable thickness (heterogeneous) models (C) with non-Newtonian mantle rheology. Velocity
arrows are not scaled. Lateral velocity vy and streak lines of particles indicated by small circle symbols.

F.A. Capitanio, M. Faccenda / Earth and Planetary Science Letters 353354 (2012) 2937

buoyancy varies laterally, the magnitude of the vorticity is larger

overall, since it is large around the internally deforming slab
(Suppl. Fig. S3B).
To better characterize the deformation in these models, we
have calculated the ratio:

 

@vy @vy
@vx @vx

8
r
@x
@z
@y
@z
indicating the ratio between the components of the velocity
gradient dening simple shear in the y direction over that in
the x direction. For the straight morphology of the slab model we
refer to the ratio of trench-parallel to trench-perpendicular
simple shear. Simple shear in constant slab buoyancy models is
negligible around subducting plates (Fig. 4A), it only increases
around the slab edges. This results in very low values of r. Instead,
in slabs of variable thickness values of r 1 are widespread in the
whole modeled upper mantle domain, indicating that (1) the
components of the simple shear are comparable, and (2) this is a
general characteristic of the ow. In this model, values of r are
larger beneath the slab between 100 and 350 km depth
(Fig. 4B), in an area 200 km wide, found beneath the trench in
our models. In this area trench-parallel simple shear is maximized, and in average r 42, reaching values of 89 closer to the
slab, so the trench-parallel simple shear is up to 89 times larger
than the trench-perpendicular component. Large values of r are
found beneath the upper plates in these models, however, this is
possibly the result of the very high viscosity in the topmost layers,
resulting in steep unrealistic velocity gradients.
3.2. Sensitivity of mantle ow to buoyancy perturbations
In Fig. 5 we summarize the results of the models 8 to 11, where
we have varied the thicker segment width, d. Measures are taken
from a plane beneath the slabs, at 450 km depth, where the ow
velocity is maximized. The dynamic pressures measured vary from
1.9 MPa in the center of constant thickness slab to 3.8 MPa beneath
the model with d2000 km. For d smaller than 2000 km, the
maximum pressure decreases, despite the buoyancy in this area
being the same. Gradients across the thickness variation are smooth
compared to the step-like buoyancy change. The subduction velocity

33

measured at depth, vsub (the magnitude of the ow velocity, outside

the slab) varies along the width of the slab, with the maximum
achieved in the center, where the buoyancy is increased. Here,
velocity ranges from 5 to  12.5 cm yr  1, from constant to varying
buoyancy slabs with d2000 km. Similar to pressure gradients,
velocity varies smoothly along the slab. Although the subduction
velocity in these models depends in the rst place on the slab
buoyancy (Capitanio et al., 2007), the velocity gradients along the
slab depend also on their stiffness, when their viscosity is nite. The
stiffness controls the propagation length of the buoyancy perturbation through the wavelength L 2phg=61=3 (Turcotte and Schubert,
1982), where g is the lithosphere/mantle viscosity contrast and h the
thickness of the lithosphere. In our models g 104 in the h 30 km
thick lithospheric core, yielding L 2200 km, and L/2 is the distance
from the step in thickness where differences between models become
small (Suppl. Fig. S4A). Although simplied, this explains why only
models where buoyancy anomaly is d L achieve a full subduction
velocity, whereas it is lower in models with d5L, despite their
buoyancy in the center being the same. Because the torque in the
mantle is applied on the plane containing the slab thickness h, slabs
accommodate sinking rate gradients, bending, readjusting dip or
migrating, resisted by a force ph. Instead, at plate margins the
torque, resulting from slab pull gradients applied on plane to the
plates width, has a larger resistance to deformation. This explains
why gradients in models surface velocity are small (Fig. 5B, thin
lines).
Trench-parallel velocity of the mantle beneath the slabs is
negligible when models buoyancy is constant, and increases with
increasing d, reaching the same velocity of the toroidal ow
around the edges in models with d 2000 km (Fig. 5C). We
explain it as the result of the subduction velocity gradients. The
pressure p beneath a slab of velocity vsub is (Stevenson and
Turner, 1977), that is the Couette ow:
p

2aZnsub
,
r

where a sin a/[(p  a)sin a] with a being the slab dip and Z is
the uid mantle viscosity and r is the downward distance along
the slab. Assuming that (1) the dip angle and, (2) the viscosity
variations along the slab are negligible, i.e. @a/@y 0 and @Z/@y0,

Fig. 4. Ratio r, Eq. (8), of trench-parallel to trench-perpendicular simple shear in constant (A) and variable thickness (B) models with non-Newtonian mantle (Fig. 1). The
cross section are taken at y 1080 km.

34

F.A. Capitanio, M. Faccenda / Earth and Planetary Science Letters 353354 (2012) 2937

understand the lateral ow component from the gradients of

subduction velocities at depth.
While this analysis holds for Newtonian mantle models as
well, where the characteristics of the ow are the same, introducing non-linearity in the mantle rheology affects the gradients of
strain patterns, i.e. strain rates. In the models, simple shear
accommodates the differential motions between the slab, where
vy  0, and the ow of largest vy, at  300 km from the slab in the
plane of the section. To illustrate this, we show the scaling of
simple shear strain rates e_ xy beneath the slab with the magnitude
of the velocity vy (Suppl. Fig. S4B).
The dissipation partitioning in these models is similar, increasing from FL 11 7 5%, in the constant buoyancy model, to
FL 17 7 4% in the model with d 2000 km. This shows that
the mantle drag force is very weakly dependent of the circulation
pattern.

4. Discussion
4.1. Comparison with previous studies and implications

Fig. 5. Comparisons of models with thickness variation along the slab width with
non-Newtonian mantle rheology. Constant thickness model mod. 8 black lines,
mod. 9 in red, mod. 10 in green and mod. 11 in light blue. Depth of section
indicated by star symbol in Fig. 2. (A) Pressure distribution in the subslab and
thickness variation along the slab width. (B) Subduction velocity in the slab (thick
lines) and velocity at the surface vsp (thin lines). (C) Lateral velocity component vy
(solid lines) increases with d and recovered velocity (dotted lines). (For interpretation of the references to color in this gure legend, the reader is referred to
the web version of this article.)

and (3) we use (9) which is for an innitely wide plate, the
pressure gradient in the direction y, perpendicular to the plane xz,
that is
@p
@p @vsub
2aZ @vsub

@y
@vsub @y
r @y

10

The pressure gradient @p/@y induces a ow of maximum

velocity vy, i.e. the Poiseuille ow, in an incompressible uid of
viscosity Z in a channel of height H this is (Turcotte and Schubert,
1982):

ny 

H2 @p
8Z @y

11

where we consider the velocity vy 0 at the slab surface. Substituting (10) in (11):
vy 

aH2 @vsub
@v
b sub ,
4r @y
@y

12

showing that mantle ow component vy develops as a result of

gradients in the subduction velocity. The lateral component of the
velocity is independent of the uid viscosity (Hall et al., 2000).
Taking r (450  h)/sin a 372 km, where a 701 is the average
dipping angle and H560 km is the subslab upper mantle, Eq.
(12) recovers comparable trench-parallel velocities beneath the
slab (Fig. 5C, dotted lines), introducing a tting constant c 2b,
despite the simplication assumed. This offers a rule of thumb to

Many studies have addressed the mantle ow associated with

subduction presenting two- and three-dimensional numerical
and laboratory models. While a review of these is not the aim
of this paper, we will focus here on the new perspective our
approach introduces.
Global models of subduction and mantle convection have
explained the pattern of plate motions, mantle ow and its
alignment with seismic anisotropy (Becker et al., 2003; Conrad
et al., 2007; Natarov and Conrad, 2012). These models derive the
large-scale buoyancy distribution of the mantle from tomographic
models, which account appropriately for the rst-order dynamics.
However, these models do not capture the smaller scale complexities of the ow around subduction zones.
Thus so far, only models where a realistic morphology of the
subduction zone has been considered can explain the complexities of the ow in the mantle wedge (Billen et al., 2003; Jadamec
and Billen, 2010; Kneller and van Keken, 2007) or as a global
feature (Natarov and Conrad, 2012). While they illustrate accurately the role of local complications in the mantle ow, making
inferences on a general ow mechanism from this approach
remains difcult.
The approach to buoyancy-driven subduction allows modeling
several of these relevant aspects at a glance, offering insight into
the fundamental mechanisms that relate motions and deformations in the mantle and plate. Variable buoyancy of plates result
in a range of subduction parameters including sinking velocity,
slab dip and radius, which are interdependent (Capitanio and
Morra, 2012; Capitanio et al., 2007). In a realistic three-dimensional setting, where these heterogeneities are present the resulting velocity and dip variations confer morphologies similar to the
observed as well as realistic continental tectonics (Capitanio et al.,
2011). Here we have shown that these heterogeneities also
account for complex mantle ow around slabs.
Other subduction zone regional features have the power to
modify the subduction parameters, such as the friction at the
plate margins (Iaffaldano et al., 2012) or the suction of upper
plates (Rodrgues-Gonzalez et al., 2012). However, these effects
act at shallow depths and cannot alter the overall sinking
dynamics (Capitanio and Morra, 2012; Capitanio et al., 2010b;
Dvorkin et al., 1993).
Along-trench buoyancy variations are required to explain
complex slab morphologies (Capitanio et al., 2011) and ow
around slabs, but also might potentially reduce the mist
between predicted and observed global plate motions. In fact,

F.A. Capitanio, M. Faccenda / Earth and Planetary Science Letters 353354 (2012) 2937

the overall mantle ow in our models is very similar to the 2D

ow usually reproduced by subduction models (e.g. Garfunkel
et al., 1986), where a convecting cell follows the downgoing
motions of the slabs. This is because the buoyancy perturbations
introduced are small enough that the fundamental dynamics is
not changed. Global models of plate motions in fact allow small
corrections to the global plate motions circuit, which have a
rather regional character (Kreemer, 2009; Stadler et al., 2010).
Although varying rheology might act to decouple plate motions
and slabs (Jadamec and Billen, 2010) to improve the t between
observed and predicted motions (Stadler et al., 2010), these still
have a regional extension. Here, we suggest that the t in plate
motions might be improved by considering the buoyancy structure of subducted slabs, rather than assuming constant thicknesses. This introduces realistic driving force gradients, and
possibly local variation to resisting forces could be less relevant.
4.2. Subduction motions, mantle ow and seismic anisotropy
Our modeling results indicate that the structure of the nonNewtonian mantle ow is very sensitive to the slab buoyancy
conguration, developing lateral ow rapidly. Yet, the relation
between mantle ow and surface plate motions is non-unique.
When stiff upper plates override subducting plates, the relative
motions of the slab are hampered, as seen here. Furthermore,
coupling gradients along the shear zone between plates might
also result in forced trench and slab motions at depth (Capitanio
et al., 2011; Iaffaldano et al., 2012). Here, we have shown that in
both cases the mantle ow at depth responds only to the buoyancy structure of the slab at depth. This can be understood by
considering that the buoyancy structure confers slab pull gradients, which can be accommodated within slabs by internal
deformations, thus stresses are not propagated to the trench. This
is the opposite case in models treating the slabs as innitely rigid
(Buttles and Olson, 1998; Hall et al., 2000), where the slab
motions at depth are strongly coupled with the plate and trench
motions at surface. Internal slab deformation as in the models
presented here is compatible with the complex morphology of the
slabs seen in the tomographic images (Grand et al., 1997), and
offers a viable process to decouple the surface plate motions with
the mantle ow, found by similar models (Becker et al., 2003;
Jadamec and Billen, 2010).
Complex mantle ow may have consequences for the interpretation of seismic anisotropy observed at subduction zones.
Here, the fast SKS component generally orients parallel to the
trench, which is explained by trench-parallel mantle circulation
around subducting plates (Long and Silver, 2008; Russo and
Silver, 1994). In our models with lateral slab density variations
the ow is effectively trench-parallel in just a 30 km thick layer
implying that, assuming a sub-parallel lattice preferred orientation (LPO) as it is done routinely by seismologists, it would not
produce substantial time delays (  0.3 s for 5% anisotropic peridotite). However, because anisotropic mantle minerals align
preferentially with the maximum stretching directions that coincides with the ow direction only during simple shear deformation (Kaminski and Ribe, 2002), the simple relation between LPO
orientation and ow direction is not always warranted, even in
the 30 km thick mantle layer with trench-parallel ow. Moreover,
in general it is difcult to predict the LPO development in the
subslab mantle where deformation is a mixture of different shear
and normal strain rate components. In order to make more robust
inferences on the resulting seismic anisotropy, an accurate
investigation that accounts for the whole deformational history
of the mantle is required (Faccenda and Capitanio, 2012).
Nevertheless, the fact that deformation occurs beneath modeled subducting slabs at a depth range of  100350 km by

35

simple shear rather than pure shear allows speculations on the

LPO development. Under conditions of trench-perpendicular vy
gradients (vy,x and vy,z) larger than trench-parallel vx gradients
(vx,y and vx,z), the LPO tends to be aligned with the trench. If the
opposite occurs the anisotropy should align perpendicular to the
trench. In our models with lateral slab buoyancy variations vy
gradients are in general comparable to/larger than vx, so that the
resulting strain ellipsoid should be strongly aligned sub-parallel
to the trench. Therefore, model results support a rather general
lateral ow component that might increase the trench-parallel
anisotropy component at the expense of the trench-perpendicular
component. Furthermore, thermal and mechanical (grains-size
dependent) rheological weakening feedbacks not considered here
act to decrease the coupling between the slab and subslab mantle
(Long and Silver, 2009), i.e. the asthenosphere, thus reducing
trench-perpendicular ow in the entrained mantle and associated
anisotropy.
Faccenda and Capitanio (2012) showed that trench-parallel
anisotropy develops beneath relatively narrow plates by pure
shear deformation and simultaneous toroidal ow when retreat
motions are maximized. In these cases, lateral pressure gradients
together with efcient subslab mantle decoupling from the overlying slab may then add to the background toroidal motions
yielding trench-parallel LPO orientation consistent with the
observed anisotropy patterns.
4.3. Implications for natural subduction zones
We have modeled the slabs negative buoyancy gradients,
mimicking similar thickness variations due to lithospheric age
variations at trenches. The validity of our outcomes can be
extended to other types of integrated buoyancy heterogeneities,
whether endogenous to the plates or exogenous.
In the subducting Nazca plate, present-day age at the trench
varies from 0, in the north and south plate divergent margins, to
 55 Ma in the center, thus varying from 0 to  90 km thickness.
The resulting Stokes velocity is thus expected to vary accordingly,
eventually driving slab deformation and upper plate tectonics
(Capitanio et al., 2011). All major subduction zones present

similar oceanic lithospheric age-gradients (Sdrolias and Muller,

2006).
Other subducting plate characteristics that have an impact on
the slab buoyancy are the nature of the lithopshere, whether
oceanic or continental, as well as the depth of subduction,
modulating the integrated slab pull. As an example, in the Tonga
subduction zone a mature oceanic plate is subducting along the
 2000 km trench. Although no thickness variation can be
invoked for this plate older than 80 Ma, towards the south
continental lithosphere is subducting in the Kermadec zone

(Sdrolias and Muller,

2006). The buoyancy gradients here might
be increased by the depth of subduction, reaching the lower
mantle beneath the Tonga but not deeper than 300 km beneath
Kermadec (Fischer et al., 1991; Lallemand et al., 2005). In this
subduction zone, convergence velocity decreases largely from
Tonga to Kermadec, from 12 to less than 1 cm yr  1. Although
speculative, large trench-parallel velocity gradients together with
rst order toroidal ow due to the maximized retreat motions
(Faccenda and Capitanio, 2012) could explain the trench-parallel
subslab seismic anisotropy measured (Long and Silver, 2008).
Typical features such as slab windows and detachments, but
also volcanic ridges, might confer slab pull gradients similar to
those idealized here, provided they are large enough to alter the
slab pull. We have shown that due to the stiffness of slabs, small
size buoyancy perturbations might not effectively alter the slab
velocity. This is compatible with the global plate motions compilations, where small ridges do not apparently alter the convergence

36

F.A. Capitanio, M. Faccenda / Earth and Planetary Science Letters 353354 (2012) 2937

motions along trenches (DeMets et al., 1990, 2010). Yet, the

simplied analysis in a previous section allows inferring that the
effect of smaller perturbations should be more evident on smaller
and/or younger plates. While in small plates the perturbations
have a larger contribution on the average buoyancy, young and
thin plates have smaller rigidity, so that perturbation should be
more effective in small and young over old and large oceanic
plates, resulting in larger deformation in young subducting slabs.

5. Conclusions
Buoyancy heterogeneities in slabs on Earth are primarily a
consequence of age-dependent thickness variations along
trenches. Three-dimensional numerical models of subduction
show that the mantle ow during subduction of heterogeneous
plates, along-trench variations of pull forces and sinking rates,
trigger mantle ow sub-parallel to the slab. For Earth-like slabs
stiffness the trench-parallel component becomes signicant only
when heterogeneities are wider than  1000 km, or less in young
subducting oceanic lithospheres. The lateral ow component does
not largely alter the convection pattern around the slab, however,
we found that (1) it introduces large trench-parallel simple shear
at depths of  100350 km and (2) at a depth of  200 km trenchparallel ow occurs. In the narrow layer of trench-parallel ow
material can be transported over large distance around very wide
slabs. Slabs internal deformation decouples the mantle ow from
motions at surface, such that on Earth trench and slab motions
might not always correspond. Therefore, we propose that heterogeneous subducting lithosphere buoyancy structures are required
to explain slabs convoluted morphologies, second-order plate
motions and complex mantle circulation. Thickness variations, as
well as slab break offs and windows, variable depth of subduction
and continental block entrainment are commonly found along
subduction zones. Because these all introduce similar buoyancy
heterogeneities along the trenches, complexities as those shown
here should have a global presence.

Acknowledgment
Discussions with T.W.Becker and L.N. Moresi are acknowledged. We thank two anonymous reviewers for insightful suggestions. This research was supported by the Australian Research
Council Discovery Projects DP0987374 and DP110101697 awarded
to F.A.C.

Appendix A. Supporting information

Supplementary data associated with this article can be found
in the online version at 10.1029/2001GC000222.

References
Afonso, J.C., Ranalli, G., Fernandez, M., 2007. Density structure and buoyancy of the
oceanic lithosphere revisited. Geophys. Res. Lett. 34, L10302, http://dx.doi.org/
10.1029/2007GL029515.

Becker, T.W., Kellogg, J.B., Ekstrom,

G., OConnell, R.J., 2003. Comparison of
azimuthal seismic anisotropy from surface waves and nite strain from global
mantle-circulation models. Geophys. J. Int. 155, 696714.
Billen, M., Gurnis, M., Simons, M., 2003. Multiscale dynamics of the Tonga
Kermadec subduction zone. Geophys. J. Int. 153, 359388.
Buttles, J., Olson, P., 1998. A laboratory model for subduction zone anisotropy.
Earth Planet. Sci. Lett. 164, 245262.
Capitanio, F.A., Faccenna, C., Zlotnik, S., Stegman, D.R., 2011. Subduction
dynamics and the origin of Andean orogeny and Bolivian Orocline. Nature
480, 8386, http://dx.doi.org/10.1038/nature10596.

Capitanio, F.A., Morra, G., 2012. The bending mechanics in a dynamic subduction
system: constraints from numerical modeling and global compilation analysis.
Tectonophysics 522523, 224234.
Capitanio, F.A., Morra, G., Goes, S., 2007. Dynamic models of downgoing platebuoyancy driven subduction: subduction motions and energy dissipation.
Earth Planet. Sci. Lett. 262, 284297.
Capitanio, F.A., Morra, G., Goes, S., Weinberg, R.F., Moresi, L., 2010a. IndiaAsia
convergence driven by the subduction of the Greater Indian continent. Nat.
Geosci. 3, 136139.
Capitanio, F.A., Stegman, D.R., Moresi, L., Sharples, W., 2010b. Upper plate controls
on deepsubduction, trenchmigrations and deformations at convergent
margins. Tectonophysics 483, 8092, http://dx.doi.org/10.1016/j.tecto.2009.
08.020.
Conrad, C.P., Behn, M.D., Silver, P.G., 2007. Global mantle ow and the development of seismic anisotropy: differences between the oceanic and continental
upper mantle. J. Geophys. Res. 112, http://dx.doi.org/10.1029/2006JB004608.
DeMets, C., Gordon, R.G., Argus, D.F., 2010. Geologically current plate motions.
Geophys. J. Int. 181, http://dx.doi.org/10.1111/j.1365-246X.2009.04491.x.
DeMets, C., Gordon, R.G., Argus, D.F., Stein, S., 1990. Current plate motions.
Geophys. J. Int. 101, 425478.
Dvorkin, J., Nur, A., Mavko, G., Ben-Avraham, Z., 1993. Narrow subducting slabs
and the origin of backarc basins. Tectonophysics 227, 6379.
Faccenda, M., Capitanio, F.A., 2012. Development of mantle seismic anisotropy
during subduction-induced 3-D ow. Geophys. Res. Lett. 39, L11305, http://dx.
doi.org/10.1029/2012GL051988.
Fischer, K.M., Creager, K.C., Jordan, T.H., 1991. Mapping the Tonga Slab. J. Geophys.
Res. 96, 1440314427.
Funiciello, F., Faccenna, C., Giardini, D., 2004. Role of lateral mantle ow in the
evolution of subduction systems: insights from laboratory experiments.
Geophys. J. Int. 157, 13931406.
Funiciello, F., Moroni, M., Piromallo, C., Faccenna, C., Cenedese, A., Bui, H.A., 2006.
Mapping mantle ow during retreating subduction: laboratory models analyzed by feature tracking. J. Geophys. Res. 111 , http://dx.doi.org/10.1029/
2005JB003792.
Garfunkel, Z., Anderson, C.A., Schubert, G., 1986. Mantle circulation and the lateral
migration of subducted slab. J. Geophys. Res. 91, 72057223.
Grand, S.P., van der Hilst, R.D., Widiyantoro, S., 1997. High resolution global
tomography: a snapshot of convection in the Earth. Geol. Soc. Am. Today 7,
17.
Hall, C.E., Fischer, K.M., Parmentier, E.M., Blackman, D.K., 2000. The inuence of
plate motions on three-dimensional back arc mantle ow and shear wave
splitting. J. Geophys. Res. 105, 2800928033.
Hoernle, K., Abt, D.L., Fischer, K.M., Nichols, H., Hauff, F., Abers, G.A., van den
Bogaard, P., Heydolph, K., Alvarado, G., Protti, M., Strauch, W., 2008. Arcparallel ow in the mantle wedge beneath Costa Rica and Nicaragua. Nature
451 , http://dx.doi.org/10.1038/nature06550.
Honda, S., Yoshida, T., 2005. Effects of oblique subduction on the 3-D pattern of
small-scale convection within the mantle wedge. Geophys. Res. Lett. 31,
L13307, http://dx.doi.org/10.1029/2005GL023106.
Iaffaldano, G., Di Giuseppe, E., Corbi, F., Funiciello, F., Faccenna, C., Bunge, H.P., 2012.
Varying mechanical coupling along the Andean margin: implications for
trench curvature, shortening and topography. Tectonophysics 526529,
1623.
Jadamec, M.A., Billen, M., 2010. Reconciling surface plate motions with rapid
three-dimensional mantle ow around a slab edge. Nature 465, 338341.
Kaminski, E., Ribe, N.M., 2002. Timescales for the evolution of seismic anisotropy
in mantle ow. Geochem. Geophys. Geosyst. 3, http://dx.doi.org/10.1029/
2001GC000222.
Kincaid, C., Grifths, R.W., 2003. Laboratory models of the thermal evolution of the
mantle during rollback subduction. Nature 425, 5862.
Kneller, E.A., van Keken, P.E., 2007. Trench-parallel ow and seismic anisotropy in
the Mariana and Andean subduction systems. Nature 450 , http://dx.doi.org/
10.1038/nature06429.
Kreemer, C., 2009. Absolute plate motions constrained by shear wave splitting
orientations with implications for hot spot motions and mantle ow. J.
Geophys. Res. 114, http://dx.doi.org/10.1029/2009JB006416.
Krien, Y., Fleitout, L., 2008. Gravity above subduction zones and forces controlling
plate motions. J. Geophys. Res. 113 , http://dx.doi.org/10.1029/2007JB005270.
Lallemand, S., Heuret, A., Boutelier, D., 2005. On the relationships between slab
dip, back-arc stress, upper plate absolute motion, and crustal nature in
subduction zones. Geochem. Geophys. Geosyst. 6, 917.
Long, M.D., Silver, P.G., 2008. The subduction zone ow eld from seismic
anisotropy: a global view. Science 319, 315319.
Long, M.D., Silver, P.G., 2009. Mantle ow insubduction systems: the subslab
ow eld and implications for mantle dynamics. J. Geophys. Res. 114,
315319, http://dx.doi.org/10.1029/2008JB006200.
Mehl, L., Hacker, B.R., Hirt, G., Kelemen, P.B., 2003. Arc-parallel ow within the
mantle wedge: evidence from the accreted Talkeetna arc, south central Alaska.
J. Geophys. Res. 108, http://dx.doi.org/10.1029/2002JB002233.
Miller, M.S., Kenneth, B.L.N., Toy, V.G., 2006. Spatial and temporal evolution of the
subducting Pacic plate structure along the western Pacic margin. J. Geophys. Res. 111, http://dx.doi.org/10.1029/2005JB003705.

Moresi, L., Dufour, F., Mulhaus,

H.B., 2003. A lagrangian integration point nite
element method for large deformation modeling of viscoelastic geomaterials.
J. Comput. Phys. 184, 476497.

F.A. Capitanio, M. Faccenda / Earth and Planetary Science Letters 353354 (2012) 2937

Morishige, M., Honda, S., 2011. Three-dimensional structure of P-wave anisotropy

in the presence of small-scale convection in the mantle wedge. Geochem.
Geophys. Geosyst. 12, Q12010, http://dx.doi.org/10.1029/2011GC003866.
Morra, G., Regenauer-Lieb, K., 2006. A coupled solid-uid method for modeling
subduction. Philos. Mag. 86, 33073323.
Morra, G., Regenauer-Lieb, K., Giardini, D., 2006. On the curvature of oceanic arcs.
Geology 34, 877880.

Muller,
R.D., Sdrolias, M., Gaina, C., Roest, W.R., 2008. Age, spreading rates, and
spreading asymmetry of the worlds ocean crust. Geochem. Geophys. Geosyst.
9, Q04006, http://dx.doi.org/10.1029/2007GC001743.
Natarov, S.I., Conrad, C.P., 2012. The role of Poiseuille ow in creating depthvariation of asthenospheric shear. Geophys. J. Int. 191, http://dx.doi.org/
10.1111/j.1365-246X.2012.05562.x.
Rodrgues-Gonzalez, J., Negredo, A.M., Billen, M., 2012. The role of the overriding
plate thermal state on slab dip variability and on the occurrence of at
subduction. Geochem. Geophys. Geosyst. 13, Q01002, http://dx.doi.org/
10.1029/2011GC003859.
Russo, R.M., Silver, P.G., 1994. Trench-parallel ow beneath the Nazca Plate from
seismic anisotropy. Science 263, 11051108.

37

Rychert, C.A., Fischer, K.M., Abers, G.A., Plank, T., Syracuse, E., Protti, J.M., Gonzalez,
V., Strauch, W., 2008. Strong along-arc variations in attenuation in the mantle
wedge beneath Costa Rica and Nicaragua. Geochem. Geophys. Geosyst. 9,
Q10S10, http://dx.doi.org/10.1029/2008GC002040.
Schellart, W.P., 2004. Kinematics of subduction and subduction-induced ow in
the upper mantle. J. Geophys. Res. 109, 2970.

Sdrolias, M., Muller,

R.D., 2006. Controls on back-arc basin formation. Geochem.
Geophys. Geosyst. 7 , http://dx.doi.org/10.1029/2005GC001090.
Sobolev, S.V., Babeyko, A.Y., 2005. What drives orogeny in the Andes? Geology 33,
617620.
Stadler, G., Gurnis, M., Burstedde, C., Wilcox, L.C., Alisic, L., Ghattas, O., 2010. The
dynamics of plate tectonics and mantle ow: from local to global scales.
Science 329, 10331038.
Stegman, D.R., Freeman, J., Schellart, W.P., Moresi, L., May, D., 2006. Inuence of
trench width on subduction hinge retreat rates in 3-D models of slab rollback.
Geochem. Geophys. Geosyst. 7, 122.
Stevenson, D.J., Turner, J.S., 1977. Angle of subduction. Nature 270, 334336.
Turcotte, D.L., Schubert, G., 1982. Geodynamics, Application of Continuum
Mechanics to Geological Problems. John Wiley & Sons, New York.