ARTICLE IN PRESS

Computers & Geosciences 30 (2004) 249–258

TEMSPOL: a MATLAB thermal model for deep subduction

zones including major phase transformations$
A.M. Negredoa, J.L. Valeraa, E. Carminatib,*
a
Department of Geophysics, Faculty of Physics, University Complutense of Madrid, Spain
b
Department of Earth Sciences, University of Roma ‘La Sapienza’, P. le Aldo Moro 5, Roma I-00185, Italy
Received 27 January 2003; received in revised form 2 January 2004; accepted 5 January 2004

Abstract

TEMSPOL is an open MATLAB code suitable for calculating temperature and lateral anomaly of density
distributions in deep subduction zones, taking into account the olivine to spinel phase transformation in a self-
consistent manner. The code solves, by means of a finite difference scheme, the heat transfer equation including
adiabatic heating, radioactive heat generation, latent heat associated with phase changes and frictional heating.
We show, with a few simulations, that TEMSPOL can be a useful tool for researchers studying seismic velocity, stress
and seismicity distribution in deep subduction zones. Deep earthquakes in subducting slabs are thought to be caused by
shear instabilities associated with the olivine to spinel phase transition in metastable olivine wedges. We investigate the
kinematic and thermal conditions of the subducting plate that lead to the formation of metastable olivine wedges.
Moreover, TEMSPOL calculates lateral anomalies of density within subducting slabs, which can be used to evaluate
buoyancy forces that determine the dynamics of subduction and the stress distribution within the slab.
We use TEMSPOL to evaluate the effects of heat sources such as shear heating and latent heat release, which are
neglected in commonly used thermal models of subduction. We show that neglecting these heat sources can lead to
significant overestimation of the depth reached by the metastable olivine wedge.
r 2004 Elsevier Ltd. All rights reserved.

Keywords: Temperature; Subduction; Phase transitions; Density anomaly; Olivine; Spinel

1. Introduction lithostatic pressures. Nevertheless, more than 20% of

earthquakes with magnitude greater than 5 occur at
Subcrustal earthquakes widely occur within subduct- depths greater than 300 km (Frohlich, 1989). Indepen-
ing slabs. Shallow and intermediate events (depths dent solutions to the paradox of the occurrence of deep
o300 km) are likely produced by brittle fracturing and earthquakes have been proposed by two groups of
frictional sliding on pre-existing faults (Scholz, 1990) scientists on the basis of deformation experiments on
and possibly enhanced by dehydration embrittlement germanate olivine and ice (Green et al., 1990; Kirby
(Meade and Jeanloz, 1991). At depths > 300 km; brittle et al., 1991). According to both groups, deep earth-
and frictional processes are likely inhibited by enormous quakes are caused by shear instabilities associated with
olivine to spinel transformation. This phase transforma-
$ tion may be kinetically hindered in a cold slab
Code available from server at http://www.iamg.org/
CGEditor/index.htm
descending into the mantle transition zone (e.g., Rubie,
*Corresponding author. Tel.: +39-6-49914950; fax: +39-6- 1984), thus creating a wedge-shaped zone of metastable
4454729. olivine that can persist down to about 660 km (Kirby
E-mail address: eugenio.carminati@uniroma1.it et al., 1996). This might explain the occurrence of deep
(E. Carminati). earthquakes and the largely recognised correlation of

0098-3004/$ - see front matter r 2004 Elsevier Ltd. All rights reserved.
doi:10.1016/j.cageo.2004.01.002
 ARTICLE IN PRESS
250 A.M. Negredo et al. / Computers & Geosciences 30 (2004) 249–258

maximum depth of seismicity with the product of the include shear or radiogenic heating or latent heat
age of the subducting lithosphere and the subduction released during phase changes. This model is commonly
velocity (e.g., Kirby et al., 1996). used to investigate the conditions for the formation of
Since olivine metastability is a temperature-controlled metastable olivine wedges (Kirby et al., 1996) and the
process, accurate and realistic thermal models of subduc- effects of buoyancy forces arising from density contrasts
tion zones are necessary for a better comprehension of related to phase transformations (e.g., Bina, 1997;
deep seismicity. Moreover, temperature and phase changes Marton et al., 1999). These models adopt a kinematic
are known to control density distribution within subduct- approach, i.e. a pre-defined velocity field is imposed. In
ing slabs and consequently to influence the dynamics of contrast, dynamic models of subduction calculate the
subduction (Bina, 1996; Schmeling et al., 1999). velocity field solving the coupled equations of conserva-
The role of subduction kinematics on the slab tion of mass, momentum and energy (e.g., Schmeling
temperature distribution has been widely investigated et al., 1999; Tezlaff and Schmeling, 2000), and include
with analytical (e.g., McKenzie, 1969) and numerical buoyancy forces due to phase transformation. Dynamic
(e.g., Minear and Toksoz, . 1970) thermal models since models are physically more self-consistent and provide a
the early discovery of plate tectonics. We aim at better understanding of subduction processes. However,
improving important aspects of existing models and they are more complex and difficult to apply to real
provide a useful tool for researchers studying deep subduction zones. Therefore, we preferred to follow a
subduction zones, with interests in seismology, meta- kinematic approach in order to create a simple code
morphism, seismic velocity distribution and dynamics of adequate to study subduction zones where basic input
subduction zones. The purpose of this work is, therefore, parameters (age of subducting lithosphere, dip and
to provide an open MATLAB code suitable for velocity of subduction) are known.
calculating temperature, density anomaly and mineral Our thermo-kinematic model improves some aspects
phases (olivine and spinel) distribution in deep subduc- of the above-mentioned works since it includes the main
tion zones, taking into account phase transformations in heat sources recognised in subduction zones: adiabatic
a self-consistent manner. Our code is intended to be heating, radioactive heat generation, phase changes and
simple and easy to apply to real subduction zones. frictional heating. Our model uses, as initial temperature
Moreover, we intend to evaluate the effects of heat for the lithosphere entering the subduction zone, the
sources such as radioactive or shear heating and latent thermal plate model GDH1 (Stein and Stein, 1992).
heat release, which are neglected in commonly used Moreover, the model calculates in a self-consistent way
thermal models (e.g., Stein and Stein, 1996). the location of the olivine to spinel phase transition, thus
permitting the calculation of an accurate distribution of
lateral anomalies of density.
2. Previous models

