- Email: [email protected]
The thermal regime of a three-dimensional subducting lithospheric slab and its electromagnetic response: a numerical model
Tectonophysics,
35
225 (1993) 35-48
Elsevier Science Publishers B.V., Amsterdam
The thermal regime of a three-dimensional subducting lithospheric slab and its electromagnetic response: a numerical model F.W. Jones, F. Pascal and M.E. Ertman Department
of Physics, University of Alberta, Edmonton,
Alta. T6G .Zll, Canada
(Received February 5, 1992; revised version accepted February 11, 1993)
ABSTRACT A three-dimensional subducting lithospheric slab is simulated by a numerical technique and its temperature regime investigated. A single-slab model or models with multiple slab sections can be considered. An example is shown of two slab sections that dip at different angles. The near-surface heat flow above the slab is considered, and compared with a two-dimensional model. The temperature changes of particular points within the slab are monitored as it subducts and compared with the two-dimensional case. Considerable differences exist between the two- and three-dimensional cases, depending on the location of the point considered. An electromagnetic model is derived from the thermal model and the perturbations of time-varying electromagnetic fields caused by the electrical resistivity variations associated with the temperature variations are investigated.
Introduction De Jonge and Wortel, in a recent paper (1990), showed a three-dimensional temperature distribution below the Mediterranean area. Their temperatures were based on a combination of approximately 300 cross-sections that were derived from a two-dimensional downgoing lithospheric slab model originally developed by Minear and Toksoz (1970). De Jonge and Wortel (1990) used their thermal models to extend near-surface geophysical, geological and paleomagnetic observations to structures in the upper mantle by converting their derived temperature information to seismic p-velocities and so modelled the velocity structure. The de Jonge and Wortel (1990) approach is important because it gives information on the characteristics of the upper mantle that is independent of the usual deep seismic study methods. However, de Jonge and Wortel (1990) noted:
“A limitation of this method of constructing a 3D model is that lateral heat flow from one 2D section to another is ignored. This heat flow can be important if large along-trench varia0040-1951/93/$06.00
tions of geometry (strong curvature, transform type displacements) or temperature are present.” It is therefore of interest to investigate the thermal regime of a downgoing slab of limited along-trench extent and the lateral heat flow associated with such a three-dimensional configuration. In addition, it is of interest to explore the temperature variations of points on the three-dimensional slab as the slab subducts, and compare these variations with the two-dimensional case that is usually modelled. An alternative way to seismic methods for probing the nature of a subducting slab is by electromagnetic methods. It is possible to relate the electrical conductivity structure to the temperature regime, and by doing so a numerical electromagnetic induction model can be used to estimate the perturbations to surface fields caused by the electrical conductivity variations. The thermal model
Minear and Toksoz (1970) and Toksoz et al. (1971) investigated the temperature field and
0 1993 - Elsevier Science Publishers B.V. All rights reserved
36
F.W. JONES
geophysical effects of a downgoing slab using a two-dimensional numerical model. They calculated the temperature fields for various models, and included the effects of viscous dissipation, adiabatic compression, phase changes and radioactive heat generation and investigated the relative effects of different parameters. Following that, Toksoz and Hsui (1978) and Bodri and Bodri (1978) investigated back-arc convection and material flow induced by descending oceanic plates using numerical models. Hsui et al. (1983) further explored the effects of induced mantle flow on the melting of subducted oceanic crust. All these models were two-dimensional. It is therefore relevant to construct a three-dimensional numerical model of a downgoing lithospheric slab, both to study the effects at depth that concerned de Jonge and Wortel (1990), who used the original Minear and Toksoz (1970) model, and to provide a basis for investigations of the effects of convection and viscous dissipation and other finer effects. In the present work, we consider only heat conduction, but the model can be adapted to include these other effects. The downgoing slab model is shown in Figure 1. The slab moves in the -y-direction and begins subduction at the ocean-continent interface. The slab can subduct at different angles as it moves
Fig. 1. The three-dimensional
downgoing
slab model.
downward, and can consist of either a single slab or a slab composed of two parts that can dip downward at different angles as shown. The quasi-dynamic model is like that of Minear and Toksoz (1970) and Toksoz et al. (19711, but modified in the manner of Sydora et al. (1978) and applied over a grid of mesh points in three dimensions. The numerical method and the values of the various parameters used here are described in Appendix 1. The lithospheric slab translation method of Sydora et al. (1978) differs from that of Minear and Toksoz (19701 and led Sydora et al. (1978) to consider the subduction mechanism. Instead of bending of the slab by flexural stress with tensional stress at the top of the bend and compressional stress at the bottom of the bend implied by the Minear and Toksoz (1970) model, the Sydora et al. (1978) model implies a downward slipping of the end of the slab together with local tension as it starts to subduct due to a normal faulting mechanism. The mechanism leads to successive normal faults propagating into the slab, parallel to the first which permit its end to descend into the mantle. This model maintains a constant vertical thickness of the lithosphere as the plate moves downward in blocks by shearing along vertical planes. Stress calculations by Lliboutry
The slab moves
in the
-y-direction.
In the models
here,
x, = 200 km, x2 = 100 km, y = 591.6 km, y, = 217 km, I = 300 km, zt = 80 km, z, = 5 km. The model is superimposed of 61 x 61 x 61 = 226,981 grid points with Ax = 10 km, Ay = 9.86 km and AZ = 5.0 km. The time intervals At = 0.35 Ma. The slab velocity depends
ET AL.
on the variable
angles of subduction
x = 600 km, on a mesh
for slab translation
(Y, p and 4.
are
THERMAL
REGIME
OF A 3-D
SUBDUCTING
LITHOSPHERIC
SLAB
(1969), first motion studies by Isacks et al. (19681, and evidence of gravity faulting given by Katsumato and Sykes (19691, Malahof (1970), Ludwig et al. (19661, Stauder (1968a,b) and Lister (1971) support such a model. The approach here does not take into account the mantle wedge flow as described by Bodri and Bodri (1978) and Toksoz and Hsui (19781, and so the model is limited in this respect. The work of Bodri and Bodri (1978) showed that the mantle wedge flow can have a significant effect on surface heat flow, and this is not accounted for here. It is clear that to produce a more realistic model of the real earth, particularly in the region near the point of subduction, the method should be extended to include such effects. The intent of the present model is to investigate the limitations expressed by de Jonge and Wortel (1990) when they applied the two-dimensional model of Minear and Toksoz (1970) to a three-dimensional situation, and in particular to explore these effects at depth.
AND
ITS ELECTROMAGNETIC
31
RESPONSE
(4
X IN KU
Thermal results Examples of the lateral heat flow associated with the three-dimensional nature of the slab are shown in Figures 2 and 3. These figures are of the heat flow in the x-direction at time 7 Ma for slabs that have subducted for 7 Ma as described. It is seen that, because of the different assumed oceanic and continental geotherms, the temperatures of the downgoing slabs early in the subduction process are greater than the surroundings and heat initially flows from the slab material to the surroundings. However, for deeper portions of the slabs heat flows from the surroundings into the slabs. These figures are for time t = 7 Ma, and so the slabs have just reached the positions illustrated in Figures 2e and 3e. If the slab motion stops at this time, they will gradually approach equilibrium with the surrounding material. This approach toward equilibrium will be accompanied by a gradual decrease in the lateral temperature gradients and heat flows. At the instant illustrated in the figures, the lateral heat flows approach 100 mWmp2. Profiles in the y-direction of surface heat flow
Fig. 2. Heat flow values (mWm-*I in the x-z plane at five positions in the y-direction for the model: C#I = 27” (0 < t d 3.5 Ma); cu=27”, p=45”(3.5 Ma
F.W. JONES
38
density at a number of positions above the subducting three-dimensional slab are compared with a two-dimensional slab model in Figure 4 for
X IN Ktl
260
I
1
___
X IN
Ktl
Fig. 3. As in Fig. 2, but for C#J = 27” (0 < t d 3.5 Ma); a = 45”, p = 56” (3.5 Ma < f G 7 Ma).
ET AL.
cases with and without shear-strain heating. Figure 4a,b is for a model without shear-strain heating. It is seen that the surface heat flow density is reduced above the subducting slab on the continental side compared with the oceanic side consistent with the different oceanic and continental geotherms. In addition, a minimum in the heat flow occurs at the point where subduction begins, the location of the oceanic trench, and results from the fact that the isotherms are drawn downward by the motion of the slab material. This heat flow minimum is not evident in the models with shear-strain heating (Fig. 4c,d). Heat is generated in these latter models, and so the isotherms are not drawn downward in the same way by the slab motion. Figure 4b and d (at 14 Ma) indicate less abrupt transitions in the heat flow profiles from the continental to oceanic sides than do Figure 4a and c which were calculated at 7 Ma, just as subduction ceased. This illustrates that heat has flowed laterally in the y-direction, with gradual approach to equilibrium after the slab motion stopped, and consequent smoothing of the profiles. The effects of the limited extents of the three-dimensional models are evident in the profiles. It is seen that the vertical heat flow is affected by lateral flow of heat near the edges of the finite slabs. This is further illustrated in Figure 5 in which the detail of the profiles of Figure 4c near the point of subduction are shown. From this figure it is seen that the shapes of the vertical heat flow density profiles show considerable detail near the point of subduction for this case. This detail is related to the shear-strain heating, and does not appear in the cases without that heat source component. The detail gradually disappears with time, as can be seen in Figure 4d. Figure 5 further shows that the vertical heat flow density values for the three-dimensional models are reduced compared with the two-dimensional case, even toward the center of the slab. Figure 6 presents the variations of temperature with time of particular points on the composite slab, i.e. the variations of temperature with time of particular volumes of material as they subduct. These are compared with the 2-D slab model. In these plots, the slabs subduct for 7 Ma, and then the temperatures are monitored for an
THERMAL
REGIME
0
80
OF A 3-D
160
SUBDUCCING
340
LITHOSPHERIC
400
320
480
SLAB
560
AND
ITS ELECTROMAGNETIC
*,r/ , 0
!
80
39
RESPONSE
/
I
160
/
/
/
340
/
/
(
/
i
,
SE0
320
Yinkm
Yinkm
Fig. 4. (a) Vertical heat flow density at the surface along profiles in the y-direction at 7 Ma for the case with no shear-strain heating; 2-D = the two-dimensional model; 3-D side = near the edge of the 27”/27* slab 3-D center = at the center of the 27”/27” slab. (b) As in (a) but for 14 Ma (subduction ceased at 7 Maf. (c) As in (a), but for the case with shear strain heating of 1.6 x lo-’ Wm-” along the tops of the slabs and 1.6 X 10e6 Wm-” along the bottoms, ends and edges of the slabs. td) As in Cc)but for 14 Ma (subduction ceased at 7 Ma).
1400
1300
Legend -2-D -----3.Dside
800[ 0
------3-D center
I 2
I 4
1 6
I 8
I 10
I 12
I 14
Time in Ma
35 / 250
1
300
,
350
/
400
,
450
1
500
Yinkm Fig. 5. The profiles of Fig. 4c near the point where subduction begins; 2-D = two-dimensional model; 3-D side = threedimensional model near the edge of the slab; 3-D center = three-dimensional model at the center of the slab.
Fig. 6. Temperature variations with time of small volumes of material with initial positions y = 375 km, z = 45 km and which subduct for the 7 Ma and then are monitored for an additional 7 Ma with no slab movement. The 2-D case is indicated (2~0, and the positions in the x-direction are: (1) near the edge of the 27’/27O slab; (2) at the center of the 27”/27” slab; (3) near the inside edge of the 27”/27” slab; (4) near the inside edge of the 27”/4_5”slab, (5) at the center of the 27”/45” slab; (6) near the outside edge of the 27”/45” slab.
F.W. JONES
40
MODEL: 27.27-27.45 Temperature at the top of slabs at timeA0.50 Htp=1.6E+,Hbt=Hend=Hside=l.6E-6 2500
Ma
1
X=24*DX
X=3l’DX __
_.
Xz35’DX X=39’DX -.-.-
Depth in km (a)
MODEL 27.27-27.45 Temperature at the top of slabs at time=lO.!X Htp=4.8E-S,Hbt=Hend=Hside=4.8E-6
MU
2500 1 2000
0
/ X=20*0X ,......,.. I.,...
a
&
X=24*0X ---
1000
X=28*0X
P
X-Jl’OX I___--
‘Oar/ I
0
50
, 100
,
,
150
200
I_ 250
Eizf; 40
Al
Depth in km 03
Fig. 7. Temperatures along the slab surfaces for a composite slab that dips at 27” for 5.25 Ma and then separates into two parts that dip at 27” and 45” up to 10.5 Ma. Two-dimensional model results are included for comparison. (a> Shear-strain heating 1.6X lo-’ Wm-3 along the tops of the slabs, and 1.6 x 10m6 Wm-” along the bottoms, ends and edges of the slabs. (b) Shear-strain heating 4.8~ lo-’ Wme3 along the tops of the slabs, and 4.8~10~~ Wrnm3 along the bottoms, ends and edges of the slabs. The solidus for 1% dry basalt is indicated (geotherm). The values are taken at 10.5 Ma, the time that slab movement ceases. The legends indicate the various positions in the x-direction where the temperature profiles along the slab surfaces are taken (DX = 10 km, therefore, e.g., 2O.DX = position x = 200 km, 24’DX = position x = 240 km, etc).
ET AL.
additional 7 Ma with no slab movement taking place. It is seen that in all cases the temperatures initially decrease more than for the 2-D case until about 4 Ma. After 4 Ma the temperatures at different positions on the 27”-dipping slab differ to some extent from the 2-D case, but of course there are great differences between the 27”-dipping Z-D case and the part of the slab that subducts at 45”. It is interesting to consider the temperature variations along the tops of the slabs, and this is found to be a strong function of the amount of shear-strain heating. Figure 7 compares the slab surface temperatures for a composite slab model for which the two parts subducted together at 27” for 5.25 Ma and then separated and subducted at 27” and 45” for an additional 5.25 Ma for two different shear-strain heating values with a 27” do-dimensional model with the same shear-strain heating values. These temperature variations are compared with the solidus for 1% dry basalt (Wyllie, 1971, figs. 6-19). It is seen that the temperatures along the tops of the slabs that subduct at 27” throughout the motion do not greatly differ from the do-dimensional case. Of course, for the slab that bends and subducts at 45” during the later half of its motion, the temperatures are quite different. It is seen that the amount of shear-strain heating has a great effect on the slab surface temperatures. For the case with larger shear-strain heating (Fig. 7b), the slab surface temperatures are greater than the solidus for 1% dry basalt, and so it is possible that partial melt could be generated along the tops of the slabs in these cases. The electromagnetic model A three-dimensional numerical electromagnetic induction model developed by Jones and Pascoe (1972) and subsequently used by Lines and Jones (1973a,b) and others has been employed here to investigate the electromagnetic response of the electrical conductivity variations associated with the three-dimensional downgoing slab. The method has been described by Jones and Vozoff (1978). The model is a semi-infinite conductor in the region z > 0 with a plane
THERMAL REGIME OF A 3-D SUBDUCTING LI’THOSPHERIC: SLAB AND ITS ELECTROMAGNETIC
41
RESPONSE
boundary and which has regions of different electrical conductivity. For the downgoing slab model, the electrical ~ondu~tivities are derived from the temperature field through the relationship:
where CQ is the limiting ~ondu~tivi~ as T approaches infinity, E is the width of the energy gap for electronic conduction, k is Boltzman’s constant, and T is absolute temperature. The values chosen here for a0 and E are 10 ohm- ’ * cm-’ and 0.7 eV, respectively. We have found that the electrical ~onductivi~ is strongly dependent on the choice of E, and small changes in E result in large changes in the conductivity. Published values of E differ greatly, depending on the mode of conduction assumed and the materiaI considered (see, e.g., MacDonald, 1959; Rikitake, 1966). The values used here are reasonable for the model, and give electrical conductivities that are consistent with those given by Duba and Lilley (1972) in their examination of the effect of an ocean ridge on geomagnetic variations. The source field is taken to be uniform and oscillating with period 27r,/w which is sufficiently long that displacement currents can be ignored, and the magnetic permeability is taken as that of free space. Maxwell’s equations are then: VXH-aE V x E = -iol~,~H where the time factor exp(iwt> is understood in all field quantities and q is the conductivity appropriate to each region. Combining these equations Jeads to: V2E - V( V * E) = i7*E where $ = wpctocr, and this can be written as three scalar equations in Cartesian co-ordinates and solved sinlultaneously by a finite difference technique over a mesh of grid points superimposed over the model of interest. The resutt gives the three components of E, from which the magnetic field components can be calculated. The eie~troma~eti~ model here is constrncted with a surface conducting iayer 10 km thick with continental conductivity 10F3 ohm- ’ . m- 1 and
0
100
zN”
200 1592
300 0
/ 100
I 200
1 300
400
I 500
600
Y (km) Fig. 8. Sections showing the temperaturefield for the model: 4 = 27” (0 < t Q 3.5 Ma); 0: = 27”, fi = 45” (3.5 Ma < f G 7 Ma); o0 = 5 x IO’.ohmcm- ’ for oceanic lithosphere and subducting slab; w()= 10 ohmW’.cmm’ for continental lithosphere; E = 0.7 eV everywhere. (a) Through the center of the 27”/27” slab. Cb) Through the center of the 27”/41” slab. Temperatures in “C,
ocean conductivity 1 ohm-’ * m-’ overlying the electrical conductivity structure derived from the thermal model. The grid that is superimposed over the whole region consists of 38 x 38 x 38 (= 54,872) mesh points. The uniform source field oscihates in the x-direction and the model is solved iteratively as described by Lines and Jones (1973a,b).
The temperature field derived from a particular model is illustrated in Figure 8, This shows a case for which the slab is composed of two sections that subduct at different angles. The responses of the electromagnetic fields to the condu~tivity anomaly produced by the slab for this mode1 are shown in Figure 9. This figure gives contour plots of the amplitudes of the six field components. The electromagnetic fields are perturbed by the conductivity variations associated with the different values of a, and the tempera-
42
F.W. JONES
120 60 0
0 0
120
240
480
360
DISTANCE
600
LKMI
I
, 360
240
120
DISTRNCE
(a)
ET AL.
1
1,co
480
IKMJ
(a
660
540 480 Y
420
E
360
=, 300 g z
240
2
t
4
120 60
0
0-l
120
240
360
DISTANCE
480
600
0
/ 120
,
,
, 360
240
DISTRNCE
[KM)
, 480
I 600
iKtll
03 iZ” 660
1
600 540 480 f Y
420
E
360
g 300 =, 2
240
trj, , , , / / o!
0
1
I
120
,
/
/
240
I
360
DISTQNCE
I
I
480
,
I
600
0
IKMI
of the electric
and magnetic (period
240
360
DISTANCE
(cl
Fig. 9. Amplitudes
I20
480
/i”
(KM1
(f)
fields at the Earth’s
surface
for the model
of Fig. 8 at frequency
= 55.5 min). (a) E,. (b) E,. (cl E,. Cd) H,. (e) H,. (f) Hz.
= 0.0003
s-l
THERMAL
REGIME
OF A 3-D
SUBDWX’ING
LITHOSPHERIC
SLAB
AND
ITS ELECTROMAGNETIC
43
RESPONSE
720 -I
660
660 600
W
360
0 z
300
fi 0
240 180
120
120
60
60 0
0 I
1
0
I
I
120
I
I
240
I
360
DISTFlNCE
I
r 480
I
IjCIO
(KM)
b Fig. 10. In-phase (a) and quadrature (b) induction arrows for the model of Fig. 7.
ture variations across the region. The outline of the downgoing slab is clearly seen in the amplitude of E,, and the deflection of the induced electric currents by the slab is exhibited by the contour plot of E,. The amplitude of E, shows that vertical electric currents are caused by the finite extent of the slab. Although evidence for the difference in dip-angles for the two parts of the composite slab is not obvious, a careful look at the amplitudes of E, and H, (Fig. 9d) shows that differences exist between the shallower and deeper parts of the slabs. The magnetic field components also reflect the presence of the ocean-continent interface. Figure 10 illustrates the in-phase and quadrature induction arrows for the model of Figure 9. The source field is in the x-direction (E,), and we see that the in-phase induction arrows indicate that the current flow is concentrated in the downgoing slabs and is deflected into the upper part of the slab that dips more abruptly. Conclusions
A numerical three-dimensional downgoing slab model has been constructed and examples given for a number of models. It is found that substan-
tial lateral heat flow can result from the finite extent of the slabs. It is also seen that the temperature histories of the slab materials can be affected by the finite dimensions of the 3-D slab compared with a two-dimensional case. The electromagnetic model indicates that the slab perturbs the surface fields, and, depending on the model, the effect of the slab can be clearly seen, though effects due to differences in the dip-angles of composite slabs appear to be small. Acknowledgements
This work has been supported by the Natural Sciences and Engineering Research Council of Canada. Appendix 1. The parameters and finite difference scheme for the thermal model
Computation of the thermal regime is accomplished using the conservation of energy equation: +g=V.(KVT)+H
(1)
F.W.lONESETAL.
44
where C,. is the specific heat (taken here as 1.3 x lo3 Jkgg’K-‘), p is the density, T is temperature, K is thermal conductivity, and H is heat generation per unit volume. In Cartesian coordinates in three dimensions eq. (1) is: 2
2
$+$+;
qp(z,T);=K(T)
aK aT
3K aT
2
1
x, Y, 2)
= C + ( 16n3ST3)/(
x,
=P(z,
(2)
77~”+ 120~a,e-~‘~~)
Y, z> T)ga(z)
T(x,
+ 21.4) + 996
- 1184O/(P + 8.6) + 1340
3K aT
where C = lattice conductivity at low temperature n = index of refraction = 1.7; = 2.5 Wm-‘K-l; S = Stefan-Boltzman constant = 5.67 X lo-’ Wmw2Ke4; lg = opacity at low temperature = . . 10-i . m-‘; a0 = electrical conducttvtty; E= width of energy gap for electronic conduction; k = Boltzman’s constant = 1.38 X 10-23JK-‘. In what follows, we will write K(T) = K, with the understanding that it depends on temperature. However, the temperature dependence of the lattice conductivity (C) is ignored. The density distribution follows the Bullen A density curve (see Sydora et al., 1978) and the heat generation term includes adiabatic and shear-strain heating as: H(T,
T = 7.99(P + 21.4)
Oceanic: T = 4.34(P + 8.6)
where the spatial and temperature dependences of K, p and H are indicated. The thermal conductivity is taken as in MacDonald (1959) and includes a constant lattice conductivity term and a radiative transfer term that depends on temperature: K(T)
Continental:
- 24914/(P
+;tllax+--+araz ay a~ +H(T,
The initial vertical temperature condition is from the relationships given by Mercier and Carter (1975):
Y, z)“z+Hs,,
where g = 9.81 m/ss2; a(z) = exp(3.58 0.0072~) x lo-” as in Toksoz et al. (1971) (z measured in km); u, = the vertical component of velocity of the downgoing slab; Hsh = the shearstrain heating, confined to 10 km thick layers along the tops, bottoms, sides and ends of the slabs, and can be adjusted according to the model desired.
where P = 0.31~ kbar km-‘, and the temperature over the surface of the Earth is held constant and equal to 0°C throughout the whole calculation, i.e., T(x,
y, 0, t) =0
t&O
The other boundary conditions are (where L,, L, represent the limits of the mesh in the X-, y- and z-directions, respectively). L,,
aT a2
=E.F=13.5mWm-’ r=L,
K’
aT aT ==0 ax Ix=o ax 1x=L,
aT aT =0 ay Iy=o ay Iy=L,= The solution proceeds by translating the slab, and then solving eq. (2) alternately in each of the three directions by the Douglas alternating direction implicit (ADI) method (Peaceman and Rachford, 1955). Using the notation: xi=(i-1)Ax
i=1,2,...,N,(N,-l)Ax=L,
y,=(j--l)Ay
j=1,2,...,N,(N,-l)Ay=L,
z,=(k-1)Az
k=1,2,...,N,(N,-l)Az=L,
T(x;, y,, zk, t> = T (i, j, k) for a fixed time t, and we consider the x-, y-, and z-directions as follows:
x-direction
Let T(i, j, k) be the temperature at t = to, and let T,,(i, j, k) be the temperature at t = to + At.
THERMAL
REGIME
OF A 3-D
SUBDUCTING
Then, the first estimate of TJi, T,(i, j, k) is given by: cIp[z,
T(i,
+H(i,
SLAB
a2T
1
aK aT +--+zz ay ay
1 aK
AyAz
- y ZAyAz
Ax
Ax
aK aT =R,T
j,k)+T,(i-1,
j, k)
T,(i + 1, j, k) - T,(i - 1, j, k)
i,k)= At -CJdz, f)T,(i, 2Ax
1
by
1
R(T)
(5)
where R,(T) = 2 Au R(T). If i = 2, then aT/ax 1X=o= 0 implies that T(1, j, k)= T(2, j, k) and also, T,(l, j, k)= T,(2, j, k).
Equation (5) then becomes: K-
-
AyAz
1 aK +?axAyAz+
2C,,P(ZJ) At
xT,(Z j,k)
I (4)
AX
+ 2 ~AYAz [T(i+l,
j, k) -2T(i,
j, k)
- 1, j, k)]
Av
1
1
1 aK
where R(T) =
+T(i
1 aK
+ ?~AYAz
(3)
(Ax)’
+zax
45
RESPONSE
Multiplying eq. (4) by 2AxAyAz = 2Au gives: K-
j, k, T)
j,k)-2T,(i,
1 aK
ITS ELECTROMAGNETIC
j, k), denoted by
If the derivatives in eq. (3) are replaced differences we have: T,(i+l,
AND
~,(i, j, k) - T(i, j, k)
j, k)]
a2T +K -+az2 [ a?
LITHOSPHERK
T,(3, j, k)
=R,(T) where in R,(T) we replace T(1, j, k) by T(2, j, k).
+--: E&[T(i+
+K &
1, j, k) - T(i-
[T(i,
j + 1, k) - 2T(i,
1, j, k)]
j, k)
If i =N, - 1, then aT/ax IXzL, = 0 implies RN, - 1, j, k) = TCN,, j, k) and T,(N, - 1, j, k) = T,(N,, j, k) and eq. (5) becomes: AyAz
+T(i,
K-
j - 1, k)]
&T(’
Ax
T,( N, - 2, j, k)
1 aK + z %AyAz
I, j, k + 1) - 2T( i, j, k)
+K(AZ)’
1
1 aK
- 2 zAyAz
+ T( i, j, k - l)] -2C,p(z,
aK 1 +--ay 2Ay [T(i, j + 1, k) - T(i, j - 1, k)] aK 1 +-a.2 2AzlT( +C,p(z,
i, j, k + 1) - T(i,
T)&T(i, i,k) +H(T)
j, k - 1)]
1
T);
I
T,(N,-
where in R,(T) as given j, k) by ZW, - I, j, k). Equation (4) therefore 2, NY- 1 and k = 2, N, equations with unknowns N,-
1.
1, j, k) =R,(T)
above we replace T(N,, represents for each j = - 1 a system of N, - 2 T,(i, j, k) where i = 2,
F.W.JONESETAL.
46
Note that the coefficients of this system contain the terms aK/ax, aK/ay, and aK/a.z, which may be computed by the formulas:
where R,(T)
Cd(T)
= -2
At
AvT,(i, j, k)
K(i + 1, j, k) - K(i - 1, j, k)
aK ax xi. y,, zq =
26x
DAY
ays z% Y,>Lk = aK
-
-_
-2T(i,
K(i, j, k + 1) - K(i, j, k - 1)
--
2h.z
a2 x,7Yj,ZA
y-direction
C,.p[ z,T(i,
j, k)]
K
a2T, =_+2 I ax2
-
2
I
j, k) + T(i, j - 1, k)] j+ 1, k)
-AxAz[T(i, 2 3Y j-
1, k)]
For j = 2, we replace in eq. (7) T,G, 1, k) by T2 (i, 2, k) and for j = I?, - 1, T,(i, NY, k) = T&i,
N, - 1, k). Equation (7) then represents, for each i = 2, N, - 1, and k = 2, N, - 1, a system with unknowns T,(i, j, k), j = 2, & - 1.
T2(i, i, k) - T( i, j, k) At
z-direction
a2T ax* I
ay2 +ayl
---+--
The third and final estimate for T,,(i, j, k), denoted by T&i, j, k) is now obtained. The operation in the z-direction is treated like that in the y-direction. Writing eq. (6) as:
1
a2T
a2T2
K
j+ 1, k)
ale
i
-T(i,
A second estimate of T,,(i, j, kl is now obtained and denoted by T,(i, j, k). In eq, (31, the terms that contain the derivatives of T with respect to y can be separated into two parts and we obtain:
f-
F[T(i,
-K
K(i, j + 1, k) - K(i, j - 1, k)
aK
I aK aT
C,.p[z,T(i,j, k)]
T,(i, j, k) - T(i, j, k) At
ay ay I aK aT
a*T
+Ka,2+
zz
+H(i, i, k, T)
(6)
Subtracting eq. (6) from eq. (3), replacing the derivatives with respect to p by differences, and rearranging gives: AxAz 1 aK - --AxAz KAY
2 ay
I
T2(i, j-
1 +y”
1, k)
i
aT,2 -+g7 az2
a2T I
I
=&(T)
17)
+H(i,j,k,T)
and subtracting eq. (8) from eq. (61, replacing the derivatives with respect to z by differences, and
THERMAL
REGIME
rearranging AxAyAz): AxAy
K-
AZ
OF A 3-D
terms
SUBDUCTING
gives
1 aK
- z ,zAxAy
I
LITHOSPHERIC
(again
SLAB
with
Au =
T3(i, j, k - 1)
j,k) 1T3(i,
1
1 aK
+ ?GAXAY
T3(i,
where At
x [T(i,
AuT2(i, j, k) - KF
j, k + 1) - 2T(i,
j, k)
+ T( i, j, k - l)]
1 aK --- 2 az AxAy[T(i, -T(i,
j, k-
j, k + 1)
l)]
Because T(x, y, 0, t ) = 0 for any t, we have T3 (i, j, 1) = 0.
From the boundary have:
condition
at z = L,,
=R,(T) andfork=N,-1: 1 aK
- 2 %AxAy
1
T3(i, j, N, - 2)
1
P(T)
+ 2C,;?Av
XT3(i, i, N, - 1) =R3(T)
-
T(i,
1, k) = T(i, 2, k),
T(i,
N,,, k) = T(i,
N,, - 1, k)
T(i,
j, N,) = T(i,
j, N, - 1) + 7
AxAy
KF
Az.F
The solution therefore proceeds by simulating the slab movement by translating the slab temperatures through the mesh and subsequently applying the above AD1 method to solve eq. (1). The numerical mesh must be such as to accommodate the required rate and angle of subduction. Movement of the slab can cease, with subsequent warming if desired. References
Equation (9) for k = 2 becomes:
AxAy
j, k) = T(N, - 1, j, k)
we
F T3(i, j, N,) = T3( i, j, N, - 1) + AZ. z
K - AZ
T(l, j, k) = T(2, j, k),
AxAy
C,P(Z~)
41
RESPONSE
Equation (9) for each i = 2, N, - 1 and j = 2, N,, - 1 represents a system of equations with unknowns T&i, j, k), k = 2, N, - 1. After completing the computations in the three directions, the temperatures at the boundaries are computed according to the boundary conditions:
T(N,,
(9)
= -2
ITS ELECTROMAGNETIC
j, k + 1)
=RdT)
R,(T)
AND
F
1 aK
+ ?%AxAy
%Az. I
Bodri, L. and Bodri, B., 1978. Numerical investigation of tectonic flow in island-arc areas. Tectonophysics, 50: 163115. De Jonge, M.R. and Wortel, M.J.R., 1990. The thermal structure of the Mediterranean upper mantle: a forward modelling approach. Terra Nova, 2: 609-616. Duba, A. and Lilley, F.E.M., 1972. Effect of an ocean ridge model on geomagnetic variations. J. Geophys. Res., 72: 7100-7105. Hsui, A.T., Marsh, B.D. and Toksoz, M.N., 1983. On melting of the subducting oceanic crust: effects of subduction induced mantle flow. Tectonophysics, 99: 207-220. Isacks, B., Oliver, J., and Sykes, L.R., 1968. Seismology and the new global tectonics. J. Geophys. Res., 73: 5855-5899. Jones, F.W. and Pascoe, L.J., 1972. The perturbation of alternating geomagnetic fields by three-dimensional conductivity inhomogeneities. Geophys. J.R. Astron. Sot., 27: 479-485. Jones, F.W. and Vozoff, K., 1978. The calculation of magnetotelluric quantities for three-dimensional conductivity inhomogeneities. Geophysics, 43: 1167-l 175. Katsumata, M. and Sykes, L.R., 1969. Seismicity and tectonics of the western Pacific, Izu-Mariana, Caroline and Ryukya-Taiwan regions. J. Geophys. Res., 74: 5923-5948.
48 Lines, L.R. and Jones, F.W., 1973a. The perturbation of alternating geomagnetic fields by three-dimensional island structures. Geophys. J.R. Astron. Sot., 32: 133-154. Lines, L.R. and Jones, F.W., 1973b. The perturbation of alternating geomagnetic fields by an island near a coastline. Can. J. Earth Sci., IO: 510-518. Lister, C.R.B., 1971. Tectonic movement in the Chile trench. Science, 173: 719-722. Lliboutry, L., 1969. Sea-floor spreading, continental drift and lithosphere sinking with an asthenosphere at melting point. J. Geophys. Res., 74: 6525-6540. Ludwig, W.J., Ewing, J.E., Ewing, M. Murauchi, S., Den, N., Asano, S., Hotta, H., Hayakawa, M., Asanuma, T., Ichikawa, K. and Noguchi, I., 1966. Sediments and structure of the Japan trench. J. Geophys. Res., 71: 2121-2137. MacDonald, G.J.F., 1959. Calculations on the thermal history of the Earth. J. Geophys. Res., 64: 1967-2~. Malahof, A., 1970. Some possible mechanisms for gravity and thrust faults under ocean trenches. J. Geophys. Res., 75: 1992-2001. Mercier, J.C. and Carter N.L., 1975. Pyroxene geotherms. J. Geophys. Res., 80: 3349-3362. Minear, J.W. and Toksoz, M.N., 1970. Thermal regime of a downgoing slab and new global tectonics. J. Geophys. Res., 75: 1397-1419.
F.W. JONES
ET AL.
Peaceman, D.W. and Rachford, Jr., H.H., 1955. The numerical solution of parabolic and elliptic differential equations. J. Sot. Ind. Appl. Math., 3: 28-41. Rikitake, T., 1966. Electromagnetism and the Earth’s Interior, Elsevier, Amsterdam, 308 pp. Stauder, W., 1968a. Mechanism of the Rat Island earthquake sequence of February 4, 1965, with relation to island arcs and sea floor spreading. J. Geophys. Res., 73: 3847-3858. Stauder, W., 1968b. Tensional character of earthquake foci beneath the Aleutian trench with relation to sea floor spreading. J. Geophys. Res., 73: 7693-7701. Sydora, L.J., Jones, F.W. and Lambert, R. St.J., 1978. The thermal regime of the descending lithosphere: The effect of varying angle and rate of subduction. Can. J. Earth Sci., 154: 626-641. Toksoz, M.N. and Hsui, A.T., 1978. Numerical studies of back-arc convection and the formation of marginal basins. Tectonophysics, 50: 177-196. Toksoz, M.N., Minear, J.W. and Julian, B.R., 1971. Temperature field and geophysical effects of a downgoing slab. J. Geophys. Res., 76: 1113-1138. Wyllie, P., 1971. The Dynamic Earth. Wiley, New York, NY, 416 pp.
35
225 (1993) 35-48
Elsevier Science Publishers B.V., Amsterdam
The thermal regime of a three-dimensional subducting lithospheric slab and its electromagnetic response: a numerical model F.W. Jones, F. Pascal and M.E. Ertman Department
of Physics, University of Alberta, Edmonton,
Alta. T6G .Zll, Canada
(Received February 5, 1992; revised version accepted February 11, 1993)
ABSTRACT A three-dimensional subducting lithospheric slab is simulated by a numerical technique and its temperature regime investigated. A single-slab model or models with multiple slab sections can be considered. An example is shown of two slab sections that dip at different angles. The near-surface heat flow above the slab is considered, and compared with a two-dimensional model. The temperature changes of particular points within the slab are monitored as it subducts and compared with the two-dimensional case. Considerable differences exist between the two- and three-dimensional cases, depending on the location of the point considered. An electromagnetic model is derived from the thermal model and the perturbations of time-varying electromagnetic fields caused by the electrical resistivity variations associated with the temperature variations are investigated.
Introduction De Jonge and Wortel, in a recent paper (1990), showed a three-dimensional temperature distribution below the Mediterranean area. Their temperatures were based on a combination of approximately 300 cross-sections that were derived from a two-dimensional downgoing lithospheric slab model originally developed by Minear and Toksoz (1970). De Jonge and Wortel (1990) used their thermal models to extend near-surface geophysical, geological and paleomagnetic observations to structures in the upper mantle by converting their derived temperature information to seismic p-velocities and so modelled the velocity structure. The de Jonge and Wortel (1990) approach is important because it gives information on the characteristics of the upper mantle that is independent of the usual deep seismic study methods. However, de Jonge and Wortel (1990) noted:
“A limitation of this method of constructing a 3D model is that lateral heat flow from one 2D section to another is ignored. This heat flow can be important if large along-trench varia0040-1951/93/$06.00
tions of geometry (strong curvature, transform type displacements) or temperature are present.” It is therefore of interest to investigate the thermal regime of a downgoing slab of limited along-trench extent and the lateral heat flow associated with such a three-dimensional configuration. In addition, it is of interest to explore the temperature variations of points on the three-dimensional slab as the slab subducts, and compare these variations with the two-dimensional case that is usually modelled. An alternative way to seismic methods for probing the nature of a subducting slab is by electromagnetic methods. It is possible to relate the electrical conductivity structure to the temperature regime, and by doing so a numerical electromagnetic induction model can be used to estimate the perturbations to surface fields caused by the electrical conductivity variations. The thermal model
Minear and Toksoz (1970) and Toksoz et al. (1971) investigated the temperature field and
0 1993 - Elsevier Science Publishers B.V. All rights reserved
36
F.W. JONES
geophysical effects of a downgoing slab using a two-dimensional numerical model. They calculated the temperature fields for various models, and included the effects of viscous dissipation, adiabatic compression, phase changes and radioactive heat generation and investigated the relative effects of different parameters. Following that, Toksoz and Hsui (1978) and Bodri and Bodri (1978) investigated back-arc convection and material flow induced by descending oceanic plates using numerical models. Hsui et al. (1983) further explored the effects of induced mantle flow on the melting of subducted oceanic crust. All these models were two-dimensional. It is therefore relevant to construct a three-dimensional numerical model of a downgoing lithospheric slab, both to study the effects at depth that concerned de Jonge and Wortel (1990), who used the original Minear and Toksoz (1970) model, and to provide a basis for investigations of the effects of convection and viscous dissipation and other finer effects. In the present work, we consider only heat conduction, but the model can be adapted to include these other effects. The downgoing slab model is shown in Figure 1. The slab moves in the -y-direction and begins subduction at the ocean-continent interface. The slab can subduct at different angles as it moves
Fig. 1. The three-dimensional
downgoing
slab model.
downward, and can consist of either a single slab or a slab composed of two parts that can dip downward at different angles as shown. The quasi-dynamic model is like that of Minear and Toksoz (1970) and Toksoz et al. (19711, but modified in the manner of Sydora et al. (1978) and applied over a grid of mesh points in three dimensions. The numerical method and the values of the various parameters used here are described in Appendix 1. The lithospheric slab translation method of Sydora et al. (1978) differs from that of Minear and Toksoz (19701 and led Sydora et al. (1978) to consider the subduction mechanism. Instead of bending of the slab by flexural stress with tensional stress at the top of the bend and compressional stress at the bottom of the bend implied by the Minear and Toksoz (1970) model, the Sydora et al. (1978) model implies a downward slipping of the end of the slab together with local tension as it starts to subduct due to a normal faulting mechanism. The mechanism leads to successive normal faults propagating into the slab, parallel to the first which permit its end to descend into the mantle. This model maintains a constant vertical thickness of the lithosphere as the plate moves downward in blocks by shearing along vertical planes. Stress calculations by Lliboutry
The slab moves
in the
-y-direction.
In the models
here,
x, = 200 km, x2 = 100 km, y = 591.6 km, y, = 217 km, I = 300 km, zt = 80 km, z, = 5 km. The model is superimposed of 61 x 61 x 61 = 226,981 grid points with Ax = 10 km, Ay = 9.86 km and AZ = 5.0 km. The time intervals At = 0.35 Ma. The slab velocity depends
ET AL.
on the variable
angles of subduction
x = 600 km, on a mesh
for slab translation
(Y, p and 4.
are
THERMAL
REGIME
OF A 3-D
SUBDUCTING
LITHOSPHERIC
SLAB
(1969), first motion studies by Isacks et al. (19681, and evidence of gravity faulting given by Katsumato and Sykes (19691, Malahof (1970), Ludwig et al. (19661, Stauder (1968a,b) and Lister (1971) support such a model. The approach here does not take into account the mantle wedge flow as described by Bodri and Bodri (1978) and Toksoz and Hsui (19781, and so the model is limited in this respect. The work of Bodri and Bodri (1978) showed that the mantle wedge flow can have a significant effect on surface heat flow, and this is not accounted for here. It is clear that to produce a more realistic model of the real earth, particularly in the region near the point of subduction, the method should be extended to include such effects. The intent of the present model is to investigate the limitations expressed by de Jonge and Wortel (1990) when they applied the two-dimensional model of Minear and Toksoz (1970) to a three-dimensional situation, and in particular to explore these effects at depth.
AND
ITS ELECTROMAGNETIC
31
RESPONSE
(4
X IN KU
Thermal results Examples of the lateral heat flow associated with the three-dimensional nature of the slab are shown in Figures 2 and 3. These figures are of the heat flow in the x-direction at time 7 Ma for slabs that have subducted for 7 Ma as described. It is seen that, because of the different assumed oceanic and continental geotherms, the temperatures of the downgoing slabs early in the subduction process are greater than the surroundings and heat initially flows from the slab material to the surroundings. However, for deeper portions of the slabs heat flows from the surroundings into the slabs. These figures are for time t = 7 Ma, and so the slabs have just reached the positions illustrated in Figures 2e and 3e. If the slab motion stops at this time, they will gradually approach equilibrium with the surrounding material. This approach toward equilibrium will be accompanied by a gradual decrease in the lateral temperature gradients and heat flows. At the instant illustrated in the figures, the lateral heat flows approach 100 mWmp2. Profiles in the y-direction of surface heat flow
Fig. 2. Heat flow values (mWm-*I in the x-z plane at five positions in the y-direction for the model: C#I = 27” (0 < t d 3.5 Ma); cu=27”, p=45”(3.5 Ma
F.W. JONES
38
density at a number of positions above the subducting three-dimensional slab are compared with a two-dimensional slab model in Figure 4 for
X IN Ktl
260
I
1
___
X IN
Ktl
Fig. 3. As in Fig. 2, but for C#J = 27” (0 < t d 3.5 Ma); a = 45”, p = 56” (3.5 Ma < f G 7 Ma).
ET AL.
cases with and without shear-strain heating. Figure 4a,b is for a model without shear-strain heating. It is seen that the surface heat flow density is reduced above the subducting slab on the continental side compared with the oceanic side consistent with the different oceanic and continental geotherms. In addition, a minimum in the heat flow occurs at the point where subduction begins, the location of the oceanic trench, and results from the fact that the isotherms are drawn downward by the motion of the slab material. This heat flow minimum is not evident in the models with shear-strain heating (Fig. 4c,d). Heat is generated in these latter models, and so the isotherms are not drawn downward in the same way by the slab motion. Figure 4b and d (at 14 Ma) indicate less abrupt transitions in the heat flow profiles from the continental to oceanic sides than do Figure 4a and c which were calculated at 7 Ma, just as subduction ceased. This illustrates that heat has flowed laterally in the y-direction, with gradual approach to equilibrium after the slab motion stopped, and consequent smoothing of the profiles. The effects of the limited extents of the three-dimensional models are evident in the profiles. It is seen that the vertical heat flow is affected by lateral flow of heat near the edges of the finite slabs. This is further illustrated in Figure 5 in which the detail of the profiles of Figure 4c near the point of subduction are shown. From this figure it is seen that the shapes of the vertical heat flow density profiles show considerable detail near the point of subduction for this case. This detail is related to the shear-strain heating, and does not appear in the cases without that heat source component. The detail gradually disappears with time, as can be seen in Figure 4d. Figure 5 further shows that the vertical heat flow density values for the three-dimensional models are reduced compared with the two-dimensional case, even toward the center of the slab. Figure 6 presents the variations of temperature with time of particular points on the composite slab, i.e. the variations of temperature with time of particular volumes of material as they subduct. These are compared with the 2-D slab model. In these plots, the slabs subduct for 7 Ma, and then the temperatures are monitored for an
THERMAL
REGIME
0
80
OF A 3-D
160
SUBDUCCING
340
LITHOSPHERIC
400
320
480
SLAB
560
AND
ITS ELECTROMAGNETIC
*,r/ , 0
!
80
39
RESPONSE
/
I
160
/
/
/
340
/
/
(
/
i
,
SE0
320
Yinkm
Yinkm
Fig. 4. (a) Vertical heat flow density at the surface along profiles in the y-direction at 7 Ma for the case with no shear-strain heating; 2-D = the two-dimensional model; 3-D side = near the edge of the 27”/27* slab 3-D center = at the center of the 27”/27” slab. (b) As in (a) but for 14 Ma (subduction ceased at 7 Maf. (c) As in (a), but for the case with shear strain heating of 1.6 x lo-’ Wm-” along the tops of the slabs and 1.6 X 10e6 Wm-” along the bottoms, ends and edges of the slabs. td) As in Cc)but for 14 Ma (subduction ceased at 7 Ma).
1400
1300
Legend -2-D -----3.Dside
800[ 0
------3-D center
I 2
I 4
1 6
I 8
I 10
I 12
I 14
Time in Ma
35 / 250
1
300
,
350
/
400
,
450
1
500
Yinkm Fig. 5. The profiles of Fig. 4c near the point where subduction begins; 2-D = two-dimensional model; 3-D side = threedimensional model near the edge of the slab; 3-D center = three-dimensional model at the center of the slab.
Fig. 6. Temperature variations with time of small volumes of material with initial positions y = 375 km, z = 45 km and which subduct for the 7 Ma and then are monitored for an additional 7 Ma with no slab movement. The 2-D case is indicated (2~0, and the positions in the x-direction are: (1) near the edge of the 27’/27O slab; (2) at the center of the 27”/27” slab; (3) near the inside edge of the 27”/27” slab; (4) near the inside edge of the 27”/4_5”slab, (5) at the center of the 27”/45” slab; (6) near the outside edge of the 27”/45” slab.
F.W. JONES
40
MODEL: 27.27-27.45 Temperature at the top of slabs at timeA0.50 Htp=1.6E+,Hbt=Hend=Hside=l.6E-6 2500
Ma
1
X=24*DX
X=3l’DX __
_.
Xz35’DX X=39’DX -.-.-
Depth in km (a)
MODEL 27.27-27.45 Temperature at the top of slabs at time=lO.!X Htp=4.8E-S,Hbt=Hend=Hside=4.8E-6
MU
2500 1 2000
0
/ X=20*0X ,......,.. I.,...
a
&
X=24*0X ---
1000
X=28*0X
P
X-Jl’OX I___--
‘Oar/ I
0
50
, 100
,
,
150
200
I_ 250
Eizf; 40
Al
Depth in km 03
Fig. 7. Temperatures along the slab surfaces for a composite slab that dips at 27” for 5.25 Ma and then separates into two parts that dip at 27” and 45” up to 10.5 Ma. Two-dimensional model results are included for comparison. (a> Shear-strain heating 1.6X lo-’ Wm-3 along the tops of the slabs, and 1.6 x 10m6 Wm-” along the bottoms, ends and edges of the slabs. (b) Shear-strain heating 4.8~ lo-’ Wme3 along the tops of the slabs, and 4.8~10~~ Wrnm3 along the bottoms, ends and edges of the slabs. The solidus for 1% dry basalt is indicated (geotherm). The values are taken at 10.5 Ma, the time that slab movement ceases. The legends indicate the various positions in the x-direction where the temperature profiles along the slab surfaces are taken (DX = 10 km, therefore, e.g., 2O.DX = position x = 200 km, 24’DX = position x = 240 km, etc).
ET AL.
additional 7 Ma with no slab movement taking place. It is seen that in all cases the temperatures initially decrease more than for the 2-D case until about 4 Ma. After 4 Ma the temperatures at different positions on the 27”-dipping slab differ to some extent from the 2-D case, but of course there are great differences between the 27”-dipping Z-D case and the part of the slab that subducts at 45”. It is interesting to consider the temperature variations along the tops of the slabs, and this is found to be a strong function of the amount of shear-strain heating. Figure 7 compares the slab surface temperatures for a composite slab model for which the two parts subducted together at 27” for 5.25 Ma and then separated and subducted at 27” and 45” for an additional 5.25 Ma for two different shear-strain heating values with a 27” do-dimensional model with the same shear-strain heating values. These temperature variations are compared with the solidus for 1% dry basalt (Wyllie, 1971, figs. 6-19). It is seen that the temperatures along the tops of the slabs that subduct at 27” throughout the motion do not greatly differ from the do-dimensional case. Of course, for the slab that bends and subducts at 45” during the later half of its motion, the temperatures are quite different. It is seen that the amount of shear-strain heating has a great effect on the slab surface temperatures. For the case with larger shear-strain heating (Fig. 7b), the slab surface temperatures are greater than the solidus for 1% dry basalt, and so it is possible that partial melt could be generated along the tops of the slabs in these cases. The electromagnetic model A three-dimensional numerical electromagnetic induction model developed by Jones and Pascoe (1972) and subsequently used by Lines and Jones (1973a,b) and others has been employed here to investigate the electromagnetic response of the electrical conductivity variations associated with the three-dimensional downgoing slab. The method has been described by Jones and Vozoff (1978). The model is a semi-infinite conductor in the region z > 0 with a plane
THERMAL REGIME OF A 3-D SUBDUCTING LI’THOSPHERIC: SLAB AND ITS ELECTROMAGNETIC
41
RESPONSE
boundary and which has regions of different electrical conductivity. For the downgoing slab model, the electrical ~ondu~tivities are derived from the temperature field through the relationship:
where CQ is the limiting ~ondu~tivi~ as T approaches infinity, E is the width of the energy gap for electronic conduction, k is Boltzman’s constant, and T is absolute temperature. The values chosen here for a0 and E are 10 ohm- ’ * cm-’ and 0.7 eV, respectively. We have found that the electrical ~onductivi~ is strongly dependent on the choice of E, and small changes in E result in large changes in the conductivity. Published values of E differ greatly, depending on the mode of conduction assumed and the materiaI considered (see, e.g., MacDonald, 1959; Rikitake, 1966). The values used here are reasonable for the model, and give electrical conductivities that are consistent with those given by Duba and Lilley (1972) in their examination of the effect of an ocean ridge on geomagnetic variations. The source field is taken to be uniform and oscillating with period 27r,/w which is sufficiently long that displacement currents can be ignored, and the magnetic permeability is taken as that of free space. Maxwell’s equations are then: VXH-aE V x E = -iol~,~H where the time factor exp(iwt> is understood in all field quantities and q is the conductivity appropriate to each region. Combining these equations Jeads to: V2E - V( V * E) = i7*E where $ = wpctocr, and this can be written as three scalar equations in Cartesian co-ordinates and solved sinlultaneously by a finite difference technique over a mesh of grid points superimposed over the model of interest. The resutt gives the three components of E, from which the magnetic field components can be calculated. The eie~troma~eti~ model here is constrncted with a surface conducting iayer 10 km thick with continental conductivity 10F3 ohm- ’ . m- 1 and
0
100
zN”
200 1592
300 0
/ 100
I 200
1 300
400
I 500
600
Y (km) Fig. 8. Sections showing the temperaturefield for the model: 4 = 27” (0 < t Q 3.5 Ma); 0: = 27”, fi = 45” (3.5 Ma < f G 7 Ma); o0 = 5 x IO’.ohmcm- ’ for oceanic lithosphere and subducting slab; w()= 10 ohmW’.cmm’ for continental lithosphere; E = 0.7 eV everywhere. (a) Through the center of the 27”/27” slab. Cb) Through the center of the 27”/41” slab. Temperatures in “C,
ocean conductivity 1 ohm-’ * m-’ overlying the electrical conductivity structure derived from the thermal model. The grid that is superimposed over the whole region consists of 38 x 38 x 38 (= 54,872) mesh points. The uniform source field oscihates in the x-direction and the model is solved iteratively as described by Lines and Jones (1973a,b).
The temperature field derived from a particular model is illustrated in Figure 8, This shows a case for which the slab is composed of two sections that subduct at different angles. The responses of the electromagnetic fields to the condu~tivity anomaly produced by the slab for this mode1 are shown in Figure 9. This figure gives contour plots of the amplitudes of the six field components. The electromagnetic fields are perturbed by the conductivity variations associated with the different values of a, and the tempera-
42
F.W. JONES
120 60 0
0 0
120
240
480
360
DISTANCE
600
LKMI
I
, 360
240
120
DISTRNCE
(a)
ET AL.
1
1,co
480
IKMJ
(a
660
540 480 Y
420
E
360
=, 300 g z
240
2
t
4
120 60
0
0-l
120
240
360
DISTANCE
480
600
0
/ 120
,
,
, 360
240
DISTRNCE
[KM)
, 480
I 600
iKtll
03 iZ” 660
1
600 540 480 f Y
420
E
360
g 300 =, 2
240
trj, , , , / / o!
0
1
I
120
,
/
/
240
I
360
DISTQNCE
I
I
480
,
I
600
0
IKMI
of the electric
and magnetic (period
240
360
DISTANCE
(cl
Fig. 9. Amplitudes
I20
480
/i”
(KM1
(f)
fields at the Earth’s
surface
for the model
of Fig. 8 at frequency
= 55.5 min). (a) E,. (b) E,. (cl E,. Cd) H,. (e) H,. (f) Hz.
= 0.0003
s-l
THERMAL
REGIME
OF A 3-D
SUBDWX’ING
LITHOSPHERIC
SLAB
AND
ITS ELECTROMAGNETIC
43
RESPONSE
720 -I
660
660 600
W
360
0 z
300
fi 0
240 180
120
120
60
60 0
0 I
1
0
I
I
120
I
I
240
I
360
DISTFlNCE
I
r 480
I
IjCIO
(KM)
b Fig. 10. In-phase (a) and quadrature (b) induction arrows for the model of Fig. 7.
ture variations across the region. The outline of the downgoing slab is clearly seen in the amplitude of E,, and the deflection of the induced electric currents by the slab is exhibited by the contour plot of E,. The amplitude of E, shows that vertical electric currents are caused by the finite extent of the slab. Although evidence for the difference in dip-angles for the two parts of the composite slab is not obvious, a careful look at the amplitudes of E, and H, (Fig. 9d) shows that differences exist between the shallower and deeper parts of the slabs. The magnetic field components also reflect the presence of the ocean-continent interface. Figure 10 illustrates the in-phase and quadrature induction arrows for the model of Figure 9. The source field is in the x-direction (E,), and we see that the in-phase induction arrows indicate that the current flow is concentrated in the downgoing slabs and is deflected into the upper part of the slab that dips more abruptly. Conclusions
A numerical three-dimensional downgoing slab model has been constructed and examples given for a number of models. It is found that substan-
tial lateral heat flow can result from the finite extent of the slabs. It is also seen that the temperature histories of the slab materials can be affected by the finite dimensions of the 3-D slab compared with a two-dimensional case. The electromagnetic model indicates that the slab perturbs the surface fields, and, depending on the model, the effect of the slab can be clearly seen, though effects due to differences in the dip-angles of composite slabs appear to be small. Acknowledgements
This work has been supported by the Natural Sciences and Engineering Research Council of Canada. Appendix 1. The parameters and finite difference scheme for the thermal model
Computation of the thermal regime is accomplished using the conservation of energy equation: +g=V.(KVT)+H
(1)
F.W.lONESETAL.
44
where C,. is the specific heat (taken here as 1.3 x lo3 Jkgg’K-‘), p is the density, T is temperature, K is thermal conductivity, and H is heat generation per unit volume. In Cartesian coordinates in three dimensions eq. (1) is: 2
2
$+$+;
qp(z,T);=K(T)
aK aT
3K aT
2
1
x, Y, 2)
= C + ( 16n3ST3)/(
x,
=P(z,
(2)
77~”+ 120~a,e-~‘~~)
Y, z> T)ga(z)
T(x,
+ 21.4) + 996
- 1184O/(P + 8.6) + 1340
3K aT
where C = lattice conductivity at low temperature n = index of refraction = 1.7; = 2.5 Wm-‘K-l; S = Stefan-Boltzman constant = 5.67 X lo-’ Wmw2Ke4; lg = opacity at low temperature = . . 10-i . m-‘; a0 = electrical conducttvtty; E= width of energy gap for electronic conduction; k = Boltzman’s constant = 1.38 X 10-23JK-‘. In what follows, we will write K(T) = K, with the understanding that it depends on temperature. However, the temperature dependence of the lattice conductivity (C) is ignored. The density distribution follows the Bullen A density curve (see Sydora et al., 1978) and the heat generation term includes adiabatic and shear-strain heating as: H(T,
T = 7.99(P + 21.4)
Oceanic: T = 4.34(P + 8.6)
where the spatial and temperature dependences of K, p and H are indicated. The thermal conductivity is taken as in MacDonald (1959) and includes a constant lattice conductivity term and a radiative transfer term that depends on temperature: K(T)
Continental:
- 24914/(P
+;tllax+--+araz ay a~ +H(T,
The initial vertical temperature condition is from the relationships given by Mercier and Carter (1975):
Y, z)“z+Hs,,
where g = 9.81 m/ss2; a(z) = exp(3.58 0.0072~) x lo-” as in Toksoz et al. (1971) (z measured in km); u, = the vertical component of velocity of the downgoing slab; Hsh = the shearstrain heating, confined to 10 km thick layers along the tops, bottoms, sides and ends of the slabs, and can be adjusted according to the model desired.
where P = 0.31~ kbar km-‘, and the temperature over the surface of the Earth is held constant and equal to 0°C throughout the whole calculation, i.e., T(x,
y, 0, t) =0
t&O
The other boundary conditions are (where L,, L, represent the limits of the mesh in the X-, y- and z-directions, respectively). L,,
aT a2
=E.F=13.5mWm-’ r=L,
K’
aT aT ==0 ax Ix=o ax 1x=L,
aT aT =0 ay Iy=o ay Iy=L,= The solution proceeds by translating the slab, and then solving eq. (2) alternately in each of the three directions by the Douglas alternating direction implicit (ADI) method (Peaceman and Rachford, 1955). Using the notation: xi=(i-1)Ax
i=1,2,...,N,(N,-l)Ax=L,
y,=(j--l)Ay
j=1,2,...,N,(N,-l)Ay=L,
z,=(k-1)Az
k=1,2,...,N,(N,-l)Az=L,
T(x;, y,, zk, t> = T (i, j, k) for a fixed time t, and we consider the x-, y-, and z-directions as follows:
x-direction
Let T(i, j, k) be the temperature at t = to, and let T,,(i, j, k) be the temperature at t = to + At.
THERMAL
REGIME
OF A 3-D
SUBDUCTING
Then, the first estimate of TJi, T,(i, j, k) is given by: cIp[z,
T(i,
+H(i,
SLAB
a2T
1
aK aT +--+zz ay ay
1 aK
AyAz
- y ZAyAz
Ax
Ax
aK aT =R,T
j,k)+T,(i-1,
j, k)
T,(i + 1, j, k) - T,(i - 1, j, k)
i,k)= At -CJdz, f)T,(i, 2Ax
1
by
1
R(T)
(5)
where R,(T) = 2 Au R(T). If i = 2, then aT/ax 1X=o= 0 implies that T(1, j, k)= T(2, j, k) and also, T,(l, j, k)= T,(2, j, k).
Equation (5) then becomes: K-
-
AyAz
1 aK +?axAyAz+
2C,,P(ZJ) At
xT,(Z j,k)
I (4)
AX
+ 2 ~AYAz [T(i+l,
j, k) -2T(i,
j, k)
- 1, j, k)]
Av
1
1
1 aK
where R(T) =
+T(i
1 aK
+ ?~AYAz
(3)
(Ax)’
+zax
45
RESPONSE
Multiplying eq. (4) by 2AxAyAz = 2Au gives: K-
j, k, T)
j,k)-2T,(i,
1 aK
ITS ELECTROMAGNETIC
j, k), denoted by
If the derivatives in eq. (3) are replaced differences we have: T,(i+l,
AND
~,(i, j, k) - T(i, j, k)
j, k)]
a2T +K -+az2 [ a?
LITHOSPHERK
T,(3, j, k)
=R,(T) where in R,(T) we replace T(1, j, k) by T(2, j, k).
+--: E&[T(i+
+K &
1, j, k) - T(i-
[T(i,
j + 1, k) - 2T(i,
1, j, k)]
j, k)
If i =N, - 1, then aT/ax IXzL, = 0 implies RN, - 1, j, k) = TCN,, j, k) and T,(N, - 1, j, k) = T,(N,, j, k) and eq. (5) becomes: AyAz
+T(i,
K-
j - 1, k)]
&T(’
Ax
T,( N, - 2, j, k)
1 aK + z %AyAz
I, j, k + 1) - 2T( i, j, k)
+K(AZ)’
1
1 aK
- 2 zAyAz
+ T( i, j, k - l)] -2C,p(z,
aK 1 +--ay 2Ay [T(i, j + 1, k) - T(i, j - 1, k)] aK 1 +-a.2 2AzlT( +C,p(z,
i, j, k + 1) - T(i,
T)&T(i, i,k) +H(T)
j, k - 1)]
1
T);
I
T,(N,-
where in R,(T) as given j, k) by ZW, - I, j, k). Equation (4) therefore 2, NY- 1 and k = 2, N, equations with unknowns N,-
1.
1, j, k) =R,(T)
above we replace T(N,, represents for each j = - 1 a system of N, - 2 T,(i, j, k) where i = 2,
F.W.JONESETAL.
46
Note that the coefficients of this system contain the terms aK/ax, aK/ay, and aK/a.z, which may be computed by the formulas:
where R,(T)
Cd(T)
= -2
At
AvT,(i, j, k)
K(i + 1, j, k) - K(i - 1, j, k)
aK ax xi. y,, zq =
26x
DAY
ays z% Y,>Lk = aK
-
-_
-2T(i,
K(i, j, k + 1) - K(i, j, k - 1)
--
2h.z
a2 x,7Yj,ZA
y-direction
C,.p[ z,T(i,
j, k)]
K
a2T, =_+2 I ax2
-
2
I
j, k) + T(i, j - 1, k)] j+ 1, k)
-AxAz[T(i, 2 3Y j-
1, k)]
For j = 2, we replace in eq. (7) T,G, 1, k) by T2 (i, 2, k) and for j = I?, - 1, T,(i, NY, k) = T&i,
N, - 1, k). Equation (7) then represents, for each i = 2, N, - 1, and k = 2, N, - 1, a system with unknowns T,(i, j, k), j = 2, & - 1.
T2(i, i, k) - T( i, j, k) At
z-direction
a2T ax* I
ay2 +ayl
---+--
The third and final estimate for T,,(i, j, k), denoted by T&i, j, k) is now obtained. The operation in the z-direction is treated like that in the y-direction. Writing eq. (6) as:
1
a2T
a2T2
K
j+ 1, k)
ale
i
-T(i,
A second estimate of T,,(i, j, kl is now obtained and denoted by T,(i, j, k). In eq, (31, the terms that contain the derivatives of T with respect to y can be separated into two parts and we obtain:
f-
F[T(i,
-K
K(i, j + 1, k) - K(i, j - 1, k)
aK
I aK aT
C,.p[z,T(i,j, k)]
T,(i, j, k) - T(i, j, k) At
ay ay I aK aT
a*T
+Ka,2+
zz
+H(i, i, k, T)
(6)
Subtracting eq. (6) from eq. (3), replacing the derivatives with respect to p by differences, and rearranging gives: AxAz 1 aK - --AxAz KAY
2 ay
I
T2(i, j-
1 +y”
1, k)
i
aT,2 -+g7 az2
a2T I
I
=&(T)
17)
+H(i,j,k,T)
and subtracting eq. (8) from eq. (61, replacing the derivatives with respect to z by differences, and
THERMAL
REGIME
rearranging AxAyAz): AxAy
K-
AZ
OF A 3-D
terms
SUBDUCTING
gives
1 aK
- z ,zAxAy
I
LITHOSPHERIC
(again
SLAB
with
Au =
T3(i, j, k - 1)
j,k) 1T3(i,
1
1 aK
+ ?GAXAY
T3(i,
where At
x [T(i,
AuT2(i, j, k) - KF
j, k + 1) - 2T(i,
j, k)
+ T( i, j, k - l)]
1 aK --- 2 az AxAy[T(i, -T(i,
j, k-
j, k + 1)
l)]
Because T(x, y, 0, t ) = 0 for any t, we have T3 (i, j, 1) = 0.
From the boundary have:
condition
at z = L,,
=R,(T) andfork=N,-1: 1 aK
- 2 %AxAy
1
T3(i, j, N, - 2)
1
P(T)
+ 2C,;?Av
XT3(i, i, N, - 1) =R3(T)
-
T(i,
1, k) = T(i, 2, k),
T(i,
N,,, k) = T(i,
N,, - 1, k)
T(i,
j, N,) = T(i,
j, N, - 1) + 7
AxAy
KF
Az.F
The solution therefore proceeds by simulating the slab movement by translating the slab temperatures through the mesh and subsequently applying the above AD1 method to solve eq. (1). The numerical mesh must be such as to accommodate the required rate and angle of subduction. Movement of the slab can cease, with subsequent warming if desired. References
Equation (9) for k = 2 becomes:
AxAy
j, k) = T(N, - 1, j, k)
we
F T3(i, j, N,) = T3( i, j, N, - 1) + AZ. z
K - AZ
T(l, j, k) = T(2, j, k),
AxAy
C,P(Z~)
41
RESPONSE
Equation (9) for each i = 2, N, - 1 and j = 2, N,, - 1 represents a system of equations with unknowns T&i, j, k), k = 2, N, - 1. After completing the computations in the three directions, the temperatures at the boundaries are computed according to the boundary conditions:
T(N,,
(9)
= -2
ITS ELECTROMAGNETIC
j, k + 1)
=RdT)
R,(T)
AND
F
1 aK
+ ?%AxAy
%Az. I
Bodri, L. and Bodri, B., 1978. Numerical investigation of tectonic flow in island-arc areas. Tectonophysics, 50: 163115. De Jonge, M.R. and Wortel, M.J.R., 1990. The thermal structure of the Mediterranean upper mantle: a forward modelling approach. Terra Nova, 2: 609-616. Duba, A. and Lilley, F.E.M., 1972. Effect of an ocean ridge model on geomagnetic variations. J. Geophys. Res., 72: 7100-7105. Hsui, A.T., Marsh, B.D. and Toksoz, M.N., 1983. On melting of the subducting oceanic crust: effects of subduction induced mantle flow. Tectonophysics, 99: 207-220. Isacks, B., Oliver, J., and Sykes, L.R., 1968. Seismology and the new global tectonics. J. Geophys. Res., 73: 5855-5899. Jones, F.W. and Pascoe, L.J., 1972. The perturbation of alternating geomagnetic fields by three-dimensional conductivity inhomogeneities. Geophys. J.R. Astron. Sot., 27: 479-485. Jones, F.W. and Vozoff, K., 1978. The calculation of magnetotelluric quantities for three-dimensional conductivity inhomogeneities. Geophysics, 43: 1167-l 175. Katsumata, M. and Sykes, L.R., 1969. Seismicity and tectonics of the western Pacific, Izu-Mariana, Caroline and Ryukya-Taiwan regions. J. Geophys. Res., 74: 5923-5948.
48 Lines, L.R. and Jones, F.W., 1973a. The perturbation of alternating geomagnetic fields by three-dimensional island structures. Geophys. J.R. Astron. Sot., 32: 133-154. Lines, L.R. and Jones, F.W., 1973b. The perturbation of alternating geomagnetic fields by an island near a coastline. Can. J. Earth Sci., IO: 510-518. Lister, C.R.B., 1971. Tectonic movement in the Chile trench. Science, 173: 719-722. Lliboutry, L., 1969. Sea-floor spreading, continental drift and lithosphere sinking with an asthenosphere at melting point. J. Geophys. Res., 74: 6525-6540. Ludwig, W.J., Ewing, J.E., Ewing, M. Murauchi, S., Den, N., Asano, S., Hotta, H., Hayakawa, M., Asanuma, T., Ichikawa, K. and Noguchi, I., 1966. Sediments and structure of the Japan trench. J. Geophys. Res., 71: 2121-2137. MacDonald, G.J.F., 1959. Calculations on the thermal history of the Earth. J. Geophys. Res., 64: 1967-2~. Malahof, A., 1970. Some possible mechanisms for gravity and thrust faults under ocean trenches. J. Geophys. Res., 75: 1992-2001. Mercier, J.C. and Carter N.L., 1975. Pyroxene geotherms. J. Geophys. Res., 80: 3349-3362. Minear, J.W. and Toksoz, M.N., 1970. Thermal regime of a downgoing slab and new global tectonics. J. Geophys. Res., 75: 1397-1419.
F.W. JONES
ET AL.
Peaceman, D.W. and Rachford, Jr., H.H., 1955. The numerical solution of parabolic and elliptic differential equations. J. Sot. Ind. Appl. Math., 3: 28-41. Rikitake, T., 1966. Electromagnetism and the Earth’s Interior, Elsevier, Amsterdam, 308 pp. Stauder, W., 1968a. Mechanism of the Rat Island earthquake sequence of February 4, 1965, with relation to island arcs and sea floor spreading. J. Geophys. Res., 73: 3847-3858. Stauder, W., 1968b. Tensional character of earthquake foci beneath the Aleutian trench with relation to sea floor spreading. J. Geophys. Res., 73: 7693-7701. Sydora, L.J., Jones, F.W. and Lambert, R. St.J., 1978. The thermal regime of the descending lithosphere: The effect of varying angle and rate of subduction. Can. J. Earth Sci., 154: 626-641. Toksoz, M.N. and Hsui, A.T., 1978. Numerical studies of back-arc convection and the formation of marginal basins. Tectonophysics, 50: 177-196. Toksoz, M.N., Minear, J.W. and Julian, B.R., 1971. Temperature field and geophysical effects of a downgoing slab. J. Geophys. Res., 76: 1113-1138. Wyllie, P., 1971. The Dynamic Earth. Wiley, New York, NY, 416 pp.