Geophys. J. Int. (2011) 186, 10951103 doi: 10.1111/j.1365-246X.2011.05113.

x
G
J
I
M
i
n
e
r
a
l
p
h
y
s
i
c
s
,
r
h
e
o
l
o
g
y
,
h
e
a
t

o
w
a
n
d
v
o
l
c
a
n
o
l
o
g
y
Subcritical dyke propagation in a host rock with
temperature-dependent viscoelastic properties
Zuan Chen
1
and Z.-H. Jin
2
1
Key Laboratory of the Study of Earths Deep Interior, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China.
E-mail: zachen@mail.igcas.ac.cn
2
Department of Mechanical Engineering, University of Maine, Orono, ME 04469, USA
Accepted 2011 June 12. Received 2011 June 12; in original form 2009 November 19
SUMMARY
In this paper, we examine the effects of temperature-dependent viscoelastic properties of
the host rock on the subcritical growth of a dyke from a magma chamber. A theoretical
relationship between the velocity of subcritical dyke growth and dyke length is established
using a perturbation solution of stress intensity factor at the dyke tip and a viscoelastic crack
growth theory in which the temperature-dependent creep properties are taken into account.
The temperature eld around the dyke is calculated using an analytic solution. The numerical
results for a dyke subcritically propagating from a magma chamber indicate that while the
general dyke growth characteristics are similar to those with constant creep properties, the
subcritical dyke growth velocity is increased by an order of magnitude by considering the
temperature dependence of the creep properties. Hence, the subcritical growth duration before
the dyke reaches the unstable growth state is signicantly shortened.
Key words Mechanics, theory, and modelling; Magma migration and fragmentation; Volcano
monitoring.
1 I NTRODUCTI ON
Magma transport through dyke propagation signicantly inuences
the geological evolution of the Earths crust and volcanism on the
Earths surface. Signicant progress has been made in understand-
ing the propagation of dykes through the upper mantle and crust
(e.g. Spence & Turcotte 1985; Spence et al. 1987; Pollard & Segall
1987; Lister &Kerr 1991; Clemens &Mawer 1992; Rubin, 1995a,b,
1998; Meriaux & Jaupart 1998; Bonafede & Rivalta 1999; Dahm
2000; Ito & Martel 2002; Rivalta & Dahm 2006; Chen et al. 2007;
Jin & Johnson 2008).
Most work on dyke propagation has focused on the critical prop-
agation of dykes, that is, the stress intensity factor (SIF) at the dyke
tip is equal to or greater than the fracture toughness of the host
rock (e.g. Weertman 1971; Spence & Sharp 1985; Spence et al.
1987; Lister 1990, 1991, 1994a,b; Rubin 1995a,b, 1998; Bonafede
&Rivalta 1999; Meriaux &Jaupart 1998; Bolchover &Lister 1999;
Dahm 2000; Menand & Tait 2002; Kuhn & Dahm 2004; Roper &
Lister 2005; Rivalta & Dahm 2006; Chen et al. 2007; Taisne &
Jaupart 2009; Traversa et al. 2010; Maccaferri et al. 2010). It is ex-
pected that small cracks form at the wall of a magma chamber. The
SIFs at the tips of those small magma-lled cracks, or dykes, may
not reach the fracture toughness and therefore will not propagate
critically. However, these small dykes may propagate subcritically
due to stress corrosion, viscous damage and other mechanisms.
Anderson & Grew (1977) presented a stress corrosion cracking
model as applied to magma fracture initiation. In their model, sub-
critical growth is due to the weakening of strained atomic bonds
at the crack tip by the chemical action of the environment agents.
Atkinson (1984) reviewed the progress on subcritical crack growth
in geological materials. He discussed common chemical mecha-
nisms responsible for the subcritical growth and presented detailed
experimental results for various rocks including granite, marble
and limestone. Rubin (1998) examined subcritical crack growth in
partially molten rocks from a surface energy point of view. Chen
& Jin (2006) investigated subcritical dyke growth from a crustal
magma chamber. Using a viscoelastic energy dissipation approach
without considering temperature-dependent viscoelastic properties,
they obtained the relationship between the subcritical propagation
velocity and dyke length, and found that dykes grow at velocities
on the order of 10
8
10
6
m s
1
.
Temperature dependence of the viscoelastic properties of the
host rock may play an important role in energy balance during
slow, subcritical dyke propagation from a magma chamber. In this
paper, we investigate subcritical dyke propagation from a magma
chamber using an energy balance approach with viscoelastic en-
ergy dissipation. The temperature dependence of the viscoelastic
properties and the inuence of heat conduction from magma to the
host rock around the dyke tip are considered. The relationship be-
tween the subcritical growth velocity and dyke length is established.
C
2011 The Authors 1095
Geophysical Journal International
C
2011 RAS
Geophysical Journal International
1096 Z. Chen and Z.-H. Jin
Numerical results are given to illustrate the effects of magma tem-
perature, magma chamber overpressure and fracture toughness on
the subcritical dyke growth velocity.
2 THEORETI CAL MODELS
2.1 Magma pressure along the dyke surfaces
Consider a dyke that subcritically propagates vertically from a
magma chamber into the host rock. The dyke length is usually
much smaller than the size of magma chamber during the subcriti-
cal propagation stage. The dyke thus may be considered as a magma
lled, plane strain edge crack in a half plane, as shown in Fig. 1,
where a = a(t) denotes the crack length and t is time. This 2-D
dyke propagation model has been adopted in a number of studies,
for example, Rubin (1995a,b), Roper & Lister (2005) and Chen &
Jin (2006).
Dyke propagation is driven by the net pressure on its surfaces
due to magma ow, lithostatic stress and tectonic stress. The uid
pressure p may be expressed as
p = p
e
(
s
gZ +), (1)
where p
e
= p
e
(Z) is the pressure due to elastic deformation of the
host rock, Z the vertical coordinate, g the gravitational acceleration,