Among existing thermal models of deep subduction 3. The model

we can distinguish two different types of approach.
Some models compute the temperature distribution only A thermo-kinematic model is calculated by solving the
within the slab, which is assumed to descend in an 2D heat equation given by
isothermal (e.g., McKenzie, 1969) or in a adiabatically


LT qb qT qT qT
heated (e.g., D.assler and Yuen, 1996; Devaux et al., 1þ þ vx þ vz
1997, 2000) mantle. The influence of the surrounding cP qT qt qx qz

2 2 
mantle is accounted for as a boundary condition. K q T q T
¼ þ 2
A second set of models (e.g., Minear and Toksoz, . rcp qx2 qz
1970; Toksoz. et al., 1971; Turcotte and Schubert, 1973;

ag LT qb H þ Ash
Stein and Stein, 1996) includes mantle–slab thermal  vz Tabs þ þ ; ð1Þ
cp cp qz rcp
interactions and allows for the calculation of accurate
temperature distributions in both slab and surrounding where LT is the latent heat due to the olivine–spinel
mantle. Minear and Toksoz . (1970) and Toksoz . et al. transformation, cP is the heat capacity, b is the fraction
(1971) developed a numerical quasi-dynamic thermal of spinel, T is the temperature, t is time, x and z are the
model that included adiabatic heating, radioactive heat coordinates, vx and vz are the horizontal and vertical
generation, phase changes and frictional heating. Stein components of the velocity, K is the thermal conductiv-
and Stein (1996) further modified this model to ity, r is the density, a is the coefficient of thermal
introduce a temperature distribution of the oceanic expansion, g is the acceleration of gravity, Tabs is the
lithosphere prior to subduction given by the plate model absolute temperature, H is the radiogenic heat produc-
GDH1 (Stein and Stein, 1992). Despite this improve- tion rate and Ash is the shear heating rate. The values of
ment, the model by Stein and Stein (1996) does not the parameters used are listed in Table 1.
 ARTICLE IN PRESS
A.M. Negredo et al. / Computers & Geosciences 30 (2004) 249–258 251

Table 1 coordinate is located. The initial thermal distribution

Definition of parameters common to all calculations of the oceanic lithosphere is calculated with the thermal
Symbol Meaning Value plate model GDH1 of Stein and Stein (1992):


K Thermal conductivity 3:2 W m1 K1 0 za þ L  z X N
2 npðza þ L  zÞ
Tlit ðz; t Þ ¼ TL þ sin
cp Specific heat 1:3  103 J K1 kg1 L n¼1
np L
a Thermal expansion 3:7  105 K1
2 2 0 
n p Kt
coefficient  exp  ; ð2Þ
H Radiogenic heat 8  107 W m3 rcp L2
production rate (for where t0 is the time elapsed since the formation of the
cont. crust)
lithosphere, za and TL indicate location and temperature
r0 Density of the 3400 kg m3
of the base of the lithospheric plate, and L is the
lithospheric mantle at
T ¼ 0 C lithospheric thickness. Our code also includes the
Drol2sp Density change due to 181 kg m3 possibility of simulating continental subduction, with
ol–sp phase transition an initial geotherm for the crust and lithospheric mantle
L Lithospheric thickness 95; 000 m given by the steady-state solution of the heat conduction
T0 Temperature at the base 1450 C equation
of the lithosphere

H TL H 2 2H
g Gravitational 9:8 m s2 Tcrust ðzÞ ¼  ðza þ L  zÞ2 þ  hc þ hc
acceleration K L KL K
0 Subduction dip angle 60  ðza þ L  zÞ; ð3Þ
where hc is the crustal thickness. This relation assumes
T=0 oC x1 constant heat production limited to the crust and a
temperature T ¼ 0 C at the surface. The initial tem-
Lithosphere Dip (θ) A perature for the lithospheric mantle in the continental
za
A' case is given by
x0

TL H 2 H
Tlit mantle ðzÞ ¼  hc ðza þ L  zÞ þ h2c : ð4Þ
re

L KL K
he

Asthenosphere
sp
ho

We assume the following adiabatic initial temperature

Lit

Qx=0 Qx=0
profile for the asthenosphere:
B ztr2

ga
TðzÞ ¼ TL exp ðza  zÞ ; ztr2 ozpza : ð5Þ
B' cp
Olivine to Spinel transition
ztr1
We also include a total temperature increase of DT;
z caused by the latent heat released during the olivine to
Transition zone
spinel phase change, which initiates at ztr2 and is
x Qz =Qb
completed at ztr1 (see Fig. 1). In this interval the
temperature increase is given by (Turcotte and Schubert,
Fig. 1. Schematic diagram illustrating model domain and 2002)
applied boundary conditions. Velocity field imposed to simulate
LT
subduction is updated as slab penetrates into surrounding DT ¼ : ð6Þ
mantle. Small points indicate nodes of finite difference grid. cp
We assume that this increase occurs linearly from a
temperature Tztr2 at ztr2 to Tztr1 at zir1 :
The terms on the left-hand side represent the rate of
DT
heat change due to the temperature change at a fixed TðzÞ ¼ ðztr2  zÞ þ Tztr2 ; ztr1 pzpztr2 : ð7Þ
ztr2  ztr1
point and advection. The first term on the right-hand
side represents heat conduction, the term containing Tabs In the underlying transition zone the initial temperature
describes the adiabatic heating, that containing LT is given by


represents the latent heat of the olivine to spinel ga
transformation, and the last term accounts for radio- TðzÞ ¼ ðTztr2 þ DTÞ exp ðztr1  zÞ ;
cp
genic and dissipative heating. For simplicity, ‘spinel’ is
0pzpztr1 : ð8Þ
here used indistinctly for wadsleyite and/or ringwoodite
because of their similar densities. In order to find the depth at which phase transforma-
The model extends vertically to the base of the upper tions begin and end, we have followed the approach
mantle (Fig. 1) where the origin of the vertical described by Schmeling et al. (1999) and Tezlaff and
 ARTICLE IN PRESS
252 A.M. Negredo et al. / Computers & Geosciences 30 (2004) 249–258

Temperature (oC) the base of the model. The latter is calculated as the
0 400 800 1200 1600 product between the vertical gradient of the initial
200 geotherm at the base of the model and the thermal
conductivity.

mantle
Olivine
cold In our thermo-kinematic approach the velocity field
1% (Fig. 1) is imposed and defined by the subduction
Depth (km)

400 velocity vs and the slab dip y: In the corner area (shaded
slab

war
in Fig. 1), where the plate begins to dip into the mantle,
Metastable mer Spinel the velocity field is given by
Olivine slab z  za
vx ¼ vs qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi;
600 99% ðx  x0 Þ2 þ ðz  za Þ2
x  x0
vz ¼ vs qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi; ð10Þ
ðx  x0 Þ2 þ ðz  za Þ2
Fig. 2. Simplified disequilibrium (kinetic) phase diagram of
olivine–spinel transition showing lines of 1% and 99% spinel where x0 is the horizontal coordinate of the subduction
fraction as simplified by Schmeling et al. (1999). Thick lines knee (see Fig. 1). In the subducting slab it is given by
represent characteristic temperature paths within a cold slab, a
vx ¼ vs cos y; vz ¼ vs sin y: ð11Þ
warmer slab and undisturbed mantle.
The modelled region is divided into a 2D grid with
horizontal and vertical spacing Dx and Dz; with Dx ¼
Dz cos y in order to force grid nodes to lie on the edges
Schmeling (2000). Instead of solving the kinetic equa-
of the slab. At each time step p, the position (B–B0 in
tions of phase transformation we adopt a simplified
Fig. 1) of the bottom edge of the slab (which was located
phase diagram (Fig. 2). We have incorporated in our
along A–A0 at t ¼ 0) is computed, and the slab velocity
modelling the simplified phase diagram of Schmeling
field (Eq. (11)) is applied to the nodes reached by
et al. (1999) from original data by Akaogi et al. (1989)
line B–B0 :
and Rubie and Ross (1994). The 1% and 99% lines
We use the following expression (Turcotte and
separate the olivine and spinel regimes. The vertical lines
Schubert, 2002) for calculating the shear heating rate:
represent the transition from metastable olivine into
vs
spinel and are only crossed by the geotherms of cold Ash ¼ t ð12Þ
slabs (Fig. 2). The shallower limit of the interval of w
phase transformation within the undisturbed mantle where t is the shear stress and w is the thickness of the
ðztr2 Þ is given by the intercept between the asthenosphere shear zone. Optionally, the code incorporates the
initial geotherm (Eq. (4)) and the shallower straight line procedure proposed by Ponko and Peacock (1995) of
of the phase diagram. ztr1 is then obtained by substitut- distributing shear heating within a 2Dx wide zone
ing Tztr2 þ DT into the equation of the phase transition centred on the interface between the slab and the
deeper straight line. This diagram is also used to overlying mantle.
compute the fraction of spinel ðbÞ and its derivatives at The finite difference formulation of the energy
any temperature and depth. For the sake of simplicity, b equation (Eq. (1)) is developed in Appendix A. We have
is assumed to change linearly from b ¼ 0 in the olivine applied the alternating direction implicit (ADI) method.
and metastable olivine field to b ¼ 1 in the spinel field. The resulting equations are solved by applying the
We assume the following relation for lateral anomaly Thomas algorithm, as shown in Appendix B. Since the
of density (with respect to an unperturbed column) due velocity field is updated as the slab penetrates into the
to thermal expansion and phase changes: mantle, the numerical scheme applied is not stable for
every combination of Dt; Dx and Dz; the time step is
Dr ¼ ðr0 þ Drol2sp DbÞð1  aðT  T0 ÞÞ; ð9Þ
internally computed by the code in order to ensure
where r0 is the density at zero temperature, Drol2sp stability.
is the density increase caused by the olivine to spinel
phase transition, Db is the contrast of fraction of
spinel with respect to an unperturbed column with a 4. Some simulations
geotherm T0 : No density changes induced by high
pressure–low temperature metamorphism are taken into We have performed several calculations to show the
account. code behaviour for a wide variety of subduction condi-
The applied thermal boundary conditions (Fig. 1) are: tions. Fig. 3 shows the calculated temperatures and the
T ¼ 0 C at the surface, null horizontal heat flux at location of the olivine–spinel (1% and 99%) transfor-
the lateral boundaries and constant vertical heat flow at mation obtained for the cases of slow subduction of a
 ARTICLE IN PRESS