s
the density of the host rock and the tectonic stress (Ru-
bin1995a,b; Rope & Lister 2005). Rope & Lister (2005) assumed
that could vary linearly in Z-direction and thus can be absorbed
into the lithostatic stress using an effective rock density
eff
.

eff
gZ =
s
gZ +. (2)
The elastic pressure along the dyke surfaces has been derived
using the lubrication theory of uid mechanics as follows (Chen
et al. 2007):
p
e
= P +g(z +a) 12V
_
z
a

2
dz, (3)
where P is the overpressure in the magma chamber, =
eff

m
,
m
is the magma density, is the magma viscosity, V is the dyke
propagation velocity, is the separation of the two dyke surfaces
and z = Z a. When z = a or Z = 0, p
e
= P.
Figure 1. A dyke in a semi-innite medium (a: dyke length, : dyke thick-
ness, p: dyke surface pressure).
2.2 Temperature eld and viscoelastic constitutive
equations with temperature-dependent creep properties
Subcritical dyke growth in the crust takes place generally near a
magma chamber where the material properties of the host rock are
greatly inuenced by the elevated temperatures. The temperature
distribution around a dyke satises the following heat conduction
equation:

2
T
Z
2
+

2
T
X
2
= 0, (4)
where T is the temperature. Here the transient effect of the temper-
ature eld is neglected as the subcritical dyke growth velocities are
very low (Chen et al. 2006). The temperatures along the magma
chamber (Z = 0) and dyke (X = 0, 0 Z a) are assumed as the
melting temperature of the magma, T
m
. By solving eq. (4) we can
get the temperature eld as follows:
T(X, Z) = T
m

T
m
T
s
h
Z
1
2
_
a
0
ln

_
X
2
+(Z

+ Z)
2
X
2
+(Z

Z)
2
_
f (Z

) dZ

,
(5)
where h is the depth of dyke base, T
s
is the temperature on the
Earths surface and f (Z

) is determined from the integral equation

(A7) in the Appendix in which the detailed derivation of eq. (5) is
given.
The relationship between creep strain rate and differential stress
may be modelled by the following power law for crustal rock under
high-temperature creep conditions (Kirby et al. 1987):
.

c
= A(
1

3
)
n
exp
_

Q
RT
_
, (6)
where
c
is creep strain,
i
(i = 1, 3) are stresses, T is the absolute
temperature, A, n and Q are constants related to rock properties and
R is the universal gas constant. Eq. (6) is generally a non-linear
relation between the strain rate and the stress. For simplicity, we use
a special case of n = 1 as follows:

c
= A

(
1

3
) exp
_

Q
RT
_
, (7)
where A

is a material constant. Eq. (7) will be used in Section 3 for

calculating the viscous dissipative energy.
2.3 SIF calculation
For the sake of simplicity, we neglect the effect of the temperature
dependence of elastic modulus in calculating the SIF at the dyke
tip. Eq. (3) indicates that the boundary condition on the dyke (or
crack) surfaces is non-linear and related to the crack opening. Using
a combined perturbation/integral equation approach, Chen et al.
(2007) obtained the expression of the SIF at the dyke tip as follows:
K
I
=
ga

a
2
[

0
(1) +

1
(1) +
2

2
(1) +
3

3
(1)], (8)
where

i
(r)(i = 0, 1, 2, 3) are dimensionless functions given by

i
(r) =

1 r
i
(r), i = 0, 1, 2, 3, (9)
and the dimensionless functions
i
(r)(i = 0, 1, 2, 3) satisfy the
C
2011 The Authors, GJI, 186, 10951103
Geophysical Journal International
C
2011 RAS
Dyke propagation with viscoelastic properties 1097
following integral equations:
_
1
1
_
1
s r
+ K(r, s)
_

0
(s) ds =
2
ga
P (1 +r),
_
1
1
_
1
s r
+ K(r, s)
_

1
(s) ds =
_
E
(1
2
)ga
_
2
_
r
1
_
H

0
_
2
ds,
_
1
1
_
1
s r
+ K(r, s)
_

2
(s) ds =
_
E
(1
2
)ga
_
2

_
r
1
_
2

0
_
H

0
_
2
_
ds,
_
1
1
_
1
s r
+ K(r, s)
_

3
(s) ds =
_
E
(1
2
)ga
_
2

_
r
1
_
_
H

0
_
2
_
3

2
1

2
0
2

0
_
_
ds,
(10)
respectively, and is a non-dimensional perturbation parameter
given by
=
12V
H
2
1
g
. (11)
In eqs (9)(11), r = 2Z/a 1, Z is the vertical coordinate as
shown in Fig. 1, K(r, s) is a known kernel, H is a length parameter
having an order of the dyke base thickness, and

i
(r) (i =0, 1, 2, 3)
are related to crack opening given by

i
(r) = a
_
1
r

i
(s) ds, i = 0, 1, 2, 3. (12)
It is noted that P is assumed as a constant as in most dyke prop-
agation studies while the buoyancy pressure ga varies with dyke
length. The integral equations in eq. (10) are solved numerically.
2.4 Analysis of subcritical crack propagation in a
viscoelastic solid with temperature-dependent creep
properties
Under isothermal conditions, the appropriate statement of the global
conservation of energy for crack extension in viscoelastic media is
(Christensen 1982)
dU
dt
+
dD
p
dt
+
dS
e
dt
= 0, (13)
where U is the elastic strain energy per unit thickness, D
p
is the
viscous dissipation energy, S
e
is the surface energy and t is time.
For the growth of an edge crack of length a, eq. (13) may also be
written as
dU
da
+
dD
p
da
+
dS
e
da
= 0, (14)
where da = Vdt with V being the growth velocity.
To study subcritical dyke propagation, we assume that an initial
short dyke of length a
0
is rst produced around a magma chamber
(Sleep 1988; Fowler 1990). The dyke then subcritically grows in a
step-by-step manner, that is, the dyke grows to a length a through
a time period Tg(a
0
) and then stops. After another time period
Tg(a), the dyke continues to grow to a new length and stops again.
Furthermore, there exists a critical length a
c
at which the SIF at the
dyke tip reaches the critical SIF. Tg(a) is equal to zero when a is
larger then a
c
, which means the termination of subcritical growth
stage and the dyke will growunstably until failure occurs. In general,
Tg(a) is the time period from the stop of the last dyke extension up
to the start of the next dyke extension. Tg(a) may be understood as
the time for unit subcritical dyke extension. The analysis process
mentioned earlier is a numerical strategy to solve our problem.
For the numerical process of subcritical crack extension described
earlier, the energy terms dU/da, dD
p
/da and dS
e
/da in eq. (14) have
been obtained in Chen (2003), Chen & Bai (2006) and Chen &
Jin (2006) for a mode I crack/dyke for plane strain. We summarize
these descriptions as follows:
For mode I crack extension problems for plane strain, the strain
energy change rate relates to the SIF by (Knott 1973)
dU
da
=
2(1
2
)
E
K
2
I
, (15)
Figure 2. The temperature distribution around a magma chamber and dyke.
C
2011 The Authors, GJI, 186, 10951103
Geophysical Journal International
C
2011 RAS
1098 Z. Chen and Z.-H. Jin
where E is Youngs modulus, is Poisson ratio and K
I
is the SIF.
Here, we neglect the dependence of Youngs modulus on the tem-
perature. The stresses around the crack tip do not vary with time
during a small crack growth step in accordance with the correspon-
dence principle (Christensen 1982). The crack tip stress eld for a
mode I crack is