A.M. Negredo et al. / Computers & Geosciences 30 (2004) 249–258 253

0 200 200
600 600
1000

20
1000

00
100 1400

10
600

0
80
1400 0
40
200

0
60

00
Depth (km)

14
00
300

12
00
10
0

0
80

140
400 1% Spinel

1600

0
00

100
16 99% Spinel

00
500

12
00
12

00
600 14

0 100 200 300 400 500 600

(a) Horizontal distance (km)

0 200 200
600 600
1000 1000

0
40
100 1400
0 1400
20

00
12
0
40

200
0
80
Depth (km)

0
60

300
0
40

400 1% Spinel
1600
1600
00 0

99% Spinel
00 100 6

500

600
14

0 100 200 300 400 500 600

(b) Horizontal distance (km)

Fig. 3. Temperature ( C; contours every 200 C) distributions calculated with (a) a model of slow subduction of a warm (30 Myr old)
oceanic lithosphere after 32 Myr and (b) fast subduction of a cold (120 Myr old) oceanic lithosphere after 6:4 Myr: Subduction
velocities are (a) 2 cm year1 and (b) 10 cm yr1 : Position of 1% and 99% olivine–spinel transformation contours are shown by thin
dashed lines.

warm (young) oceanic lithosphere and fast subduction heating or radiogenic production are considered in both
of a cold (old) oceanic lithosphere, after 32 and 6:4 Myr cases.
of subduction, respectively. Values assumed for the Both models show (Fig. 3) lower temperature in the
velocity of subduction are 2 and 10 cm yr1 and ages are slab interior with respect to the slab boundaries.
30 and 120 Myr; respectively. The values of the thermal Moreover, slabs with high thermal parameters show
parameter f; calculated as the product of the vertical lower temperatures at fixed depths. Fig. 3 clearly shows
component of the subduction rate and the age of the the influence of subduction kinematics and pre-subduc-
subducting lithosphere, are 520 and 10; 392 km; respec- tion thermal state of the lithosphere on the position of
tively. A latent heat of phase transformation of LT ¼ the olivine to spinel phase change. In slow and warm
90 kJ kg1 (Turcotte and Schubert, 2002), and no shear subduction zones (Fig. 3a), the phase change is
 ARTICLE IN PRESS
254 A.M. Negredo et al. / Computers & Geosciences 30 (2004) 249–258

40
100

40

0
12
200

40
Depth (km)

300

80
120
400 240 1% Spinel

0
16
200
500 40 99% Spinel

600 80

0 100 200 300 400 500 600

(a) Horizontal distance (km)

100
0
0 6
80 12 1
40

200
Depth (km)

40

300
0 0
16 24

400 1% Spinel
0
12
40

0
500 12 99% Spinel
80

600
40

0 100 200 300 400 500 600

(b) Horizontal distance (km)

Fig. 4. Lateral anomaly of density (kg m3 contours every 40 km m3 ) with respect to right boundary of model domain calculated
with (a) a model of slow subduction of a warm (30 Myr old) oceanic lithosphere after 32 Myr and (b) fast subduction of a cold
(120 Myr old) oceanic lithosphere after 6:4 Myr: Subduction velocities are (a) 2 cm yr1 and (b) 10 cm yr1 : Position of 1% and 99%
olivine–spinel transformation contours are shown by thin dashed lines.