i j
=
K
I

2r
f
i j
(), (16)
where r and are the polar coordinates centred at the crack tip at
the end of the step,
ij
is stress tensor and f
ij
() are the angular
distributions of stresses.
The viscoelastic constitutive equations for creep may be written
as follows:
e
i j
= s
i j
G(t ), e
v
= pK(t ),
e
i j
=
i j
e
v
, e
v
=
i i
/3,
s
i j
=
i j
p
v
, p
v
=
i i
/3,
(17)
where G(t) and K(t) are the deviatoric and volumetric compliances,
respectively, s
ij
is the deviatoric stress tensor, p
v
is the hydrostatic
stress, e
ij
is the deviatoric strain tensor, e
v
is the volumetric strain
and
ij
is the strain tensor.
In a crack extension/stop step, the dissipative energy density Q
e
at a point of the medium may be calculated as
Q
e
=
_
Tg
0

i j
.

i j
d , (18)
where Tg is the time period from the stop of the last dyke extension
up to the start of the next dyke extension and a dot over
ij
denotes
strain rate.
It is assumed that the dissipation of energy is due to the stress
singularity at the crack tip. Combining eqs. (16) and (18), the total
dissipative energy dD
p
per unit thickness in the step can be calcu-
lated by integrating Q
e
over a crack tip area surrounded by a small
circle around the crack tip as follows:
dDp =
_
da
0
_

Q
e
rddr = K
2
I
[B(Tg) B(0)] da, (19)
where B(Tg) is dened by
B(Tg) =
_
25
24

4
3

_
G(Tg) +
4
3
(1 +)K(Tg), (20)
In eq. (20), G(Tg) and K(Tg) are the deviatoric and volumetric
compliances of the viscoelastic solid.
The surface energy rate is given by
dSe
da
= 2, (21)
where represents the surface energy per unit area and is related
to the fracture toughness K
Ic
by (Jaeger and Cook, 1976)
=
(1
2
)
E
K
2
I c
. (22)
Substituting eqs. (15), (19) and (21) into eq. (14), we have
B(Tg) =
4
K
2
I
_
_
K
2
I c
K
2
I
_
(1
2
)
E
_
+ B(0). (23)
Figure 3. Dyke subcritical growth velocity versus the magma temperature for various dyke lengths (P = 3.0 Mpa, K
IC
= 52 MPa m
1/2
).
C
2011 The Authors, GJI, 186, 10951103
Geophysical Journal International
C
2011 RAS
Dyke propagation with viscoelastic properties 1099
The above eq. (23) must be satised by Tg(a).
Next we evaluate the function B(Tg) for the following
temperature-dependent creep functions (Miannay 2001):
G(Tg) = A

Tg exp
_

Q
RT
_
,
K(Tg) = 0.
(24)
Eq. (24) may be obtained from eq. (7) for creep process with the
rock constants determined by creep experiments. From the numer-
ical result of the temperature distribution around a dyke (as shown
in Fig. 2), the temperature in eq. (24) may be taken as the magma
temperature. Substituting eqs (20) and (24) into (23) and solving
for Tg = Tg(a), we have
Tg(a) =
288 (1
2
)
A

E (75 96)
_
_
K
I C
K
I
_
2
1
_
exp
_
Q
RT
_
. (25)
In accordance with the denition of Tg(a) (time period for unit
subcritical dyke extension), the relation between Tg(a) and the sub-
critical dyke growth velocity, V, can be established as
V =
1
Tg(a)
. (26)
3 RESULTS AND DI SCUSSI ON
In our calculations, we use the following typical properties for the
host rock and magma (Rubin 1995a,b; Kirby et al. 1987): E = 50
GPa, = 0.25, = 300 kg m
3
, = 50 Pa-s, g = 9.8 m s
2
.
Thermally activated Newtonian creep (eq. 7) of crustal rocks corre-
sponds to a viscosity decreasing with temperature. The temperature-
dependent viscosity of rock
rock
may be expressed as follows:

rock
=
1
A

exp
_
Q
RT
_
, (27)
Using A

=0.0001 MPa
1
s
1
and Q=137 kJ mol
1
(Turcotte &
Schubert 2002) yields
rock
=4.049 10
14
, 4.175 10
15
, 1.569
10
18
and 4.272 10
20
Pa-s at temperatures of 1280

C, 1000

C,
600

C and 400

C, respectively. We perform a parametric study
on the effects of physical parameters (fracture toughness, magma
temperature and magma chamber pressure) on the subcritical dyke
growth velocity. The SIFs are calculated using eq. (8).
In this work, we only consider steady-state temperature distri-
bution by thermal conduction process based on the fact that dyke
subcritical growth is very slow. Fig. 2 shows the isothermal curves
(constant temperature) calculated from eq. (5). The magma tem-
perature, dyke length, dyke base depth and surface temperature
are assumed as T
m
= 1020

C, a = 500 m, h = 20000 m and
T
s
= 0, respectively. It can be seen that the temperature near the
dyke is close to the magma temperature and decreases with increas-
ing distance from the dyke. The numerical results suggest that the
magma temperature may be used in eq. (24) for determining the
temperature-dependent creep properties used for crack tip energy
balance calculations.
The dyke subcritical growth velocity V versus magma temper-
ature behaviour is shown in Figs 35. The effect of initial dyke
length on the subcritical growth behaviour is examined in Fig. 3
where the magma chamber overpressure is P = 3.0 MPa and
Figure 4. Dyke subcritical growth velocity versus the magma temperature for various dyke lengths (P = 3.0 Mpa, K
IC
= 100 MPa m
1/2
).
C
2011 The Authors, GJI, 186, 10951103
Geophysical Journal International
C
2011 RAS
1100 Z. Chen and Z.-H. Jin
fracture toughness K
IC
= 52 MPa m
1/2
. It is seen from Fig. 3 that
at a given dyke length, V increases with an increase in T
m
. The
subcritical velocity is on the order of 10
7
10
4
m s
1
, much
smaller than the unstable propagation velocities in the range of
0.0110 m s
1
(Rubin 1995a,b). The subcritical growth velocity
increases rapidly when the dyke is approaching the critical length
(71 m) at which the SIF reaches the fracture toughness. The results
indicate signicant effects of temperature-dependent creep proper-
ties of the host rock on the subcritical growth velocity. For example,
at a dyke length of 70 m, V is 1.0 10
4
m s
1
at T
m
= 1050

C.
V increases to 5 10
4
m s
1
when T
m
is 1250

C. Clearly, high
temperature-induced thermal dissipation in the host rock results in
earlier unstable dyke propagation than that without considering the
effects of temperature-dependent creep properties of the host rock.
Fig. 4 shows the subcritical growth velocity versus magma temper-
ature for various dyke length at an increased fracture toughness of
K
IC
=100 MPa m
1/2
. This fracture toughness value may be near the
high end of the range for basaltic oceanic crust as the toughness of
crustal rock generally increases with increasing conning pressure
(Jin & Johnson 2008). The magma chamber overpressure is still
P = 3.0 MPa. Similar subcritical dyke growth behaviour can be
observed. However, the dyke growth velocity is signicantly lower
compared with that for the case of K
IC
= 52 MPa m
1/2
shown in
Fig. 3 under the same magma temperature and dyke length. We note
that the critical dyke length is now 220 m and longer than that for
K
IC
= 52 MPa m
1/2
.
The effect of magma chamber overpressure P on the subcriti-
cal dyke growth behaviour is examined in Fig. 5. The initial dyke
length is now xed at 85 m and the fracture toughness is K
IC
=
100 MPa m
1/2
. It is seen that the subcritical growth velocity in-
creases with increasing magma chamber pressure for a given magma
temperature. Hence, higher P will result in early unstable dyke
propagation.
Fig. 6 shows the dyke growth duration versus dyke length under
the conditions of P = 3 MPa, K
IC
= 52 MPa m
1/2
and T
m
=
1000

C, 1100

C and 1280