shallower in the slab than in the surrounding mantle, isotherm, which controls the initiation of metastable
due to depressed isotherms in the slab and to the slope of olivine to spinel transformation, reaches maximum
the equilibrium curve (Fig. 2), and reaches shallowmost depths of about 545 km and the conditions for the
depths in the core of the slab. Colder and faster slabs formation of a wedge of metastable olivine are attained.
(Fig. 3b) also show a shallowing of the transformation Deep seismicity is therefore expected to occur for
depth within the slab, but the innermost portions of the subduction processes with these kinematic and initial
slab show a deepening of the transition due to the fact thermal conditions.
that the temperature path in this region enters the The distributions of lateral anomalies of density
metastable olivine field (Fig. 2). In this case, the 600 calculated for both models are shown in Fig. 4. Down
 ARTICLE IN PRESS
A.M. Negredo et al. / Computers & Geosciences 30 (2004) 249–258 255

0
300
50

50
350

20
20
20
100

-5
400

80
Depth (km)

50
80

Depth (km)
150 450 LT = 90 kJ kg -1

50
Ash = 9.5 µW m-3

500
200 LT = 90 kJ kg-1
Ash = 0 µW m-3
20 550
250 LT = 0 kJ kg-1
Ash = 0 µW m-3
600
300
650
100 150 200 250 300 350 400 450
300 350 400 450 500 550 600
Distance (km)
Horizontal distance (km)
Fig. 6. Influence of latent heat release and shear heating during
Fig. 5. Lateral anomaly of density (kg m3 contours every
fast oceanic subduction on depth reached by metastable olivine
30 km m3 ) with respect to right boundary of model domain,
after 4 Myr of continental subduction with a velocity of wedge. Position of 1% olivine–spinel transformation ðb ¼ 0:01Þ
subduction of 5 cm yr1 : is shown.

to depths of about 300 km; slow slabs show lower geometry of the metastable olivine wedge, for the model
densities contrasts than fast slabs, due to higher of fast subduction of a cold slab (same parameters as
temperatures. At depths greater than 300–350 km; the Figs. 3b and 4b). We have considered a shear stress t ¼
shallowing of the olivine to spinel transition within the 30 MPa and a thickness of the shear zone w ¼ 2Dx; to
slab produces a dramatic increase of density in the case obtain with Eq. (12) a shear heating rate of Ash ¼
of slow subduction of a young plate (Fig. 4a). A negative 9:5 mW m3 : We show that neglecting these heat sources
buoyancy force arises from this positive density con- can lead to significant overestimation (about 85 km in
trast. In contrast, it may be noted that the presence of a the case shown in Fig. 5) of the depth reached by the
metastable olivine wedge in cold and fast slabs induces a metastable olivine wedge, and thereby of the predicted
remarkable lowering of density in the innermost maximum depth of seismicity.
portions of the slabs with respect to the surrounding
spinel phase. The resulting positive buoyancy can
significantly counteract the negative buoyancy men- 5. Discussion
tioned before, thus producing the so-called ‘parachute
effect’ (Kirby et al., 1996; Schmeling et al., 1999; Tezlaff We have presented here a code suitable for calculating
and Schmeling, 2000). temperature, mineral phases and density anomaly
Fig. 5 shows the lateral anomalies of density resulting distribution in deep subduction zones taking into
from subduction of continental lithosphere, after 4 Myr account the olivine–spinel phase transformation in a
with a velocity of subduction of 5 cm yr1 : A constant self-consistent manner. With respect to commonly used
density of 2800 kg m3 has been assumed for a 30 km thermal models, which neglect shear or radiogenic
thick continental crust. For simplicity, we have not heating or latent heat released during phase changes,
included density changes due to metamorphism, but we show that neglecting these heat sources can lead to
common phase diagrams can be easily incorporated in significant overestimation of the predicted maximum
TEMSPOL. The high negative density anomaly in the depth of seismicity. Further suitable refinements of the
subducted crust causes a positive buoyancy force that present model are including more complicated phase
opposes subduction, which is then only expected to diagrams and a variable thermal conductivity (Hauck
continue if the continental slab is pulled by a deeper et al., 1999).
oceanic portion of the slab. Modelling outputs can be used as inputs for numerical
Heat sources such as shear heating and latent heat simulations which require the knowledge of temperature
release are often neglected in thermal models of and density anomaly distributions of subducting zones.
subduction zones. TEMSPOL permits to include such For instance, the computed temperature can be used to
heat sources and to evaluate their contribution to the calculate the brittle and ductile yield stress distribution,
thermal state of subduction zones. Fig. 6 shows the assuming a nonlinear rheology, and correlate the
influence of latent heat and shear heating on the obtained stress distribution with observed seismicity
 ARTICLE IN PRESS
256 A.M. Negredo et al. / Computers & Geosciences 30 (2004) 249–258

patterns in oceanic and continental slabs (e.g., Carmi- equation is

nati et al., 2002). Another possible application relevant pþ1 pþ1
Ti;jpþ1  Ti;jp Tiþ1;j  Ti1;j
for tomography studies is to use the temperature Rpi;j þ vpxi;j
anomalies and temperature derivatives of seismic Dt 2Dx
!
velocities to obtain the seismic velocity perturbation p p
Ti;jþ1  Ti;j1
distribution (e.g., Sleep, 1973; De Jonge et al., 1994; þ vpzi;j
2Dz
Deal et al., 1999) induced by subduction processes.
pþ1
These synthetic tomographic images can be compared K Tiþ1;j  2Ti;jpþ1 þ Ti1;j
pþ1