C. When the temperature of magma
is 1000

C, dyke subcritical growth is very slow and it needs about
1.48 yr for the dyke to reach 70 m, which is almost equal to the
critical dyke length 71 m. When the temperature of magma is 1100

C, dyke subcritical growth becomes faster and it needs about 0.58

yr for the dyke to reach 70 m. When the temperature of magma is
1280

C, dyke subcritical growth speeds up and it only needs about
0.14 yr for the dyke to reach 70 m. Dyke growth becomes unstable
when the dyke approaches 71 m.
According to Turcott & Schubert (2002), the solidication time
t
s
of a dyke may be calculated as follows:
t
s
=
b
2
4
2
2
, (28)
where b is the dyke width, is the thermal diffusivity and
2
is a parameter depending on the difference between the magma
temperature and the host rock temperature. Taking a dyke width
b = 1 m (which is approximately equal to the opening at the dyke
base in our calculations), thermal diffusivity = 0.5 mm
2
s
1
and

2
= 0.2 (which corresponds a temperature difference T
m
T
0
=
100

C, see Fig. 2), we get t
s
=144.7 d, or 0.406 yr. The calculation
above is based on a stationary dyke and hence is conservative. The
actual solidication time may be longer for a growing dyke as new
Figure 5. Dyke subcritical growth velocity versus the magma temperature for various magma chamber pressures (a = 85 m, K
IC
= 100 MPa m
1/2
).
C
2011 The Authors, GJI, 186, 10951103
Geophysical Journal International
C
2011 RAS
Dyke propagation with viscoelastic properties 1101
Figure 6. The relationship between the dyke growth duration and its length (T
m
= 1000

C, 1100

C and 1260

C, P = 3.0 Mpa, K
IC
= 52 MPa m
1/2
).
magmas ow from the magma chamber into the dyke during prop-
agation. Hence, the dyke may solidify during its subcritical growth
when the magma temperatures are 1000

C and 1100

C. When
the magma temperature is 1280

C, the dyke can reach the critical
state of unstable propagation before complete solidication occurs.
In our analysis, the elastic modulus E is assumed to be constant
and does not depend on the temperature. Based on the experimental
results of Liu et al. (2000), the elastic modulus of granite reduces
to about 77 per cent of its room temperature value at a temperature
of 1280

C. Assuming a temperature-independent elastic modu-
lus will likely cause a relative error on the order of 20 per cent
in the calculated subcritical growth velocity. Consideration of the
temperature-dependent creep properties, however, will increase
the subcritical velocity by an order of magnitude as indicated in
Figs 35 as compared with our previous results with temperature-
independent creep properties. Hence, temperature dependence of
creep properties is more dominant and the effect of temperature
dependence of elastic modulus on the subcritical dyke propagation
may be neglected.
4 . CONCLUSI ON
In this paper, we obtain the temperature eld around a dyke subcrit-
ically growing froma magma chamber and investigate the effects of
temperature-dependent creep properties of the host rock on the dyke
subcritical growth behaviour. An energy balance approach is used
which considers energy dissipations due to viscoelasticity including
temperature effect. The numerical results indicate that the energy
dissipation due to higher temperature of the host rock at the dyke
tip results in higher subcritical propagation velocity compared with
that without considering the effects of temperature-dependent creep
properties of the host rock. The subcritical propagation velocity in-
creases with an increase in the magma temperature. We also nd
that the fracture toughness K
Ic
and the magma chamber overpres-
sure P have profound effects on the subcritical dyke propagation
velocity. Subcritical growth of dyke is a generally slowprocess with
the velocity in the range of 10
7
10
4
m s
1
before the dyke length
reaches its critical length. A subcritically growing dyke may or may
not reach the critical state of unstable propagation depending on
the magma solidication rate and the temperature in the host rock
around the magma chamber. Inclusion of the effects of temperature-
dependent creep properties of the host rock will increase the chance
for the dyke to reach the unstable state.
ACKNOWLEDGMENTS
This work is supported by the National Natural Science Foundation
of China (40874046). The authors would also like to thank the editor
and two reviewers for helpful discussions and comments.
REFERENCES
Anderson, O.L. & Grew, P.C., 1977. Stress corrosion theory of crack prop-
agation with applications to geophysics, Rev. Geophys., 15, 77103.
Atkinson, B.K., 1984. Subcritical crack growth in geological materials, J.
geophys. Res., 89(B6), 40774114.
Bolchover, P. &Lister, J.R., 1999. The effect of solidication on uid-driven
fracture, Proc. R. Soc. London A, 455, 23892409.
Bonafede, M & Rivalta, E., 1999. On tensile cracks close to and across the
interface between two welded elastic half-space, Geophys. J. Int., 138,
410434.
C
2011 The Authors, GJI, 186, 10951103
Geophysical Journal International
C
2011 RAS
1102 Z. Chen and Z.-H. Jin
Chen, Z. 2003. Analysis of a microcrack model and constitutive equa-
tion for time-dependent dilatancy of rock, Geophys. J. Int., 155, 601
608.
Chen, Z. & Bai, W., 2006. Fault creep growth model and its relationship
with occurrence of earthquakes, Geophys. J. Int., 165, 272278.
Chen, Z. & Jin, Z.-H., 2006. Magma-driven subcritical crack growth and
implication for dike initiation from a magma chamber, Geophys. Res.
Lett., 33, L19307, doi:10.1029/2006GL026979.
Chen, Z., Jin, Z.-H. & Johnson, S.E., 2007. A perturbation solution for dike
propagation in an elastic medium with graded density, Geophys. J. Int.,
169, 348356.
Christensen, R. M., 1982. Theory of Viscoelasticity: An Introduction, 2nd
edn, Academic, New York.
Clemens, J.D. & Mawer, C.K., 1992. Granitic magma transport by fracture
propagation, Tectonophysics, 204, 339360.
Dahm, T., 2000. Numerical simulations of the propagation and arrest of
uid-lled fractures in the Earth, Geophys. J. Int., 141, 623638.
Fowler, A.C., 1990. A compaction model for melt transport in the Earths
asthenosphere, Part II, in Magma Transport and Storage, pp.1632, ed.
Ryan, M.P., Wiley, Chichester.
Ito, G. & Martel, S.J., 2002. Focusing of magma in the upper
mantle through dike interaction. J. geophys. Res., 107(B10), 2223,
doi:10.1029/2001JB00251.
Jaeger, T.C. & Cook, N.G.W., 1976. Fundamentals of Rock Mechanics, 2nd
edn, Chapman & Hall, London.
Jin, Z.-H. & Johnson, S.E., 2008. Magma-driven multiple dike propagation
and fracture toughness of crustal rocks, J. geophys. Res., 113, B03206,
doi:10.1029/2006JB004761.
Kirby, S.H., Kronenberg, A.K., 1987. Rheology of lithosphere: selected
topics, Rev. Geophys., 25, 12191244.
Knott, J.F., 1973. Fundamentals of Fracture Mechanics, Butterworths, Lon-
don.
Kuhn, D. & Dahm, T., 2004. Simulation of magma ascent by dykes in the
mantle beneath mid-ocean ridges. J. Geodyn., 38, 147159.
Lister, J.R., 1990. Buoyancy-driven uid fracture: the effects of mate-
rial toughness and of low-viscosity precursors, J. Fluid Mech., 210,
263280.
Lister, J.R., 1991. Steady solutions for feeder dikes in a density-stratied
lithosphere, Earth planet. Sci. Lett., 107, 233242.
Lister, J.R. 1994a. The solidication of buoyancy-driven ow in a exible-
walled channel, part 1 constant-volume release, J. Fluid Mech., 272,
2144.
Lister, J.R. 1994b. The solidication of buoyancy-driven ow in a exible-
walled channel, part 2 continual release, J. Fluid Mech., 272, 4565.
Liu, Q. & Xu, X. 2000. Damage analysis of brittle rock at high temperature,
Chinese J. Rock Mech. Eng., 19(4), 408411.
Maccaferri, F., Bonafede, M. & Rivalta, E., 2010. A numerical model of
dyke propagation in layered elastic media, Geophys. J. Int., 180, 1107
1123.
Menand, T. & Tait, S.R., 2002. The propagation of a buoyant liquid-lled
ssure from a source under constant pressure: an experimental approach,
J. geophys. Res., 107(B11), 2306, doi:10.1029/2001JB000589.
Meriaux, C. &Jaupart, C., 1998. Crack propagation through an elastic plate,
J. geophys. Res., 103(B8), 18 29518 314.
Miannay, D.P., 2001. Time-Dependent Fracture Mechanics, Springer-Verlag,
New York.
Pollard, D.D. and Segall, P., 1987. Theoretical displacements and stresses
near fracture in rock: with applications to faults, joints, veins, dikes,
and solution surfaces, in Fracture Mechanics of Rock, pp. 277349, ed.
Atkinson, B.K., Academic, San Diego, CA.
Rivalta, E. & Dahm T., 2006. Acceleration of buoyancy-driven fractures
and magmatic dikes beneath the free surface, Geophys. J. Int., 166,
14241439.
Roper, S.M. & Lister, J.R., 2005. Buoyancy-driven crack propagation from
an over-pressure source, J. Fluid Mech., 536, 7998.
Rubin, A.M., 1995a. Propagation of magma-lled crack, Ann. Rev. Earth
planet. Sci., 23, 287336.
Rubin, A.M., 1995b. Getting granite dikes out of the source region, J.
geophys. Res., 100(B4), 59115929.
Rubin, A.M., 1998. Dike ascent in partially molten rock, J. geophys. Res.
103(B9), 20 90120 919.
Sleep, N.H., 1988. Tapping of melt by veins and dikes, J. geophys. Res., 93,
10 25510 272.
Spence, D.A. & Turcotte, D.L., 1985. Magma-driven propagation of cracks,
J. geophys. Res., 90, 575580.
Spence, D.A., Sharp, P. & Turcotte, D.L., 1987. Buoyancy-driven crack
propagation: a mechanism for magma migration, J. Fluid Mech., 174,
135153.
Taisne, B. & Jaupart, C., 2009. Dike propagation through layered rocks, J.
geophys. Res., 114, B09203, doi:10.1029/2008JB006228.
Traversa, P., Pinel, V. & Grasso, J.R., 2010. A constant inux model for dike
propagation: implications for magma reservoir dynamics, J. Geophys.
Res., 115, B01201, doi:10.1029/2009JB006559.
Turcotte, D. & Schubert, G., 2002. Geodynamics, Cambridge University
Press, New York, NY.
Weertman, J. 1971. Theory of water-lled crevasses in glaciers applied to
vertical magma transport beneath oceanic ridges, J. geophys. Res., 76,
11711183.
APPENDI X: DETAI LED DERI VATI ON
OF TEMPERATURE FI ELD
The boundary conditions of the steady-state thermal conduction
problem for a dyke at a magma chamber (modelled as a magma
lled edge crack in a half-space) are
T = T
m
Z = 0,
T = T
s
Z = h,
T = T
m
X = 0, 0 Z a
T
X
= 0, X = 0, Z a.
(A1)
The temperature in the host rock can be expressed as
T(X, Z) = T
(1)
(Z) + T
(2)
(X, Z), (A2)
where T
(1)
(Z) is the temperature without considering the dyke and
is given as
T
(1)
(Z) = T
m