with tomographic models of subduction zones. ¼

rpi;j cp Dx2
The resulting lateral anomaly of density distribution !
p
can be introduced into dynamic models assuming either Ti;jþ1  2Ti;jp p
þ Ti;j1
þ
an elastic or viscoelastic rheology (e.g., Bina, 1997; Dz2
Wittaker et al., 1992; Giunchi et al., 1996; Negredo et al.,
 p
!
1999) to evaluate buoyancy forces and the resulting LT qb p agTabs i;j
 vpzi;j þ
stress regime and compare it with seismotectonic cp qz i;j cp
observations in different subduction zones. Other Hi;j þ ti;j ðvs =wÞ
dynamic effects of the presence of metastable olivine þ ; ðA:1Þ
rpi;j cp
that can be investigated with dynamic models are the
change in the velocity of subduction and in the where
topography of the trench.

LT qb p
Rpi;j ¼ 1 þ :
cp qT i;j

Acknowledgements After grouping terms Eq. (A.1) becomes

pþ1
Ai Ti1;j þ Bi Ti;jpþ1 þ Ci Tiþ1;j
pþ1
¼ Di ; ðA:2Þ
The authors are very grateful to P. Vacher and an
anonymous reviewer for their constructive revision of
where
the manuscript. C. Doglioni is thanked for encourage-
ment and for fruitful discussions. MIUR funding Rpi;j p K
Ai ¼  v  ;
(C. Doglioni) is acknowledged. Spanish research pro- 2Dx xi;j rpi;j cp Dx2
! y Cajal, REN2001-3868-C03-02/MAR, and
jects Ramon
BTE2002-02462 have financially supported this study.
Rpi;j 2K
Italian MURST ‘Progetto Giovani Ricercatori’ Bi ¼ þ ;
(E. Carminati) supported this study. Dt rri;j cp Dx2

Rpi;j p K
Ci ¼ v 
Appendix A. Finite difference formulation of the energy 2Dz zi;j rpi;j cp Dx2
equation
and
The two-dimensional ADI method is a scheme in Rpi;j p Rpi;j p
which the dependent variable is solved for implicitly in Di ¼ Ti;j  v ðT p  Ti;j1
p
Þ
Dt 2Dz zi;j i;jþ1
the one direction and then in the other direction at time K
steps p þ 1 and p þ 2; respectively. þ p ðT p  2Ti;jp þ Ti;j1
p
Þ
ri;j cp Dz2 i;jþ1
We use the following convention for the grid system:
 !
p
LT qb p agTabsi;j Hi;jp þ tpi;j ðvs =wÞ
 vpzi;j þ þ :
x ¼ iDx; i ¼ 1; 2; y; I; cp qz i;j cp rpi;j cp
z ¼ jDz; j ¼ 1; 2; y; J;
An equation similar to Eq. (A.2) is obtained at each
t ¼ pDt; p ¼ 1; 2; y; P; point in the jth row, yielding a set of I equations. This
Ti;jp ¼ TiDx;jDz
pDt
: set of equations gives the temperature at p þ 1 at each
point of the jth row in terms of variables known at the
time step p: We apply the Thomas’ algebraic scheme to
Consider the finite difference expression of Eq. (1), in solve this set of equations (see Appendix B).
which temperature derivatives with respect to x are We now consider the finite difference formulation of
evaluated at time step p þ 1; and derivatives with respect Eq. (1) in which temperature derivatives with respect to
to z at time step p: At grid point i; j the resulting x are evaluated at time p þ 1; and derivatives with
 ARTICLE IN PRESS
A.M. Negredo et al. / Computers & Geosciences 30 (2004) 249–258 257

respect to z at time p þ 2: Appendix B. Application of Thomas’ algorithm to the

pþ1 pþ1
ADI scheme of the energy equation
Ti;jpþ2  Ti;jpþ1 Tiþ1;j  Ti1;j
Rpi;j þ vpxi;j
Dt 2Dx We consider the system of equations obtained at time
pþ2 pþ2
! step p þ 1:
Ti;jþ1  Ti;j1
þ vpzi;j Ai Ti1;j þ Bi Ti;j þ Ci Tiþ1;j ¼ Di ;
2Dz
pþ1 pþ1 pþ1 i ¼ 1; 2; y; I  1; ðB:1Þ
K Tiþ1;j  2Ti;j þ Ti1;j
¼ p where we are omitting the superscripts. The Thomas
ri;j cp Dx2
! algorithm operates by reducing this system of equations
pþ2 pþ2 pþ2
Ti;jþ1  2Ti;j þ Ti;j1 to upper triangular form, by eliminating Ti1 in each of
þ
Dz2 the equations. Consider that the first k equations (B.1)

 pþ1
! have been reduced to
p LT qb pþ1 agTabsi;j
 vzi;j þ Ti;j  wi Tiþ1;j ¼ gi ; i ¼ 1; 2; y; k; ðB:2Þ