T
m
T
s
h
Z. (A3)
The basic equation and the boundary conditions for T
(2)
(Z) are

2
T
(2)
Z
2
+

2
T
(2)
X
2
= 0,
T
(2)
= 0 Z = 0, X 0,
T
(2)
=
T
m
T
s
h
Z, X = 0, 0 < Z < a,
T
(2)
X
= 0 X = 0, Z > a.
(A4)
C
2011 The Authors, GJI, 186, 10951103
Geophysical Journal International
C
2011 RAS
Dyke propagation with viscoelastic properties 1103
The solution of the above problem can be obtained using Fourier
transform as follows;
T
(2)
(X, Z) =
_

0
a()e
X
sin( Z) d,
a() =
2

_
a
0
f (Z

) sin( Z

) dZ

,
(A5)
where f (Z) is dened by
T
(2)
X
= f (Z) X = 0, 0 < Z < a, (A6)
and satises the following integral equation

_
a
0
f (Z

) ln

+ Z
Z

dZ

=
T
m
T
s
h
Z 0 Z a. (A7)
The nal expression for T
(2)
(Z) has the following form:
T
(2)
(X, Z) =
1
2
_
a
0
ln
_
X
2
+(Z

+ Z)
2
X
2
+(Z

Z)
2
_
f (Z

) dZ

. (A8)
C
2011 The Authors, GJI, 186, 10951103
Geophysical Journal International
C
2011 RAS