cp qz i;j cp
Hi;jp þ tpi;j ðvs =wÞ the last of these equations is
þ ðA:3Þ
rpi;j cp Tk;j  wk Tkþ1;j ¼ gk ðB:3Þ

which, grouping terms, gives and the next equation in the original form is

pþ2 Akþ1 Tk;j þ Bkþ1 Tkþ1;j þ Ckþ1 Tkþ1;j ¼ Dkþ1 : ðB:4Þ

Aj Ti;j1 þ Bj Ti;jpþ2 þ Cj Ti;jþ1
pþ2
¼ Dj ; ðA:4Þ
We eliminate Tk;j from (B.3) and (B.4) and rearrange
where terms
Rpi;j p K Ckþ1 Dkþ1  Akþ1 gk
Aj ¼  v  ; Tkþ1;j  Tkþ2;j ¼ ðB:5Þ
2Dz xi;j rpi;j cp Dz2 Bkþ1 þ Akþ1 wk Bkþ1 þ Akþ1 wk
comparison with Eq. (B.2) shows that the coefficients wi
and gi can be obtained from the following recurrence
Rpi;j 2K
Bj ¼ þ ; relations:
Dt rpi;j cp Dz2
Ci Di  Ai gi1
wi ¼ gi ¼
Bi þ Ai wi1 Bi þ Ai wi1
Rpi;j p K
Cj ¼ v  i ¼ 2; 3; y; I: ðB:6Þ
2Dz zi;j rpi;j cp Dz2
The initial values can be obtained by applying in
and Eq. (B.2) the condition of zero horizontal heat flow at
the lateral boundaries, T1;j ¼ T2;j ; then
Rpi;j pþ1 Rpi;j p
Dj ¼ T  v ðT pþ1  Ti1;jpþ1
Þ w1 ¼ 1; g1 ¼ 0: ðB:7Þ
Dt i;j 2Dx xi;j iþ1;j
K We use these initial values to recursively find all the
þ p ðT pþ1  2Ti;jpþ1 þ Ti1;j
pþ1
Þ
ri;j cp Dx2 iþ1;j coefficients. Then we can calculate the last value of T by

 pþ1
! applying again the boundary condition TI;j ¼ TI1;j in
LT qb pþ1 agTabs i;j
 vpzi;j þ Eq. (B.2):
cp qz i;j cp
gI1
Hi;jp þ tpi;j ðvs =wÞ TI;j ¼ : ðB:8Þ
1  wI1
þ :
rpi;j cp
Therefore, beginning from the known value of TI;j ;
Eq. (B.2) gives the values TI1;j ; TI2;j ; in order,
An equation of form (A.4) can be obtained for each
finishing with T1;j : This computation is stable if jwi jo1:
point of the ith column, yielding J simultaneous
In the following time step, p þ 2; the temperature
equations. The I sets of J equations obtained at time
derivatives with respect to z are evaluated at p þ 2; and
p þ 2 are solved as before applying the Thomas
derivatives with respect to x at time p þ 1: We now apply
algorithm (Appendix B). Therefore, the temperature at
the boundary conditions at the base and surface of the
time p þ 2 is obtained from temperature at time p in two
model. We assume constant basal heat flow Qb to obtain
time steps. An equation implicit in x (Eq. (A.2)) is used
for the first time step and an equation implicit in z Qb
Ti;2 ¼ Ti;1  Dz: ðB:9Þ
(Eq. (A.4)) for the second. K
 ARTICLE IN PRESS
258 A.M. Negredo et al. / Computers & Geosciences 30 (2004) 249–258

We then apply Eq. (B.2) Kirby, S.H., Stein, S., Okal, E., Rubie, D.C., 1996. Metastable

 mantle phase transformations and deep earthquakes in
Qb
Ti;1  w1 Ti;1  Dz ¼ g1 ; ðB:10Þ subducting oceanic lithosphere. Reviews of Geophysics 34,
K 261–306.
so if w1 ¼ 1; then g1 ¼ ðQb =KÞDz: We then apply Marton, F.C., Bina, C.R., Stein, S., Rubie, D.C., 1999. Effects
Eq. (B.6) to find recursively the rest of the coefficients. of slab mineralogy on subduction rates. Geophysical
Finally, with the condition of constant surface tempera- Research Letters 26, 119–122.
McKenzie, D.P., 1969. Speculations on the consequences and
ture of 0 C; so Ti;J ¼ 0; we can apply Eq. (B.2) to
causes of plate motions. Geophysical Journal of the Royal
calculate Ti;J1 ; y; Ti;1 :
Astronomical Society 188, 1–32.
Meade, C., Jeanloz, R., 1991. Deep-focus earthquakes an
References recycling of water into the Earth’s mantle. Science 252,
68–72.
Akaogi, M., Ito, E., Navrotsky, L., 1989. Olivine-modified Minear, J.W., Toksoz, . M.N., 1970. Thermal regime of a
spinel–spinel transitions in the system Mg2 SiO4  Fe2 SiO4 : downgoing slab and new global tectonics. Journal of
calorimetric measurements, thermochemical calculation and Geophysical Research 75, 1397–1419.
geophysical application. Journal of Geophysical Research Negredo, A.M., Carminati, E., Barba, S., Sabadini, R., 1999.
94, 15671–15685. Dynamic modelling of stress accumulation in Central Italy.
Bina, C.R., 1996. Phase transition buoyancy contributions to Geophysical Research Letters 26, 1945–1948.
stresses in subducting lithosphere. Geophysical Research Ponko, S.F., Peacock, S.M., 1995. Thermal modelling of the
Letters 23, 3563–3566. southern Alaska subduction zone: insight into the petrology
Bina, C.R., 1997. Patterns of deep seismicity reflect buoyancy of the subducting slab and overlying mantle wedge. Journal
stresses due to phase transitions. Geophysical Research of Geophysical Research 100, 22117–22128.
Letters 24, 3301–3304. Rubie, D.C., 1984. The olivine–spinel transformation and
Carminati, E., Giardina, F., Doglioni, C., 2002. Rheological the rheology of subducting lithosphere. Nature 308,
control on subcrustal seismicity in the Apennines subduc- 505–508.
tion (Italy). Geophysical Research Letters 29, doi:10.1029/ Rubie, D.C., Ross, C.R., 1994. Kinetics of the olivine–spinel
2001GL014084. transformation in subducting lithosphere: experimental
D.assler, R., Yuen, D.A., 1996. The metastable olivine wedge in constraints and implications for deep slab process. Physics
fast subducting slabs; constraints from thermo-kinetic of Earth and Planetary Interiors 86, 223–241.
coupling. Earth and Planetary Science Letters 137, 109–118. Schmeling, H., Monz, R., Rubie, D.C., 1999. The influence of
Deal, M., Nolet, G., Van der Hilst, R.D., 1999. Slab olivine metastability on the dynamics of subduction. Earth
temperature and thickness from seismic tomography. and Planetary Science Letters 165, 55–66.
Journal of Geophysical Research 104, 28789–28802. Scholz, C.C.H., 1990. The Mechanics of Earthquakes and
De Jonge, M.R., Wortel, M.J.R., Spakman, W., 1994. Regional Faulting, Cambridge University Press Cambridge, UK,
scale tectonic evolution and the seismic velocity structure of 439pp.
the lithosphere and upper mantle: the Mediterranean Sleep, N.H., 1973. Teleseismic P-wave transmission through
region. Journal of Geophysical Research 99, 12091–12108. slabs. Bulletin of the Seismological Society of America 63,
Devaux, J.P., Schubert, G., Anderson, C., 1997. Formation of a 1349–1373.
metastable olivine wedge in a descending slab. Journal of Stein, C.A., Stein, S., 1992. A model for the global variation in
Geophysical Research 102, 24627–24637. oceanic depth and heat flow with lithospheric age. Nature
Devaux, J.P., Fleitout, L., Schubert, G., Anderson, C., 2000. 359, 123–129.
Stresses in a subducting slab in the presence of a meta- Stein, S., Stein, C.A., 1996. Thermo-mechanical evolution of
stable olivine wedge. Journal of Geophysical Research 105, oceanic lithosphere: implications for the subduction process
13365–13373. and deep earthquakes. Subduction from top to bottom.
Frohlich, C., 1989. The nature of deep focus earthquakes. Geophysical Monograph 95, 1–17.
Annual Review of Earth and Planetary Sciences 17, 227–254. Tezlaff, M., Schmeling, H., 2000. The influence of olivine
Giunchi, C., Sabadini, R., Boschi, E., Gasperini, P., 1996. metastability on deep subduction of oceanic lithosphere.
Dynamic models of subduction: geophysical and geological Physics of Earth and Planetary Interiors 120, 29–38.
evidence in the Tyrrhenian Sea. Geophysical Journal . M.N., Minear, J.W., Julian, B.R., 1971. Temperature
Toksoz,
International 126, 555–578. field and geophysical effects of a downgoing slab. Journal of
Green, H.W., Young, T.E., Walker, D., Scholz, C.H., 1990. Geophysical Research 76, 1113–1138.
Anticrack-associated faulting at very high pressure in Turcotte, D.L., Schubert, G., 1973. Frictional heating of the
natural olivine. Nature 348, 720–722. descending lithosphere. Journal of Geophysical Research
Hauck, S.A. II, Phillips, R.J., Hofmeister, A.M., 1999. Variable 78, 5876–5886.
conductivity: effects on the thermal structure of slabs. Turcotte, D.L., Schubert, G., 2002. Geodynamics, 2nd Edition.
Geophysical Research Letters 26, 3257–3260. Cambridge University Press, Cambridge, UK, 456pp.
Kirby, S.H., Durham, W.B., Stern, L.A., 1991. Mantle phase Wittaker, A., Bott, M.H.P., Waghorn, G.D., 1992. Stress and
changes and deep-earthquake faulting in subducting litho- plate boundary forces associated with subduction plate
sphere. Science 252, 216–225. margins. Journal of Geophysical Research 97, 11933–11944.