- Email: [email protected]
On the dynamics of coupled grain size evolution and shear heating in lithospheric shear zones
Accepted Manuscript On the dynamics of coupled grain size evolution and shear heating in lithospheric shear zones Bradford J. Foley PII: DOI: Reference:
S0031-9201(18)30026-8 https://doi.org/10.1016/j.pepi.2018.07.008 PEPI 6174
To appear in:
Physics of the Earth and Planetary Interiors
Received Date: Accepted Date:
17 January 2018 13 July 2018
Please cite this article as: Foley, B.J., On the dynamics of coupled grain size evolution and shear heating in lithospheric shear zones, Physics of the Earth and Planetary Interiors (2018), doi: https://doi.org/10.1016/j.pepi. 2018.07.008
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
On the dynamics of coupled grain size evolution and shear heating in lithospheric shear zones Bradford J. Foley∗ Department of Geosciences, Pennsylvania State University, University Park PA
Abstract The formation of weak, localized shear zones in the lithosphere via tectonic forces is a requirement for the operation of plate tectonics on Earth. However, how shear zones form is not well understood. Shear heating and grain size reduction are two commonly cited mechanisms, and both processes are likely to be active in lithospheric shear zones. However, shear heating and grain size reduction can potentially act against each other, as elevated temperatures increase grain growth rates, and fast grain growth can eliminate weak shear zones formed via grain size reduction. Here, an effectively one-dimensional simple shear model that considers shear heating and grain size evolution is used to determine how grain size reduction and shear heating interact. The results indicate that for constant velocity boundary conditions at typical tectonic strain rates and mid- to lower-lithosphere temperatures, grain size reduction suppresses shear heating by weakening the shear zone. As a result shear heating at these conditions is small (∼ 10 K or less), unless grain size reduction via damage is significantly less effective than the models assume. ∗
Corresponding author Email address: [email protected] (Bradford J. Foley)
Preprint submitted to Physics of the Earth and Planetary Interiors
July 17, 2018
For constant stress boundary conditions, grain size reduction enhances shear heating and lowers the stress threshold required to trigger thermal runaway instabilities, as previous studies have found. Finally, combining shear heating and grain size reduction allows weak shear zones to form over a wider range of conditions than would result from either mechanism acting in isolation; coupling between shear heating and grain size reduction does not impede weak shear zone formation. Keywords:
1
1. Introduction
2
Why Earth is alone among the terrestrial planets of our solar system in
3
operating under a plate tectonic, rather than a stagnant lid, mode of man-
4
tle convection has long been an open question. In order for plate tectonics
5
to operate, mantle convective forces must be capable of forming weak, lo-
6
calized shear zones in the mid- to lower-lithosphere that effectively act as
7
plate boundaries (e.g. Tackley, 2000a; Bercovici et al., 2015). Without such
8
shear zones, the lithosphere is too strong to sink under its own weight and
9
form subduction zones. In particular, lithospheric strength is highest in the
10
mid- to lower-lithosphere, where deformation is predominantly ductile (e.g.
11
Kohlstedt et al., 1995); here the high viscosity of the lithosphere, due to the
12
temperature-dependence of rock rheology, produces a strong viscous resis-
13
tance to plates sliding past each other, or subducting beneath each other.
14
Forming low viscosity, ductile shear zones in this region is thus essential for
15
plate tectonics to operate. Moreover, localization of deformation into narrow
16
shear zones, that leave the interiors of plates stable and non-deforming, is
2
17
a defining characteristic of plate tectonics. Indeed, localized, ductile shear
18
zones are commonly seen in exhumed sections of the lower crust or litho-
19
sphere on Earth (e.g. Poirier, 1980; White et al., 1980; Hanmer, 1988; Jin
20
et al., 1998; Warren & Hirth, 2006; Skemer et al., 2010; Linckens et al., 2011).
21
However, how such shear zones form is not well understood.
22
Weakness and localization are linked, and both are necessary for the op-
23
eration of plate tectonics, but there are different requirements for shear zone
24
weakening and shear localization. Weak ductile shear zones can form via any
25
mechanism that allows deformation to induce a reduction in rock viscosity,
26
while localization additionally requires that stress decreases with increas-
27
ing strain rate (Mont´esi & Zuber, 2002; Bercovici & Karato, 2003). This
28
study focuses primarily on the processes that control shear zone strength.
29
Many mechanisms have been proposed for forming weak shear zones, in-
30
cluding shear heating (e.g. Yuen et al., 1978; Fleitout & Froidevaux, 1980;
31
Kaus & Podladchikov, 2006; Thielmann & Kaus, 2012), the development
32
of aligned mineral fabrics (e.g. Tommasi et al., 2009; Mont´esi, 2013), and
33
grain size reduction (e.g. Braun et al., 1999; Bercovici & Karato, 2003; War-
34
ren & Hirth, 2006; Landuyt & Bercovici, 2009; Bercovici & Ricard, 2012;
35
Platt, 2015; Bercovici & Ricard, 2016; Mulyukova & Bercovici, 2017). In
36
the above mechanisms deformation either causes dissipative heating, which
37
lowers viscosity through the temperature-dependence of mantle and crustal
38
rheology (e.g. Karato & Wu, 1993; Hirth & Kohlstedt, 2003), grain size re-
39
duction, which leads to weakening when deformation occurs in a grain size
40
sensitive creep regime such as diffusion creep or grain boundary sliding (e.g.
41
Hirth & Kohlstedt, 2003; Mont´esi & Hirth, 2003), or formation of mineral
3
42
alignment fabrics or foliations that cause viscous anisotropy, where viscos-
43
ity is lower when shearing in the direction minerals are aligned along (e.g.
44
Mont´esi, 2007; Hansen et al., 2012b).
45
Ultimately, in order to understand the origin and evolution of plate tec-
46
tonics on Earth, these shear zone weakening mechanisms must be incorpo-
47
rated into mantle convection models. While previous models have typically
48
invoked failure of the lithosphere as a result of pseudoplastic yielding for
49
weak shear zone formation (e.g. Moresi & Solomatov, 1998; Tackley, 2000b;
50
van Heck & Tackley, 2008; Foley & Becker, 2009), more sophisticated mecha-
51
nisms such as shear heating (Thielmann & Kaus, 2012) and grain size reduc-
52
tion (e.g. Landuyt et al., 2008; Foley et al., 2012; Bercovici & Ricard, 2014)
53
have been recently considered. However, most previous work has focused
54
on individual weakening mechanisms operating in isolation, when in reality
55
multiple mechanisms are likely to be active simultaneously, and thus may
56
interact in potentially unexpected ways.
57
In particular, this study focuses on coupling between grain size reduction
58
and shear heating, because these processes can potentially offset each other,
59
and thus impede the formation of weak shear zones (e.g. Kameyama et al.,
60
1997). Grain size reduction is a well supported mechanism for plate bound-
61
ary formation, as it is commonly seen in exhumed shear zones in the form
62
of mylonites and ultra-mylonites (e.g. White et al., 1980; Warren & Hirth,
63
2006). However, grains grow faster at warmer temperatures than cooler tem-
64
peratures (e.g. Evans et al., 2001), meaning that any shear heating induced
65
by flow in the shear zone may lead to faster grain growth, larger grains, and
66
viscous strengthening of the shear zone, if it deforms via grain size sensitive
4
67
creep. Thus, determining how temperature and grain size co-evolve in a shear
68
zone is essential for determining the conditions under which mylonites form,
69
and more generally what processes control the strength of plate boundaries
70
on Earth.
71
The combined effects of shear heating and grain size evolution have been
72
studied previously. Kameyama et al. (1997) shows that shear heating largely
73
eliminates grain size reduction during simple shear, unless grain growth in
74
the lithosphere is slower than experiments on pure olivine samples indicate.
75
However, grain growth in the lithosphere likely is slower than what is observed
76
in monominerallic olivine experiments, due to grain boundary pinning effects
77
from secondary phases (e.g. Evans et al., 2001; Hiraga et al., 2010). Moreover,
78
secondary phases not only impede grain growth, but also enhance grain size
79
reduction (Bercovici & Ricard, 2012; Farla et al., 2013; Linckens et al., 2015;
80
Cross & Skemer, 2017). Thielmann et al. (2015) & Thielmann (2017) also
81
studied coupled shear heating and grain size reduction in the context of deep
82
earthquake formation, and included some of the effects of secondary phases
83
on grain size evolution. However, as Thielmann et al. (2015) & Thielmann
84
(2017) focus on deep earthquakes induced by runaway shear heating instabil-
85
ities, they do not study how shear heating and grain size reduction interact
86
in stably flowing, ductile shear zones, a major goal of this paper. Moreover,
87
this paper more explicitly includes the effects of secondary phases, using a
88
grain size evolution theory based on Bercovici & Ricard (2012), and uses dy-
89
namic recrystallization experiments to constrain the key model parameters
90
in the theory, as in Mulyukova & Bercovici (2017).
91
In particular, I use an effectively one-dimensional model of a shear zone
5
92
undergoing simple shear to determine the conditions under which shear heat-
93
ing is negligible, and those where it is significant, in terms of the applied stress
94
or velocity driving deformation in the shear zone, the ambient lithospheric
95
temperature, shear zone width, and the parameters governing grain size evo-
96
lution. When shear heating is negligible, grain size reduction is unaffected
97
by heating and controls the shear zone strength, while when shear heating is
98
significant it impacts the grain size in the shear zone and the overall strength.
99
Scaling laws for the steady-state grain size and temperature are developed
100
for both constant stress and constant velocity boundary conditions, and used
101
to determine the conditions where shear heating becomes significant. The
102
models are also used to determine how coupled shear heating and grain size
103
evolution influences shear zone strength, in comparison to the cases where
104
grain size reduction or shear heating are considered in isolation.
105
The paper is organized as follows: §2 gives the theoretical formulation of
106
the shear zone model, §3 gives steady-state model results and scaling anal-
107
yses for both constant velocity and constant stress boundary conditions, §4
108
shows model results and scaling analysis for the time evolution of shear zone
109
development, §5 discusses the implications of the results for plate boundary
110
formation and the possible effects of uncertainties in the model formulation,
111
and conclusions are summarized in §6.
112
2. Theory
113
2.1. Governing equations
114
A two-dimensional shear zone is modeled (Figure 1). The center lies at
115
x = 0 and boundary at x = L/2, where L is the total width of the shear 6
116
zone, L/2 is the half width, and x measures distance from the shear zone
117
center (or core). Only half of the shear zone is modeled due to symmetry.
118
The shear zone boundary has a fixed temperature, equal to the background
119
temperature in the lithosphere of Tl , and is held at either a fixed velocity, U ,
120
or stress, τl . In reality, different thermal boundary conditions could apply, as
121
the point where temperature converges to the background temperature could
122
lie beyond the edge of the deforming shear zone. As discussed in §5.3, the
123
model setup may thus underestimate the degree of shear heating. The shear
124
zone is assumed to have no pressure gradients and no flow in the x-direction,
125
and to be infinitely long in the y-direction. As a result, the velocity profile is
126
also assumed to be uniform in y, and thus the shear zone model is effectively
127
one-dimensional (i.e. there is no variation in temperature, grain size, or
128
velocity the y-direction). With these assumptions, as well as assuming the
129
viscosity is isotropic (see §2.2), conservation of momentum reduces to ∂τxy = 0, ∂x
(1)
130
where τxy is the deviatoric shear stress in the shear zone. Equation (1)
131
implies that shear stress is constant in x in the shear zone. Using the viscous
132
constitutive law τxy = 2µε˙xy = µ
∂w , ∂x
(2)
133
where µ is viscosity, ε˙xy is the strain rate, and w is the velocity in the shear
134
zone in the y-direction. Combining (1) and (2) gives ∂ ∂x
∂w µ = 0. ∂x 7
(3)
y"
At"x"="0:" w"="0" dT/dx"="0"
0"
At"x"="L/2:"" T"="Tl" w"="U"or" τ"="τl"
L/2"
x"
Figure 1: Sketch of the shear zone model setup, including a schematic velocity profile.
135
Both temperature and grain size vary within the shear zone. Evolution of
136
temperature, T , is governed by energy conservation, where T is assumed to
137
vary only in the x-direction, resulting in no advective heat flow in the shear
138
zone as a result of the assumptions for velocity outlined above, and with
139
no heat sources aside from shear heating. In reality grain growth converts
140
surface energy into heating, and acts as a heat source, but this effect is found
141
to be small for grains larger than 10−8 m (Thielmann, 2017), and is neglected
142
here. With these assumptions temperature evolves as (1 − f )Ψ ∂T ∂ 2T =κ 2 + ∂t ∂x ρCp
(4)
143
where t is time, κ is thermal diffusivity, Ψ is the rate of deformational work,
144
f is the fraction of deformational work that goes into grain size reduction
145
(and thus 1 − f is the fraction available to drive shear heating), ρ is density,
146
and Cp is the specific heat (see Table 1 for a list of all variables and Table 2 8
147
for a list of parameters and their assumed values). The deformational work
148
rate is defined as Ψ = τij ε˙ij . For a one-dimensional shear zone undergoing
149
simple shear with isotropic viscosity, as modeled here, the only non-zero
150
components of τij are τxy = τyx , and for ε˙ij are ε˙xy = ε˙yx . As a result
151
Ψ = 2τxy ε˙xy , or Ψ = 2τ ε, ˙ where τ and ε˙ are the second invariants of the
153
deviatoric stress and strain-rate tensors, respectively; these quantities are p p defined as τ = (1/2)τij τij and ε˙ = (1/2)ε˙ij ε˙ij .
154
2.2. Rheology
152
155
A purely viscous, composite rheology that includes both dislocation and
156
diffusion creep is assumed. The purely viscous rheology is chosen as this study
157
focuses on ductile shear zone formation in the lower lithosphere. The model
158
therefore ignores elasticity, which can lead to thermal runaway instabilities
159
during shear zone formation, even with constant velocity boundary condi-
160
tions, as built up elastic stress is released (Regenauer-Lieb & Yuen, 1998;
161
Kelemen & Hirth, 2007; Thielmann et al., 2015). However, at temperatures
162
greater than ∼ 1000 K, as are expected in the mid- to lower-lithosphere, and
163
at typical geologic strain rates, elasticity is negligible (Thielmann, 2017), and
164
thus the purely viscous rheology used in this study is justifiable. Note also
165
that by including models with constant stress boundaries, thermal runaway
166
instabilities are still captured in the purely viscous models presented here,
167
and thus some first order comparisons to viscoelastic models can be made
168
(see §3.1.2 & 5.2).
169
The composite rheology used in this study is given by: ε˙ = Bdis τ n + Bdif Am τ 9
(5)
170
where Bdis is the compliance for dislocation creep, Bdif the compliance for
171
diffusion creep, A is the fineness, or inverse grain size, and n and m are
172
constants. I use n = 3, rather than n = 3.5 as given in Hirth & Kohlstedt
173
(2003), as an analytical expression for stress as a function of strain rate can
174
not be obtained with n = 3.5, and such an expression is necessary for the
175
scaling analysis performed in §3.2 and Appendix B. Using n = 3 instead
176
of n = 3.5 does not significantly change the strain rate obtained from the
177
dislocation creep flow law at a stress of 1 MPa, which is in the middle of the
178
range of stresses covered in this study. I also use m = 3 as appropriate for
179
Coble creep (e.g. Hirth & Kohlstedt, 2003). The dislocation and diffusion
180
creep compliances are temperature dependent, following Bdis
181
Edis = Bdis0 exp − RT
(6)
and
Bdif
Edif = Bdif 0 exp − RT
(7)
182
where Edis and Edif are the activation energies for dislocation and diffusion
183
creep, respectively, and R is the universal gas constant. I use Edis = 530 kJ
184
mol−1 and Edif = 375 kJ mol−1 (Hirth & Kohlstedt, 2003). The composite
185
rheology will be controlled by whichever mechanism, diffusion or dislocation
186
creep, leads to the larger strain rate. The field boundary between deformation
187
controlled by diffusion creep versus that controlled by dislocation creep can
188
be found by equating the strain rates from the two mechanisms. Writing this
10
189
boundary in terms of fineness gives the field boundary fineness, Af ield , as Af ield =
Bdis0 τ n−1 Bdif 0
m1
exp
Edif − Edis . mRT
(8)
190
When fineness is larger than Af ield deformation is in the diffusion creep
191
regime, and in the dislocation creep regime when fineness is lower than Af ield .
192
2.3. Damage theory
193
The evolution of fineness, A, is modeled using a simplified version of
194
the two-phase damage theory developed in Bercovici & Ricard (2012) and
195
extended in Bercovici & Ricard (2016). Bercovici & Ricard (2012) shows
196
that deformation of a mixture of phases, predominately olivine and pyroxene
197
in mantle peridotite, quickly enters a “pinned” state, where the grain sizes
198
of both phases are controlled by the curvature of the interface between the
199
phases, or the interface roughness. In the pinned state the grain size is
200
equal to the interface roughness multiplied by the factor π/2; the factor
201
(π/2)m is thus incorporated into Bdif 0 in the diffusion creep compliance in
202
this study. With the above assumptions, the fineness, A, evolves as (e.g.
203
Foley & Bercovici, 2014; Foley et al., 2014) dA fΨ = − hAp dt γ
(9)
204
where γ is the surface tension, h is the grain growth, or healing, rate, and p
205
is a constant (p = 5 is used, as estimated for the pinned state by Bercovici
206
& Ricard (2012)). Note that p in the above formulation is analogous to q in
207
the two-phase theory of Bercovici & Ricard (2012); p and q are related as 11
208
p = q +1. The healing rate is temperature-dependent, following (e.g. Karato,
209
1989; Evans et al., 2001) Eh h = hn exp − , RT
(10)
210
where hn is a pre-exponential constant and Eh is the activation energy for
211
grain growth. The first term on the right hand side of (9) represents grain
212
size reduction (fineness growth) as a result of partitioning a fraction, f , of the
213
deformational work into surface energy (f is referred to here as the damage
214
partitioning fraction). The second term represents grain growth (fineness
215
reduction) which acts to reduce the surface energy in the rock. A steady-state
216
fineness exists where grain size reduction and grain growth are in balance.
217
Note that in this formulation, where grain size is assumed to be pinned by
218
the interface between the two mineral phases that make up the bulk rock,
219
both dislocation and diffusion creep contribute to the deformational work
220
that drives grain size reduction. While dynamic recrystallization can only
221
occur via dislocation creep (e.g. Karato, 2008), damage to, or warping of,
222
the interface between phases can occur via any solid-state creep mechanism
223
(Bercovici & Ricard, 2012).
224
The bulk grain growth, or healing, rate in a rock in the pinned state is
225
calibrated to grain growth experiments. The pre-exponential constant, hn ,
226
is calculated in terms of a standard grain growth pre-exponential constant,
227
h∗n , valid for grain growth without pinning, from the following relationship
228
(Bercovici & Ricard, 2013): h∗n (p − 1) ∗ p−3 hn = r . 500 12
(11)
Table 1: Model Variables
Variable τxy x µ ε˙xy w T t Ψ τ ε˙ Bdis Bdif A h
Meaning Shear stress Distance from shear zone center Viscosity Strain rate Velocity Temperature Time Deformational work Second invariant of stress tensor Second invariant of strain rate tensor Dislocation creep compliance Diffusion creep compliance Fineness (inverse grain size) Healing rate
Units Pa m Pa s s−1 m s−1 K s W m−3 Pa s−1 −n Pa s−1 mm Pa−1 s−1 m−1 p−1 m s−1
Equation (1) (1) (2) (2) (2) (4) (4) (4) below (4) below (4) (6) (7) (9) (10)
229
Here r∗ is the interface roughness where the pinned state is reached; for a
230
mixture of olivine and pyroxene typical of peridotite r∗ = 10−6 m (Bercovici
231
& Ricard, 2013). As in Foley & Bercovici (2014) and Foley et al. (2014), h∗n
232
is calculated for a given grain growth activation energy from the standard
233
grain growth rate of 3 × 10−15 m2 s−1 at T = 1630 K given in Bercovici &
234
Ricard (2012), yielding h∗n =
3 × 10−15 . −Eh exp R1630
(12)
235
In the pinned state, f represents the fraction of deformational work
236
that partitions into damage along the interface between phases, while Eh
237
is, strictly speaking, the activation energy for coarsening of the interface be-
238
tween phases; as pining then controls grain size, these quantities effectively
239
govern grain size evolution. However, values of these parameters for dam-
240
age to the grains themselves, when pinning is absent, could be different. I
241
thus use f ∗ , h∗ , and Eh∗ to represent damage partitioning fraction, healing 13
242
rate, and grain growth activation energy for grains in the absence of pinning.
243
Values for f and Eh are not well known, but dynamic recrystallization ex-
244
periments can be used to constrain f ∗ or Eh∗ (Rozel et al., 2011; Mulyukova
245
& Bercovici, 2017). In particular, laboratory experiments show a trend in
246
recrystallized grain size with temperature, which constrains the temperature
247
dependence of the ratio f ∗ /h∗ , but does not constrain f ∗ and Eh∗ uniquely
248
(see Appendix A). With a chosen Eh∗ , f ∗ can be inferred from experimental
249
results and is generally a function of temperature, such that the temperature
250
dependence of f ∗ /h∗ matches the experiments (Rozel et al., 2011; Mulyukova
251
& Bercovici, 2017). In this study I chose a value of Eh∗ that results in f ∗ in-
252
dependent of temperature, in order to simplify the model. This is found at
253
Eh∗ ≈ 430 kJ mol−1 , consistent with the value used in Mulyukova & Bercovici
254
(2017). With this Eh∗ and h∗n as given in (12), f ∗ ≈ 10−4 (see Appendix A).
255
The simplest assumption is then that f = f ∗ and Eh = Eh∗ , and I make this
256
assumption for Eh in this study. Mulyukova & Bercovici (2017) suggests f
257
could be lower than f ∗ , finding that f ∼ 10−5 f ∗ can fit field observations of
258
exhumed mylonites. I thus test a range of f = 10−4 − 10−9 in this study.
259
The chosen value of Eh gives hn ≈ 1.4 × 10−15 m4 s−1 .
260
2.4. Numerical solution scheme
261
Equations (4) and (9) are solved to obtain the one-dimensional profiles of
262
temperature and fineness through the shear zone as functions of time. The
263
shear zone is discretized using 401 grid points in the x−direction. Solution
264
of (9) is straightforward because fineness grid points are independent of each
265
other in the grain size evolution formulation employed in this study, where
266
fineness evolution is assumed to be only a function of the local strain rate, 14
Parameter L Tl U τl κ f ρ Cp n m Bdis0 Edis Bdif 0 Edis R γ p hn ε˙0
Table 2: Model Parameters and Assumed Values Meaning Assumed value Total shear zone width 2-200 km, reference value: 20 km Background lithosphere temperature 900-1300 K, reference value: 1000 K Imposed shear zone boundary velocity Calculated from ε˙0 (m s−1 ) Imposed shear zone boundary stress 104 − 109 Pa −6 Thermal diffusivity 10 m s−1 (ref. 4) −9 −4 Damage partitioning fraction 10 − 10 , reference value: 10−4 (ref. 2) Lithosphere density 2800 kg m−3 (ref. 4) Heat capacity 1000 J kg−1 K−1 (ref. 4) Stress exponent 3 (ref. 3) grain size sensitivity exponent 3 (ref. 1) Dislocation creep compliance constant 1.1 × 10−13 Pa−n s−1 (ref. 1) Dislocation creep activation energy 530 kJ mol−1 (ref. 1) Diffusion creep compliance constant 1.5 × 10−15 mm Pa−1 s−1 (ref. 1) Diffusion creep activation energy 375 kJ mol−1 (ref. 1) Universal gas constant 8.314 J mol−1 (ref. 4) Surface tension 1 J m−2 (ref. 5) grain growth exponent 5 (ref. 2) grain growth pre-exponential constant 1.4 × 10−15 m4 s−1 (ref. 2) Initial strain rate 10−18 − 10−9 s−1
References for chosen parameter values: ref. 1) Hirth & Kohlstedt (2003); ref. 2) derived in this paper; ref. 3) n = 3.5 in Hirth & Kohlstedt (2003), but n = 3 greatly simplifies model analysis and is chosen here. Using n = 3 instead of n = 3.5 does not significantly change the values of Bdis0 and Edis (Mulyukova & Bercovici, 2017); ref. 4) Turcotte & Schubert (2002); ref. 5) Mulyukova & Bercovici (2017)
15
Equation above (1) above (1) above (1) above (1) (4) (4),(9) (4) (4) (5) (5) (6) (6) (7) (7) (7) (9) (9) (10) below (13)
267
temperature, and stress. Thus (9) is solved at each grid point using MAT-
268
LAB’s ordinary differential equation (ODE) solver “ode15s.” The equation
269
for temperature is more difficult, as it contains a diffusion term. To solve
270
(4) the method of lines is used (Schiesser, 2012), and the resulting system
271
of ODEs is also solved with ode15s in MATLAB. As introduced in §2.1, the
272
boundary at x = 0 corresponds to the shear zone center, and at x = L/2 to
273
the shear zone edge. Due to symmetry in temperature at the shear zone cen-
274
ter, ∂T /∂x = 0 at x = 0. The shear zone edge is held at a fixed temperature
275
Tl , and thus T = Tl at x = L/2.
276
277
For a constant stress boundary condition at the shear zone edge, Ψ can be written in a straightforward manner as Ψ = 2τl ε˙ = 2Bdis τln+1 + 2Bdif Am τl2 .
(13)
278
However, for a constant velocity boundary condition, the stress varies with
279
time as temperature and fineness evolve. For constant velocity boundary
280
conditions, the initial strain rate in the shear zone, ε˙0 , is specified, and the
281
constant edge velocity U = Lε˙0 /2. For a given stress, strain rate as a function
282
of x can be obtained for a known temperature and fineness profile in the shear
283
zone from (5). Integrating ε(x) ˙ from x = 0 to x = L/2 gives the velocity at
284
the edge of the shear zone for a given stress. Thus, stress can be determined
285
for the constant velocity shear zones by finding the stress that satisfies Z
x=L/2
ε(x)dx ˙ − U = 0.
(14)
x=0
286
The MATLAB root finding function “fzero” is used to find τ that satisfies 16
287
(14) at each timestep, and Ψ is then determined from (13). The shear zone
288
is assumed to initially have a uniform temperature equal to Tl and uniform
289
fineness equal to A0 , the reference fineness; A0 = 100, corresponding to a
290
reference grain size of 1 cm, is used.
291
When high temperatures are reached as a result of shear heating, melting
292
can occur. Melting is not explicitly included in the models, which assume a
293
solid state rheology. A reasonable estimate for the melting temperature in
294
the lower lithosphere is ≈ 1600 K, which is the solidus temperature of olivine
295
at ≈ 50 km depth based on the melting curve of Hirschmann (2000). In
296
the figures shown below, models that yield temperatures above this nominal
297
melting temperature are highlighted.
298
3. Results: Steady-state fineness and temperature
299
3.1. Numerical model results
300
3.1.1. Constant velocity shear zones
301
I first consider results from shear zones with constant velocity boundary
302
conditions. Figure 2 shows typical fineness and temperature profiles through
303
shear zones at steady-state. These steady-state profiles are found to be stable
304
to perturbations in fineness (perturbations up to 100 times the steady-state
305
fineness over regions extending up to 0.1L from the shear zone center were
306
tested). At low initial strain rates, fineness is constant spatially across the
307
shear zone, and the temperature rise due to shear heating is negligible (Fig-
308
ures 2A & B). As a result, strain rate and deformational work rate are also
309
constant spatially, and there is no localization in the shear zone (Figures 2C
310
& D). At steady-state and when diffusion creep is the dominant deforma17
311
tion mechanism (as Figures 2 & 3 show is the case for the majority of the
312
model results), fineness is a monotonic function of stress (e.g. Foley et al.,
313
2012; Foley & Bercovici, 2014; Foley & Driscoll, 2016). As stress is constant
314
spatially as required by (1), fineness will also be constant and localization
315
will not occur, when shear heating is negligible. Note that although local-
316
ization does not occur for the simple damage formulation used in this study,
317
more sophisticated grain size reduction theories do allow for localization (see
318
§5.3). Higher initial strain rates then increase the deformational work rate,
319
and hence lead to higher fineness in the shear zone (Figure 2B).
320
Eventually, initial strain rates are high enough for shear heating to be-
321
come significant. In this situation the peak temperature (Tp ) occurs at the
322
shear zone center, away from the cooler fixed temperature boundaries, while
323
the opposite is true for fineness (Figures 2A & B). High temperatures en-
324
hance grain growth and decrease fineness; thus the smallest fineness occurs
325
at the shear zone center, while fineness grows towards the shear zone edge as
326
temperatures cool. The deformational work rate is also larger at the shear
327
zone core than at the margins (Figure 2D), but the effect of temperature is
328
stronger and results in the documented grain coarsening at the shear zone
329
center. When shear heating is significant, localization can occur when the
330
temperature rise is on the order of the background temperature (e.g. Bercovici
331
& Karato, 2003), as seen in the case with an initial strain rate of 10−10 s−1
332
(Figures 2C & D). Furthermore, the fineness throughout the shear zone is
333
smaller at ε˙0 = 10−10 s−1 than at ε˙0 = 10−12 s−1 , as a result of shear heating
334
inducing grain growth.
18
1200 1000 800
Strain Rate [s-1]
Fineness [m-1]
1400
A
105 104 103 102 101 100
C
10-10 10-12 10-14 10-16 10-18
0 1 2 3 4 5 6 7 8 9 10
Distance [km]
Def. Work Rate [Pa s-1]
Temperature [K]
106 1600
B
D
10-4 10-6 10-8 10-10 10-12 10-14
0 1 2 3 4 5 6 7 8 9 10
Distance [km] Strain Rate = 10-18 s-1 Strain Rate = 10-15 s-1 Strain Rate = 10-12 s-1 Strain Rate = 10-10 s-1
Figure 2: Temperature (A), fineness (B), strain rate (C), and deformational work rate (D) as a function of distance from the shear zone center (x). L/2 = 10 km, Tl = 1000 K, f = 10−4 , and other parameters as given in Table 2. Four different initial strain rates, ε˙0 , are used: ε˙0 = 10−18 s−1 (solid line), ε˙0 = 10−15 s−1 (dashed line), ε˙0 = 10−12 s−1 (dot-dashed line), ε˙0 = 10−10 s−1 (dotted line). Blue lines give the field boundary fineness for each initial strain rate, where fineness values larger than the field boundary result in deformation in the diffusion creep regime.
19
Strain Rate [s-1]
Strain Rate [s-1]
10-16 10-14 10-12 10-10
10-16 10-14 10-12 10-10
107
10-6 10-8
100 10-2 10-4 L/2 = 1 km L/2 = 10 km L/2 = 100 km
10-6 10-8 104
102
107
100 10-2 10-4 -4
10-8
105 103
f = 10 f = 10-6 f = 10-9
10-16 10-14 10-12 10-10 -1
Strain Rate [s ]
106
L/2 = 1 km L/2 = 10 km L/2 = 100 km
105 104 C
f = 10-4 f = 10-6 f = 10-9
106 105 104
103
G
108
J
104
L/2 = 1 km L/2 = 10 km L/2 = 100 km
107
103
105
102 101
106
108
F
106 104
Tl = 1100 K Tl = 1000 K Tl = 900 K
107
103
101
I
102
10-6
102
107
E
102
Tl = 1100 K Tl = 1000 K Tl = 900 K
103 101
Fineness [m-1]
104
Stress [Pa]
Tl = 1100 K Tl = 1000 K Tl = 900 K
105 104
Stress [Pa]
10-4
10
108
B
6
Stress [Pa]
100 10-2
Fineness [m-1]
A
102
107 106 105 104 103
10-16 10-14 10-12 10-10
-4
K
f = 10 f = 10-6 f = 10-9
10-16 10-14 10-12 10-10
Strain Rate [s-1]
-1
Strain Rate [s ]
Def. Work Rate [Pa s-1] Def. Work Rate [Pa s-1] Def. Work Rate [Pa s-1]
Strain Rate [s-1] 10-16 10-14 10-12 10-10
Fineness [m-1]
Temperature Rise [K] Temperature Rise [K] Temperature Rise [K]
104
Strain Rate [s-1] 10-16 10-14 10-12 10-10
Tl = 1100 K Tl = 1000 K Tl = 900 K
10-2 10-4 10-6 10-8 10-10 10-12
D
L/2 = 1 km L/2 = 10 km L/2 = 100 km
10-2 10-4 10-6 10-8 10-10 10-12
H
10-2 L 10-4 10-6 10-8 10-10 10-12
f = 10-4 f = 10-6 f = 10-9
10-16 10-14 10-12 10-10
Strain Rate [s-1]
Figure 3: Steady-state temperature rise (Tp − Tl ), fineness, stress, and deformational work rate at the shear zone core from the full numerical solution (circles) and scaling laws for the steady-state solutions, developed in §3.2.1, (lines). Background temperature is varied for a constant f = 10−4 and L/2 = 10 km (A, B, C, & D), shear zone width varied for constant f = 10−4 and Tl = 1000 K (E, F, G, & H), and f is varied for constant L/2 = 10 km and Tl = 1000 K (I, J, K, & L). Symbols representing numerical model results and lines representing the scaling laws are colored to represent Tl , L, or f as indicated in the legend for each figure. Open symbols and dotted lines denote model results where shear heating leads to peak temperature above the nominal melting temperature. The dashed lines in the figures for fineness as a function of strain rate mark the field boundary between dislocation and diffusion creep, where fineness larger than the field boundary means deformation takes place in the diffusion creep regime.
20
335
Thus, the temperature rise (Tp − Tl ) and fineness at the shear zone cen-
336
ter in steady-state show two clear trends with initial strain rate; one at low
337
strain rates where shear heating is negligible, and one at higher strain rates
338
where shear heating is significant (Figure 3). Fineness increases with strain
339
rate until a critical threshold is reached, where the high temperatures caused
340
by shear heating counteract damage, and steady-state fineness begins to de-
341
crease with increasing strain rate. Eventually temperatures become so high
342
that melting would occur, meaning further grain growth and temperature
343
increases beyond this point would not take place in nature. Temperature
344
rise increases with increasing strain rate for all strain rates, but the slope
345
becomes smaller at the initial strain rate threshold where fineness begins
346
to decrease. Here, significant shear heating causes the viscosity to decrease
347
sharply with increasing strain rate, such that stress also decreases and the
348
deformational work rate becomes only a weakly increasing function of strain
349
rate (Figure 3). As the temperature rise is controlled by Ψ, temperature rise
350
also becomes a weakly increasing function of strain rate when shear heating
351
is significant. That stress decreases with increasing strain rate when shear
352
heating is significant also explains why localization can happen in this case,
353
as Mont´esi & Zuber (2002) shows that stress being inversely proportional to
354
strain rate is a key requirement for shear localization.
355
Changing the shear zone background temperature has a minor effect on
356
temperature rise, but a large influence on steady-state fineness when shear
357
heating is negligible (Figure 3A & B); when shear heating is significant the
358
background temperature has no strong influence on both temperature rise
359
and shear zone core fineness. When shear heating is negligible, lower back-
21
360
ground temperatures suppress healing and therefore allow smaller steady-
361
state grain sizes (larger fineness) to be reached. However, when shear heat-
362
ing is significant the temperature in the shear zone is controlled by shear
363
heating rather than the background temperature, and thus grain sizes ap-
364
proach the same value for different background temperatures. Varying shear
365
zone width has a strong influence on the temperature rise, and, when shear
366
heating is significant, on the shear zone core fineness as well (Figure 3E &
367
F). The larger L, the larger the temperature rise because the temperature
368
gradient across the shear zone is smaller. As a result, when shear heating is
369
significant, fineness is lower for wider shear zones, because the higher tem-
370
peratures enhance grain growth. When shear heating is negligible, fineness
371
is unaffected by shear zone width because in these models the shear zone
372
boundary velocity is varied with L to produce a given initial strain rate (i.e.
373
varying L does not change the initial strain rate). Increasing the amount
374
of deformational work that goes into grain size reduction, f , leads to larger
375
fineness in the shear zone. As a result, the effective viscosity of the shear
376
zone is decreased. As Ψ = 2τ ε˙0 and τ = 2µε˙0 , where µ is the effective vis-
377
cosity, decreasing µ decreases Ψ. Increasing f therefore also decreases the
378
temperature rise (Figures 3I-L). Grain size reduction thus suppresses shear
379
heating by weakening the shear zone.
380
3.1.2. Constant stress shear zones
381
With a constant stress boundary condition, a runaway shear heating in-
382
stability is possible above a critical driving stress, where heating leads to
383
higher strain rates in the shear zone, which induce further heating and in-
384
crease strain rates even more; the instability continues until the shear zone 22
385
experiences melting (e.g. Schubert & Yuen, 1978; Yuen & Schubert, 1979;
386
Kameyama et al., 1999; Turcotte & Schubert, 2002; Kaus & Podladchikov,
387
2006). As shown by Thielmann et al. (2015), grain size reduction lowers
388
the critical stress where runaway heating occurs, because grain size reduc-
389
tion weakens the shear zone, thereby enhancing deformational work rate and
390
shear heating. The results of this study confirm that grain size reduction low-
391
ers the threshold stress where thermal runaway occurs, as shown in §3.3.2.
392
In this section model results that reach a steady-state are shown; any mod-
393
els that undergo thermal runaway have no steady-state solution and are not
394
shown. As with the constant strain rate shear zones, tests of the stability of
395
the steady-states reached were performed. Steady-states are stable as long as
396
perturbations are not so large that they induce thermal runaway by weaken-
397
ing the shear zone. Any perturbations that do not trigger thermal runaway
398
decay and the unperturbed steady-state is recovered.
399
Consistent with the constant velocity case, shear localization does not
400
occur when shear heating is negligible (Figure 4). In fact, significant local-
401
ization is not seen in any steady-state solution, because large degrees of shear
402
heating are only found in models that experience thermal runaway, and thus
403
do not reach a steady-state. When thermal runaway is ongoing very high val-
404
ues of temperature and fineness in the shear zone core do lead to significant
405
shear localization.
406
Results for steady-state temperature rise and fineness at the shear zone
407
core (Figure 5) are similar to those from the constant velocity models at low
408
initial strain rates. Increasing stress leads to higher fineness and temperature
409
at the shear zone core, while lower background temperatures and higher
23
A
1010 1000
106
Fineness [m-1]
Temperature [K]
1020
C
10-12 10-14 10-16 -18
0 1 2 3 4 5 6 7 8 9 10
Distance [km]
Def. Work Rate [Pa s-1]
Strain Rate [s-1]
10-10
10-20
104 103
B
102 101 100
990
10
105
10-4
D
10-6 10-8 10-10 10-12 10-14 10-16
0 1 2 3 4 5 6 7 8 9 10
Distance [km] Stess = 104 Pa Stress = 105 Pa Stress = 4 × 105 Pa Stress = 7 × 105 Pa
Figure 4: Profiles of temperature (A), fineness (B), strain rate (C), and deformational work rate (D) across a shear zone with L/2 = 10 km, Tl = 1000 K, f = 10−4 , and other parameters as given in Table 2. Four different stresses, τl , are used: τl = 104 Pa (solid line), τl = 105 Pa (dashed line), τl = 4 × 105 Pa (dot-dashed line), τl = 7 × 105 Pa (dotted line). Blue lines give the field boundary fineness for each stress.
24
10-8
Tl = 1100 K Tl = 1000 K Tl = 900 K
10-10
100
E
10-2 10-4 10-6 10-8
L/2 = 1 km L/2 = 10 km L/2 = 100 km
10-10 102 100 10-2 10-6 10-8
f = 10-4 f = 10-6 f = 10-9
10-10 104
105
106
107
Stress [Pa]
105 Tl = 1100 K Tl = 1000 K Tl = 900 K
103 102
106
F
10-10 G 10-12 10-14 10-16 10-18 10-20 10-22 10-24
105 104
L/2 = 1 km L/2 = 10 km L/2 = 100 km
103 102 101
I
10-4
B
104
104 10-10 C 10-12 -14 10 10-16 10-18 10-20 10-22 10-24
107
101
Fineness [m-1]
102
106
106
f = 10-4 f = 10-6 f = 10-9
J
105 104 103 102 101
104
105
106
107
Stress [Pa]
10-10 K 10-12 10-14 10-16 10-18 10-20 10-22 10-24 104
105
106
107
Tl = 1100 K Tl = 1000 K Tl = 900 K
L/2 = 1 km L/2 = 10 km L/2 = 100 km
f = 10-4 f = 10-6 f = 10-9
105
106
107
Stress [Pa]
Stress [Pa] Def. Work Rate [Pa s-1] Def. Work Rate [Pa s-1] Def. Work Rate [Pa s-1]
10-6
106
105
Strain Rate [s-1]
10-4
104
Strain Rate [s-1]
10-2
Stress [Pa]
Stress [Pa]
107
Strain Rate [s-1]
106
Fineness [m-1]
105
Fineness [m-1]
Temperature Rise [K]
Temperature Rise [K]
Temperature Rise [K]
Stress [Pa] 104 102 A 100
104 10-4 D 10-6 10-8 10-10 10-12 10-14 10-16 10-18 10-20
105
107
Tl = 1100 K Tl = 1000 K Tl = 900 K
10-4 H 10-6 10-8 10-10 10-12 10-14 10-16 10-18 10-20 10-4 L 10-6 10-8 10-10 10-12 10-14 10-16 10-18 10-20 104
106
L/2 = 1 km L/2 = 10 km L/2 = 100 km
f = 10-4 f = 10-6 f = 10-9
105
106
107
Stress [Pa]
Figure 5: Same as Figure 3, but for constant stress boundary conditions rather than constant velocity. Scaling laws for the steady-state shear zone center temperature and fineness are developed in §3.2.2, and strain rate (C,G, & K) is strain rate at the shear zone core.
25
410
damage partitioning fractions lead to larger fineness. The effect of shear zone
411
width on temperature rise and fineness is also consistent with the constant
412
velocity models. However, one major difference is that with a constant stress
413
boundary condition, a lower shear zone viscosity leads to higher rates of
414
deformational work, rather than lower. At fixed τl , lower viscosity causes a
415
higher strain rate, and thus the product of stress and strain rate is higher. As
416
a result, increasing background temperature or f causes the temperature rise
417
to increase, rather than decrease as in the constant velocity models (Figure
418
5A, D, I & L). Thus, with constant stress boundary conditions, damage
419
enhances, rather than hinders, shear heating.
420
3.2. Scaling analysis
421
3.2.1. Constant velocity shear zones
422
The model results can be understood by a simple scaling analysis. The
423
full numerical models demonstrate that shear localization does not occur,
424
and thus deformational work rate, Ψ, is constant in x, when shear heating
425
is negligible (Figure 2). I thus assume Ψ is constant in space for this scal-
426
ing analysis. Localization can occur when shear heating is significant, but
427
the scaling laws derived here are found to provide a good approximation of
428
the numerical model results even when shear heating is important. With a
429
constant Ψ, (4) can be integrated to give the peak temperature at the shear
430
zone core, Tp : Tp = Tl +
(1 − f )ΨL2 . 8k
(15)
431
Deformational work is then defined at the peak temperature and steady-state
432
fineness in the shear zone core. The model results nearly all lie within the 26
433
diffusion creep regime as a result of extensive grain size reduction. I thus
434
neglect dislocation creep in the following derivation; scaling laws derived
435
using the full composite rheology are given in Appendix B. For constant
436
velocity boundary conditions, Ψ is written in terms of strain rate as Ψ = 2τ ε˙0 =
2ε˙20 . Bdif Am
(16)
437
The steady-state fineness at the shear zone center can then be found from
438
(16), and (7) - (10) as, As =
439
2f ε˙20 γBdif 0 hn
1 p+m
exp
Eh + Edif . RTp (p + m)
(17)
The equation for peak temperature, combining (15) - (17), and (7) is (1 − f )L2 Tp − Tl = 4k
2f γhn
m − p+m
ε˙20 Bdif 0
p p+m
exp
pEdif − mEh . (18) RTp (p + m)
440
(18) is a transcendental equation for Tp that can be solved numerically (in
441
this work the fzero function in MATLAB is used). With Tp found, the shear
442
zone core steady-state fineness is directly computed from (17).
443
The scaling laws for Tp and As are shown in Figure 3, and match the
444
model results well. In order to match the full range of results obtained from
445
the numerical shear zone models, the scaling laws derived in Appendix B with
446
a full composite rheology are shown in Figure 3. However, scaling laws with
447
a composite rheology only differ from those that neglect dislocation creep,
448
(17) & (18), for the limited range of conditions where deformation is in the
449
dislocation creep regime (e.g. at ε˙0 > 10−12 s−1 for the f = 10−9 models 27
450
in Figure 3J). At low strain rates, where shear heating is negligible, the
451
scaling laws and numerical model results match exactly. At high strain rates
452
the scaling laws underestimate Tp and As , because modest shear localization
453
in the shear zone core causes higher rates of deformational work in this
454
region than the scaling analysis predicts. Nevertheless, (17)-(18) provide a
455
good first order estimate of the shear zone core temperature and fineness at
456
steady-state.
457
When shear heating is negligible Tp ≈ Tl , and thus Tp can be replaced
458
with Tl in the right hand sides of (17) - (18), and Tp − Tl and As can be
459
computed directly. The trends in Tp and As with imposed strain rate can
460
then be easily understood from the scaling laws. Fineness scales with strain
461
rate as As ∼ ε˙0
462
rise scales as Tp −Tl ∼ ε˙0
463
strain rate than fineness, increasing strain rate will eventually cause shear
464
heating to become significant and impede grain size reduction, as the model
465
results show.
2/(p+m)
when shear heating is negligible, while the temperature 2p/(p+m)
. As temperature increases more sharply with
466
When Tp − Tl ∼ 10, the shear zone compliance, and hence effective
467
viscosity, begins to be significantly affected by shear heating, and Tp is no
468
longer approximately equal to Tl . The weak increase in Tp − Tl seen in
469
the numerical models at these high strain rate conditions is captured by the
470
scaling law for Tp . From (18), there is a competition between increasing strain
471
rate, which acts to increase shear zone temperature as described above, and
472
the resulting increase in temperature which acts to dampen shear heating by
473
decreasing viscosity and in turn stress. For parameters used in this study,
474
the result of this competition is the weak increase in temperature rise with
28
475
strain rate seen from the numerical models. Likewise the decrease in fineness
476
with strain rate when shear heating is significant is captured by (17) through
477
a similar competition; increasing strain rate acts to increase fineness, but
478
higher temperatures act to decrease fineness as a result of both faster grain
479
growth and a lower viscosity. In the end the temperature effect wins out,
480
resulting in the decreasing fineness trend seen in the models.
481
The scaling laws also demonstrate how other model parameters and shear
482
zone characteristics control the steady-state temperature and fineness. For
483
example, Tp − Tl scales as L2 , because a wider shear zone results in slower
484
diffusive heat transport, and thus a higher temperature. Fineness is not
485
directly influenced by L, which does not enter (17), but will be indirectly
486
influenced by variations in shear zone temperature caused by variations in L.
487
Shear zone center fineness increases with the damage to healing ratio, f /hn ,
488
as As ∼ (f /hn )1/(p+m) , because a higher damage to healing ratio means
489
more deformational work goes into reducing grain size and/or grain growth
490
is slower. Damage to healing ratio also enters into the scaling law for tem-
491
perature rise, as increasing fineness acts to dampen shear heating; (18) shows
492
that Tp − Tl ∼ (f /hn )−m/(p+m) .
493
3.2.2. Constant stress shear zones
494
Scaling laws for constant stress boundary conditions can be derived in a
495
similar fashion to those in §3.2.1. Again assuming that Ψ is constant in x
496
(Figure 4) and diffusion creep dominates, Ψ = 2Bdif Am τl2 .
29
(19)
497
Defining Ψ at the shear zone core temperature and fineness, As =
498
2f Bdif 0 τl2 γhn
1 p−m
exp
Eh − Edif RTp (p − m)
(20)
and (1 − f )L2 Tp = Tl + 4k
2f γhn
m p−m
τl2 Bdif 0
p p−m
exp
mEh − pEdif . RTp (p − m)
(21)
499
For most conditions where a steady-state solution exists, Tp ≈ Tl . Only at
500
stresses just before thermal runaway takes place does significant shear heating
501
of 10-100 K occur. Thus for the discussion of steady-state shear zones at low
502
stress in this section, Tp = Tl in the right hand side of (20) - (21) is assumed.
503
The scaling laws show that As ∼ τl
504
the increasing trends in As and Tp with stress seen in the numerical models.
505
Steady-state shear zone core fineness is an increasing function of f /hn , as in
506
the constant velocity case, scaling as As ∼ (f /hn )1/(p−m) . Furthermore, peak
507
temperature increases with increasing damage to healing ratio, scaling as
508
Tp ∼ (f /hn )m/(p−m) , because lowering shear zone viscosity via more extensive
509
damage increases the rate of deformational work (see §3.1.2). The preceding
510
reasoning also explains why grain size reduction lowers the threshold stress
511
where thermal instability occurs, as shown in Thielmann et al. (2015). The
512
dependences of fineness and peak temperature on background temperature
513
are the result of two competing effects. For fineness, increasing temperature
514
acts to decrease viscosity, and thus increase Ψ, but also increases the grain
515
growth rate. Thus the difference between the grain growth and diffusion
516
creep viscosity activation energies controls which of these two competing
2/(p−m)
30
2p/(p−m)
and Tp ∼ τl
, explaining
517
affects dominates, as (20) demonstrates. In the case of peak temperature,
518
the same competing effects of varying Tl are at play, and the sign of the term
519
mEh − pEdif in (21) determines which competing effect wins out. For the
520
baseline parameter values chosen here, mEh − pEdif is negative, and thus
521
increasing Tl increases Tp .
522
3.3. Steady-state solutions and the onset of significant shear heating
523
3.3.1. Constant velocity shear zones
524
Using the scaling laws for temperature and fineness at the shear zone
525
core, (17) and (18), steady-state fineness and temperature for constant ve-
526
locity shear zones are determined for a range of strain rates, background
527
temperatures, and two different damage partitioning fractions (Figure 6).
528
The results show the same trends seen in the numerical models, and criti-
529
cally, map out the conditions under which shear heating becomes significant.
530
Shear heating significantly affects shear zone rheology when fineness begins
531
to decrease with increasing strain rate, and thus the onset of significant shear
532
heating is represented by the peak in fineness with increasing strain rate, for
533
a given background temperature, seen in Figure 6. The fineness peak oc-
534
curs at lower values of ε˙0 as Tl decreases, though the influence of Tl is weak.
535
Lowering f shifts the fineness peak to lower ε˙0 , as smaller strain rates are
536
needed to cause significant shear heating when damage is less effective. The
537
fineness peak is found to occur at a constant temperature rise of ≈ 20 − 30
538
K. The onset of significant shear heating can thus be described by defining
539
a critical temperature rise, ∆Tcrit , above which shear heating is important.
540
Setting Tp − Tl = ∆Tcrit in (18) one can solve for the critical strain rate, ε˙crit ,
31
-14 -16
-10 2
-12
1 0 -1 -2 -3 -4 -5 -6
4
-18 900 1000 1100 1200 1300 Background Temperature [K]
B
3
-14 -16
f = 10-9 -10
C
-12
2
1 0 -1 -2 -3 -4
-18 900 1000 1100 1200 1300 Background Temperature [K]
-7 -6 -5 -4 -3 -2 -1 0 1 2 3 log10(Temperature Rise [K])
-14 4
log10(Strain Rate [s-1])
-12
-16
-18 900 1000 1100 1200 1300 Background Temperature [K]
f = 10-9 -10
-14
log10(Strain Rate [s-1])
-12
A
-16
3
D
-18 900 1000 1100 1200 1300 Background Temperature [K]
2
3
4
5
log10(Strain Rate [s-1])
-10
f = 10-4
5
log10(Strain Rate [s-1])
f = 10-4
6
log10(Fineness [m-1])
Figure 6: Steady-state temperature rise (A) and fineness (B) at the shear zone core for f = 10−4 , and for f = 10−9 (C & D). Dashed line marks the onset of significant shear heating, as given by (22), and the dot-dashed line in C & D shows where temperatures reach the nominal melting temperature (i.e. melting would occur above this line). Other model parameters are given in Table 2. The perceptually-uniform color map “roma” is used in this study to prevent visual distortion of the data (Crameri, 2018b,a).
32
Half Width = 104 m
-14
4.5
-12
3.5
3.0
-11
A -9
-8
-7
-6
-5
-4
-3
B
1200 -12
-13
4.0
1300
1100
1000
900
-13
log10(Half Width [m])
-15
-11
5.0
Background Temperature [K]
Background Temp. = 1000 K
-9
-8
-7
log10(f)
-6
-5
-4
-3
log10(f)
-15
-14
-13
-12
-11
-10
log10(Strain Rate [s-1]) Figure 7: Critical strain rate marking the onset of significant shear heating, ε˙crit , as a function of f and shear zone half width, L/2 (A), and f and background temperature, Tl (B).
541
marking the fineness peak as ε˙crit =
542
2∆Tcrit k (1 − f )L2
p+m 2p
4f γhn
m 2p
1 2
(2Bdif 0 ) exp
mEh − pEdif . (22) 2pR(Tl + ∆Tcrit )
Equation (22) is plotted in Figure 6 assuming ∆Tcrit = 20 K.
543
Tradeoffs between how damage and shear zone parameters control whether
544
significant shear heating occurs can be assessed from (22). For example,
545
ε˙crit ∼ L−(p+m)/p , demonstrating that a narrow shear zone requires a higher
546
driving strain rate in order for shear heating to be important, due the influ-
547
ence of shear zone width on heat conduction (e.g. see Figure 3F). Likewise
548
ε˙crit ∼ (f /hn )m/(2p) , as the higher the damage to healing ratio, the more shear 33
549
heating is suppressed by grain size reduction, and thus the larger strain rate
550
must be for significant shear heating.
551
Figure 7A shows ε˙crit for variable f and L. For the baseline value of f =
552
10−4 , the critical strain rate for shear heating to become significant is larger
553
than typical tectonic strain rates. For a 5 cm yr−1 boundary velocity, tectonic
554
strain rates range from 1.5×10−12 −1.5×10−14 s−1 for L/2 = 1−100 km. Thus
555
tectonic strain rates are insufficient to induce significant shear heating when
556
f = 10−4 in the purely viscous models presented here. With the lower bound
557
of f = 10−9 , shear zones of 1 km half width or smaller still require strain
558
rates higher than typical tectonic strain rates for significant shear heating,
559
though wider shear zones can experience shear heating at tectonic strain
560
rates; however, even in this case, the temperature rises are typically low, on
561
the order of 10 K. Only shear zones on the order of 100 km wide, as wide as
562
an entire plate boundary (e.g. Bird, 2003), would experience shear heating
563
of 100 K or more as a result of lower lithosphere ductile flow based on these
564
models. The background temperature also influences the onset of significant
565
shear heating, though with a lesser influence than shear zone width. As
566
Tl decreases ε˙crit also decreases, and thus shear zones are more likely to
567
experience appreciable shear heating at lower background temperatures, as
568
also seen in Figures 3 & 6.
569
The constant velocity shear zone results clearly demonstrate that grain
570
size reduction dampens shear heating, by lowering the effective shear zone
571
viscosity, and hence lowering the deformational work rate. As a result, an
572
important implication is that ignoring grain size reduction in shear zone
573
models leads one to significantly over predict the magnitude of shear heating.
34
-10 -12 -14 -16
A
Damage, f = 10-9 3
Melt
2
1 0 -1 -2 -3 -4 -5 -6
No Damage 3
B
Melt
-10 -12
2
2
1
-14
1
0
0
-1
-2
-3 -4
-18 900 1000 1100 1200 1300 Background Temperature [K]
C
-5
900 1000 1100 1200 1300 Background Temperature [K]
-1 -2 -3 -4
-16
-18 900 1000 1100 1200 1300 Background Temperature [K]
log10(Strain Rate [s-1])
log10(Strain Rate [s-1])
Damage, f = 10-4
-7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 log10(Temperature Rise [K])
Figure 8: Temperature rise caused by shear heating when grain size evolution is considered with f = 10−4 (A) and with f = 10−9 (B), and for a constant grain size of 1 cm (i.e. with no damage considered) (C). Shear zone half width is L/2 = 10 km, and other parameters are given in Table 2. Above the dot-dashed line melting would occur in the shear zone.
574
Figure 8 compares the steady-state temperature rise for shear zones where
575
grain size reduction is taken into account, and those where grain size is
576
assumed to be constant. Shear heating causes a temperature rise of 100 K at
577
strain rates of ε˙ ∼ 10−12 −10−14 s−1 for background temperatures of 1300-900
578
K, respectively, when grain size is fixed, while with grain size reduction shear
579
heating at such conditions is suppressed. A temperature rise of 100 K is not
580
seen until strain rates reach ∼ 10−11 s−1 with f = 10−4 or ∼ 10−12 − 10−13
581
with f = 10−9 , due to the effects of grain size reduction.
582
3.3.2. Constant stress shear zones
583
For constant stress shear zones, steady-state fineness and temperature
584
rise are calculated from (20) and (21) and shown in Figure 9. The region
585
where thermal runaway occurs, and no steady-state solution exists, is also
586
shown. As in the constant velocity case, the trends in Tp and As with Tl , 35
f = 10-4
A
7
8
B
Thermal Runaway
Thermal Runaway
7
6
6 2 0
5
5
-2 -4 -6 -8
5
log10(Stress [Pa])
log10(Stress [Pa])
8
f = 10-4
4
4 4 900 1000 1100 1200 1300 900 1000 1100 1200 1300 Background Temperature [K] Background Temperature [K]
log10(Stress [Pa])
8
Thermal Runaway
7 6 5
f = 10-9 Thermal Runaway
C
2 0
8
D
7
4
-2 -4 -6 -8 -10 -12 -14
6
3
5
2
log10(Stress [Pa])
f = 10-9
4 4 900 1000 1100 1200 1300 900 1000 1100 1200 1300 Background Temperature [K] Background Temperature [K]
-14-12-10 -8 -6 -4 -2 0 2 log10(Temperature Rise [K])
1
2
3
4
5
6
log10(Fineness [m-1])
Figure 9: Steady-state temperature rise (A) and fineness (B) at the shear zone core for constant stress boundary conditions and f = 10−4 . Where thermal runaway occurs (dark red region) no steady-state solution exists. Steady-state temperature rise (C) and fineness (D) are also shown for f = 10−9 . Dashed line denotes where thermal runaway occurs. Shear zone half width is L/2 = 10 km, and other parameters are given in Table 2.
36
587
τl , and f are the same as previously discussed. A salient feature of Figure 9
588
is that lowering f significantly increases the stress required to drive thermal
589
runaway. The critical stress needed to induce thermal runaway can be found
590
as follows. First, (21) is rewritten as (1 − f )L2 F 2 = Tl + 4k
2f γhn
m p−m
τl2 Bdif 0
p p−m
exp
mEh − pEdif RTp (p − m)
− Tp , (23)
591
where F2 = 0 when a steady-state solution for Tp is found. When steady-
592
state solutions exist, there are two solutions, one at low Tp that matches the
593
numerical model results, and one at higher Tp . Between these two steady-
594
state solutions lies a minimum in F2 , that is less than 0. With increasing
595
τl , the minimum of F2 increases, and eventually reaches 0; at this point
596
both steady-state solutions occur at the same Tp . For τl any higher than
597
this critical point, no steady-state solutions exist and the system goes into
598
thermal runaway. Assume that Tpcrit is a solution to (21) (i.e. it satisfies
599
F2 = 0) at the critical stress for the onset of thermal runaway, τcrit . Tpcrit
600
will also be a minimum of F2 , and thus τcrit can be found from dF2 /dTp = 0,
601
evaluated at Tp = Tpcrit :
τcrit =
−1 Bdif20
2 k 4R(p − m)Tpcrit (pEdif − mEh )(1 − f )L2
p−m 2p
2f hn γ
m − 2p
exp
pEdif − mEh . 2pRTpcrit (24)
602
Tpcrit is then found by plugging (24) into (21), at Tp = Tpcrit and τl = τcrit ,
603
and solving for Tpcrit : Tpcrit =
pEdif − mEh 2R(p − m)
4RTl (p − m) 1− 1− pEdif − mEh 37
21 ! .
(25)
604
Together, (24) and (25) give τcrit as a function of Tl , L, and parameters
605
for diffusion creep rheology and grain size evolution. Equations (24) and
606
(25) demonstrate how grain size evolution alters the critical stress threshold
607
where runaway thermal instability occurs. As previously explained, decreas-
608
ing grain size softens the shear zone, thereby enhancing the rate of deforma-
609
tional work and causing thermal runaway to occur at lower stresses than it
610
would without grain size reduction (see also Thielmann et al., 2015). Equa-
611
tion (24) captures this effect, showing that τcrit ∼ (f /hn )−m/(2p) , such that
612
with increasing damage to healing ratio thermal runaway occurs at a lower
613
driving stress. Likewise increasing Bdif 0 , which would lower the bulk shear
614
zone viscosity, also lowers τcrit . The critical stress is also inversely propor-
615
tional to L, as a wider shear zone allows for more heating and thus can more
616
easily experience thermal runaway.
617
3.4. Implications for shear zone strength and the operation of plate tectonics
618
A critical motivating question of this work is whether coupling between
619
shear heating and grain size reduction impedes weak shear zone formation
620
when both effects are significant. Figure 10 shows the shear zone core viscos-
621
ity and stress (assuming constant velocity boundary conditions) as a function
622
of strain rate, for a shear zone with only grain size reduction, one with only
623
shear heating, and one where shear heating and grain size reduction are com-
624
bined. At low strain rates, grain size reduction produces lower stresses and
625
viscosities, i.e. a weaker shear zone, than shear heating. Likewise at high
626
strain rates the opposite is true, and shear heating leads to greater weakening
627
than grain size reduction. With combined shear heating and grain size reduc-
628
tion, shear zone strength is controlled by whichever mechanism provides the
38
100 10-2 10-4 10-6
Stress [Pa]
109
A C
108 107 106 105 104 10-18 10-16 10-14 10-12 10-10 Strain Rate [s-1]
Fineness [m-1]
102
Strain Rate [s-1] 10-18 10-16 10-14 10-12 10-10 107
B
106 105 104 103
Viscosity [Pa s]
Temperature Rise [K]
Strain Rate [s-1] 10-18 10-16 10-14 10-12 10-10 104
1026 D 1024 22 10 1020 1018 1016 1014 1012 10-18 10-16 10-14 10-12 10-10 Strain Rate [s-1]
Shear heating and damage scaling law Damage only scaling law Shear heating only scaling law Shear heating and damage model Damage only model Shear heating only model
Figure 10: Temperature rise caused by shear heating (A), shear zone core fineness (B), stress (C), and viscosity (D) as a function of strain rate. Numerical shear zone model results are shown as symbols, and scaling laws as lines. A case where only shear heating is considered (diamonds and dotted line), where only grain size reduction (damage) is considered (squares and dashed line), and one where damage and shear heating are combined (cirles and solid line). Boundary conditions are constant velocity, f = 10−4 , L/2 = 10 km and Tl = 1000 K.
39
629
greatest degree of weakening, and hence the least resistance to deformation.
630
At low strain rates a shear zone where shear heating and grain size reduction
631
is combined has identical strength to a shear zone with only grain size reduc-
632
tion, and at high strain rates it has identical strength to one where only shear
633
heating is considered. At the transition point where shear heating becomes
634
significant, and thus both grain size evolution and shear heating materially
635
impact the rheology, the combined model is significantly weaker than either
636
end-member case, as grain size reduction and shear heating both contribute
637
to weakening the rock. Thus, rather than impeding each other and limiting
638
weakening, combined shear heating and grain size reduction allow weak plate
639
boundaries to form over a wider range of conditions than either mechanism
640
acting in isolation. Specifically at low strain rates grain size reduction leads
641
to significant weakening, well beyond that which would be provided by shear
642
heating, while at high strain rates shear heating induces significantly more
643
weakening than grain size reduction alone can provide.
644
4. Results: Time evolution
645
4.1. Constant stress shear zones
646
In addition to the final steady-state a shear zone will reach, I also con-
647
strain the time to steady-state to assess whether the timescales for shear zone
648
temperature and fineness evolution are geologically reasonable. I first con-
649
sider constant stress boundary conditions (§3.1.2), where the timescale for
650
reaching steady-state, or thermal runaway, can be well described by a sim-
651
ple scaling law. Time evolution results for a range of conditions are shown
652
in Figure 11; models that experience thermal runaway are stopped as shear
653
zone temperature reaches the melting temperature. 40
Temperature Rise [K]
102 100 10-2 10-4 10-6 10-8 10-10 10-12
C
102 100 10-2 10-4 10-6 10-8 10-10 10-12
E
Time [Myrs] 10-5 100 105 1010 1015 106 Fineness [m-1]
Temperature Rise [K]
A
B Tl = 900 K Tl = 1000 K Tl = 1100 K Tl = 1200 K Numerical Model Scaling Law
105 104 103 102
Fineness [m-1]
106
D τl = 105 Pa τl = 106 Pa τl = 107 Pa τl = 108 Pa Numerical Model Scaling Law
105 104 103 102
102 G 100 10-2 10-4 10-6 10-8 10-10 10-12 10-5 100 105 1010 1015 Time [Myrs]
Fineness [m-1]
106
F
105
L/2 = 1 km L/2 = 10 km L/2 = 100 km Numerical Model Scaling Law
104 103 102 106
Fineness [m-1]
Temperature Rise [K]
102 100 10-2 10-4 10-6 10-8 10-10 10-12
Temperature Rise [K]
Time [Myrs] 10-5 100 105 1010 1015
H
105
f = 10-4 f = 10-6 f = 10-9 Numerical Model Scaling Law
104 103 102 10-5 100 105 1010 1015 Time [Myrs]
Figure 11: Time evolution of temperature rise at the shear zone center, Tp − Tl , and fineness at the shear zone center. In panels A & B, Tl is varied for a constant τl = 107 Pa, f = 10−4 , and L/2 = 10 km. In panels C & D, τl is varied for a constant Tl = 1000 K, f = 10−4 , and L/2 = 10 km. Panels E & F vary L/2 with a constant τl = 107 Pa, f = 10−4 , and Tl = 1000 K, and G & H vary f for constant τl = 107 Pa, L/2 = 10 km, and Tl = 1000 K. Blue lines give the time to reach steady-state from the scaling laws developed in §4.1.
41
654
With a constant stress, fineness and temperature show a rapid rise to
655
their steady-state values, or into a thermal runaway when no steady-state
656
exists, after a long period of slow change.
657
and fineness occur at the same time, as grain size reduction enhances shear
658
heating, and vice-versa, for constant stress boundary conditions. The timing
659
of this rapid rise can be understood from the equation for fineness (9) ne-
660
glecting the grain growth term, as during the rapid rise in fineness grain size
661
reduction is the dominant process:
The rapid rises in temperature
dA fΨ = . dt γ
(26)
662
As the initial fineness is small, deformation typically begins in the dislocation
663
creep regime. I first determine the time to reach the field boundary fineness
664
by integrating (26) from A0 to Af ield for a pure dislocation creep rheology: tf ield
(Af ield − A0 )γ = exp 2f τln+1 Bdis0
Edis . RTl
(27)
665
If A0 > Af ield , meaning deformation begins in the diffusion creep regime,
666
then tf ield is set to 0. In the diffusion creep regime, (26) can be integrated
667
to give A=
668
1 1−m (1 − m)2f τl2 Bdif 0 t Edif exp − + A1−m , f ield γ RTl
(28)
which has a singularity at γA1−m f ield ts = exp (m − 1)2f τl2 Bdif 0 42
Edif . RTl
(29)
669
Runaway grain size reduction occurs as this singularity is approached, and
670
continues until grain growth becomes important and steady-state is reached,
671
or until the shear zone melts during thermal runaway. The total time for
672
rapid grain size reduction to begin is then given by ttot = tf ield + ts , which
673
matches the numerical model results well (Figure 11). Specifically, the scal-
674
ing laws accurately capture the increase in time needed for the onset of rapid
675
grain size reduction (fineness growth) caused by decreasing Tl , τl , or f . Low-
676
ering Tl , τl , or f all slow grain size reduction by either lowering the rate of
677
deformational work (Tl and τl ), or the amount of deformational work that
678
drives grain size reduction (f ). Changing L has no discernible effect on the
679
time to rapid grain size reduction, because deformational work rate is not
680
dependent on L when the driving stress is fixed.
681
Using (27) and (29), the timescale for the fineness rise is calculated for
682
a range of driving stresses and background temperatures (Figure 12). Only
683
at high stress, high temperature conditions (from above 1000 K at τl ≈ 109
684
Pa to above 1300 K at τl ≈ 107 Pa) can significant grain size reduction and
685
heating occur on geologically reasonable timescales of ∼ 100 Myrs or less,
686
for the constant stress boundary condition models employed here. Lower
687
stress and temperature conditions produce improbably long timescales for
688
significant grain size reduction to occur, or for thermal runaway to develop
689
when there is no steady-state solution. Thus, although thermal runaway is
690
predicted to occur at stresses as low as 1 MPa in this study, it may take
691
too long for the instability to develop to be seen in natural settings at these
692
conditions. However, note that additional deformation mechanisms, such as
693
geometric softening due to viscous anisotropy (e.g. Tommasi et al., 2009;
43
1200 2 2
4
1100
4
6
1000
1000 6
10
12
5
8
900
1100
6 7 8 log10(Stress [Pa])
9
-2 0 2 4 6 8 10 12 14 log10(Time [Myrs])
-20 -18 -16 -14 -12 -10 log10(Strainrate [s-1])
0
2
4
6
8
900
10 12
log10(Time [Myrs])
Figure 12: Time for rapid grain size reduction to begin as a function of stress, τl , and background temperature, Tl , as calculated from the scaling laws in §4.1. Shear zone width of L/2 = 10 km and f = 10−4 are assumed.
44
Background Temperature [K]
0
1200
8
Background Temperature [K]
1300
-2
1300
694
Hansen et al., 2012b), could be important during the early stages of shear
695
zone formation, and may speed up the development of fine grained mylonite
696
zones, or the onset of thermal runaway.
697
Dynamic recrystallization experiments constrain the damage to healing
698
ratio, but not f uniquely, so there is significant uncertainty in the exact
699
value of f . As the timing of rapid grain size reduction scales with 1/f , this
700
uncertainty maps into the timescale analysis presented in this section. For
701
example, using f = 0.1 would result in significant grain size reduction devel-
702
oping 1000 times faster than shown in Figure 12, and allow weak shear zones
703
to form in 100 Myrs or less for a broader range of stresses and background
704
temperatures. The influence of f on the timing of rapid grain size reduction
705
likely plays a major role in why Mulyukova & Bercovici (2018) find shorter
706
timescales for significant grain size reduction to occur (generally on the order
707
of 10-100 Myrs) then what is presented in Figure 12; they use f = 10−2 −10−1
708
at lithospheric temperatures and stresses on the order of 108 . Note also that
709
Mulyukova & Bercovici (2018) involves a driving stress that increases with
710
time, as they model passive margin evolution, instead of the constant stress
711
used in this study.
712
4.2. Constant velocity shear zones
713
With a constant velocity boundary condition, timescales to reach steady
714
state are similar to the constant stress case. However, unlike the constant
715
stress case, there is no grain size reduction singularity, and grain size re-
716
duction takes place more gradually from the initial grain size to the final,
717
steady-state grain size (Figure 13). As a result, a simple scaling law for the
718
time when steady-state grain size is reached is more difficult to derive, and 45
100 10-2 10-4
106 Fineness [m-1]
A
Time [Myrs] 10-5 100 105 1010 1015
102
C
100 10-2 10-4 10-6 10-8
E
102 100 10-2 10-4
D ε˙ 0 = 10-17 s-1 ε˙ 0 = 10-15 s-1 ε˙ 0 = 10-13 s-1 ε˙ 0 = 10-11 s-1
105 104 103
F
105
L/2 = 1 km L/2 = 10 km L/2 = 100 km
104 103 102
G
100 10-2 10-4 10-6
103
106
10-6 102
104
102 Fineness [m-1]
104
Tl = 900 K Tl = 1000 K Tl = 1100 K Tl = 1200 K
105
106 Fineness [m-1]
104
B
102
10-6
10-5 100 105 1010 1015 Time [Myrs]
106 Fineness [m-1]
Temperature Rise [K] Temperature Rise [K] Temperature Rise [K] Temperature Rise [K]
102
Time [Myrs] 10-5 100 105 1010 1015
H
105
f = 10-4 f = 10-6 f = 10-9
104 103 102
10-5 100 105 1010 1015 Time [Myrs]
Figure 13: Time evolution of temperature rise at the shear zone center, Tp − Tl , and fineness at the shear zone center. In panels A & B, Tl is varied for a constant ε˙0 = 10−15 s−1 , f = 10−4 , and L/2 = 10 km. In panels C & D, ε˙0 is varied for a constant Tl = 1000 K, f = 10−4 , and L/2 = 10 km. Panels E & F vary L/2 with a constant ε˙0 = 10−15 s−1 , f = 10−4 , and Tl = 1000 K, and G & H vary f for constant ε˙0 = 10−15 s−1 , L/2 = 10 km, and Tl = 1000 K.
46
719
hence not attempted here. However, the numerical results give guidance on
720
the typical timescales for grain size reduction and shear heating. The models
721
show that there are two important time scales: one where significant grain
722
size reduction occurs and there is a peak in fineness, and a longer timescale
723
to reach the final steady-state. Unlike the constant stress case, background
724
temperature has only a minor influence on the time for significant grain size
725
reduction to occur (Figure 13B), because of two competing effects: Lowering
726
the background temperature increases the deformational work rate and hence
727
rate of grain size reduction, but at the same time smaller grain sizes must be
728
reached before grain growth becomes important, at least when shear heating
729
is minor. The time to reach the final steady-state does increase modestly with
730
decreasing background temperature. Increasing the imposed strain rate, ε˙0 ,
731
has a strong effect on the time for significant grain size reduction to develop,
732
and, intuitively, lowers this timescale (Figure 13D). However, the time to
733
reach steady-state is not strongly sensitive to ε˙0 . Shear zone width has very
734
little influence on the time to reach steady-state, while damage fraction, f ,
735
does; increasing f leads to faster grain size reduction and thus shortens the
736
time for both significant grain size reduction to occur and for steady-state to
737
be reached (Figures 13F & H).
738
Overall, timescales are unreasonably long, greater than 1 Gyrs, for most
739
conditions modeled. Only at high imposed strain rates, ε˙0 ∼ 10−13 − 10−11
740
s−1 , can shear zones form on 1 − 10 Myr timescales. Although such strain
741
rates are not unreasonable for 1-10 km wide shear zones, driving a shear
742
zone at ε˙0 ∼ 10−13 − 10−11 s−1 before significant grain size reduction and/or
743
shear heating has developed requires high stresses. Specifically stresses of
47
744
τ ∼ 1 − 10 GPa are required to drive a shear zone at strain rates of ε˙0 ∼
745
10−13 − 10−11 s−1 , with a background temperature of 1000 K; such stresses
746
are higher than typical tectonic stresses (e.g. Behr & Platt, 2011). Thus,
747
as discussed in §4.1, additional weakening mechanisms that are active early
748
in shear zone formation may be necessary for mylonite zones to form on
749
geologically reasonable timescales. Note that localization might also hasten
750
shear zone development, by increasing strain rates in the localized zone.
751
Therefore more sophisticated damage theories that allow for localization, as
752
discussed in §5.3, may also be important for shear zone formation.
753
5. Discussion
754
5.1. Implications for observations of natural shear zones
755
Overall, the results demonstrate that in stably flowing shear zones that
756
can be described with constant velocity boundary conditions, like plate bound-
757
aries, and where elasticity is negligible, shear heating is diminshed by grain
758
size reduction for most of the conditions explored. Grain size reduction
759
weakens the shear zone and thus lowers the deformational work that drives
760
shear heating. At the mid- to lower-lithosphere temperatures explored in this
761
study either weakly effective damage or very high strain rates are needed for
762
shear heating of ∼ 10 K, which is sufficient to have an impact on the shear
763
zone grain size and rheology. Shear heating of 100 K or more is only seen
764
for an even narrower range of conditions. Mont´esi (2013) argues that shear
765
heating of greater than 100 K is needed for significant localization, so other
766
mechanisms, such as fabric development or feedbacks between grain size re-
767
duction and phase mixing (see §5.3), are likely to be more important for 48
768
shear localization in lower lithosphere shear zones. However, with a constant
769
stress boundary condition the results indicate that grain size reduction sig-
770
nificantly enhances the likelihood of thermal runaway instabilities, consistent
771
with previous studies (see §5.2), by lowering the stress threshold necessary
772
for triggering thermal runaway. Thus transient pulses of significant shear
773
heating at shallower depths, where elasticity is important, may be common,
774
so long as instability can develop on geologically reasonable timescales.
775
Whether significant shear heating is observed in natural ductile shear
776
zones, however, has been debated. Metamorphic assemblages in exhumed
777
fault zones and shear zones have been interpreted as showing shear heating
778
in some studies (Graham & England, 1976; England & Molnar, 1993; Leloup
779
& Kienast, 1993; Camacho et al., 2001), though alternative explanations
780
that do not require shear heating have also been proposed for many of these
781
localities (Leloup et al., 2001; Gilley et al., 2003; Kidder & Ducea, 2006;
782
Chapman et al., 2011; Kohn, 2014). If shear heating in the deeper lithosphere
783
is generally muted, then suppression of shear heating via grain size reduction
784
provides a possible explanation. This study also predicts that for constant
785
strain rate shear zones when shear heating is non-negligible, the temperature
786
distribution in the shear zone leads to grain growth in the shear zone core.
787
As a result, grain sizes are largest in the shear zone center and become
788
progressively smaller towards the shear zone edge, where temperatures are
789
cooler. Such a grain size distribution within a shear zone has not been
790
observed in nature to my knowledge, but observing such a trend would be
791
an indicator that non-negligible levels of shear heating took place.
49
792
5.2. Comparison to previous studies
793
This work reaffirms the results of Thielmann et al. (2015) & Thielmann
794
(2017), that grain size reduction lowers the critical stress needed to trigger
795
thermal runaway. In this study, grain size reduction is more extensive than
796
in Thielmann et al. (2015) & Thielmann (2017), as the pinned state formu-
797
lation for grain damage both suppresses grain growth and allows grain size
798
reduction to continue in the diffusion creep regime. As a result, I find the
799
stress threshold for thermal runaway can be lowered even further than what
800
is shown in Thielmann et al. (2015) & Thielmann (2017). However, the long
801
time scales needed for instability to develop mean high stresses of 100 MPa
802
or more may still be necessary for thermal runaway to be seen in natural
803
shear zones, as explained in §4.1.
804
Hansen et al. (2012a) compared rock deformation experiments performed
805
at constant strain rate and constant stress conditions, and found that the
806
constant stress conditions enhance shear localization and lead to unstable
807
deformation, where strain accelerates during deformation. This is similar
808
to what is found here, where constant velocity boundary conditions due not
809
allow for a runaway instability, while constant stress conditions do. However,
810
there are important differences to note as well. In Hansen et al. (2012a),
811
localization occurs due to the combined effects of grain size reduction and
812
fabric development, while here fabric development is not modeled (viscosity
813
is assumed to be isotropic) and grain size reduction alone does not lead to
814
localization, as explained in §3.1.1. Localization seen in Hansen et al. (2012a)
815
could be a transient effect that would be erased at higher strain, or a result
816
of effects not included in the grain size formulation used in this study, such
50
817
as fabric development or the influence grain size itself has on the effectiveness
818
of grain size reduction (discussed further in §5.3).
819
5.3. Model limitations
820
As the goal of this study is to constrain the first order effects of coupled
821
shear heating and grain size reduction on ductile shear zone strength, a sim-
822
ple model is used that can be easily understood and analyzed. Namely an
823
effectively one-dimensional shear zone model with a constant temperature at
824
the shear zone edge, which allows steady-state solutions to exist, is employed.
825
Other possibilities are a no heat flux boundary at the shear zone edge or far
826
from the shear zone edge (i.e. as x goes to infinity), neither of which allow
827
for steady-state solutions. Ultimately vertical heat transport to the cold sur-
828
face and to the base of the lithosphere will limit the amount of shear heating
829
that can take place, and allow a steady-state thermal profile to be reached,
830
in a way that is analogous to the simple one-dimensional model used in this
831
study. Nevertheless, it is likely that the models shown here underestimate
832
the peak temperature reached in the shear zone core and skew the onset of
833
significant shear heating to higher driving stress or strain rate, as a result
834
of the constant temperature condition applied to the shear zone edge. How-
835
ever, the overall results, that shear heating is insignificant and grain size
836
reduction is extensive at low stress or strain rate, and shear heating is sig-
837
nificant while grain size reduction is hampered at higher τl or ε˙0 , is unlikely
838
to be significantly affected by the temperature boundary condition. Likewise
839
the result illustrated in §3.4, that combined grain size reduction and shear
840
heating allows weak shear zones to form over a wider range of conditions
841
than either mechanism acting in isolation, should also hold true regardless of 51
842
the temperature boundary condition, as it is caused by the tendency of rock
843
deformation to occur via whichever mechanism allows for the highest strain
844
rate.
845
The effectively one-dimensional, simple shear model setup used also means
846
any possible vertical flow in the shear zone is not captured. If vertical flow
847
is significant it would change the deformational work rate, and hence grain
848
size and temperature evolution, in the shear zone, but ultimately this issue is
849
beyond the scope of this study as it involves a different shear zone geometry.
850
The rheology assumed is also purely viscous, when in reality elasticity can
851
be important at low temperatures. In particular, elasticity can drive thermal
852
runaway even with constant velocity boundary conditions, by releasing stress
853
that builds up during shearing. The models presented here therefore cannot
854
capture these thermal runaway scenarios, which would result in lower grain
855
sizes and higher temperatures in shear zones than shown in this study. How-
856
ever, as explained in §2.2, elasticity is negligible for most of the conditions
857
explored in this paper.
858
The models also use a constant damage partitioning fraction, f . While the
859
choice of Eh in this study eliminates any possible temperature dependence
860
for f , at least as inferred from dynamic recrystallization experiments (see
861
Appendix A), f may also be a function of grain size (Bercovici & Ricard,
862
2016; Mulyukova & Bercovici, 2017). Damage is most efficient when pyroxene
863
and olivine grains are well mixed, allowing Zener pinning to both suppress
864
grain growth and enhance grain size reduction (Bercovici & Ricard, 2012).
865
As mixing is thought to be more efficient at smaller grain sizes than larger
866
ones (Bercovici & Skemer, 2017; Cross & Skemer, 2017), the pinned state is
52
867
not reached until grains are reduced to small enough sizes for the phases to
868
easily mix. Thus there are two limits to the effective value of f ; a smaller
869
value of f characteristic of rock that is not in the pinned state when grain
870
size is large, and a larger value, characteristic of the pinned state, when grain
871
sizes are small. In this study, the pinned state is assumed to prevail at all
872
conditions, meaning that the value of f used is the value for the pinned state.
873
Allowing f to vary with grain size would allow shear localization, forming
874
a fine-grained region at the shear zone core where damage is operating in
875
the pinned state, with the rest of the shear zone seeing coarser grains and
876
low strain rates (Bercovici & Ricard, 2016). How shear heating would be
877
affected with a grain size dependent f is not known without running full
878
models, but it is likely that when significant grain size reduction and shear
879
localization occurs, shear heating would be suppressed for constant velocity
880
boundary conditions, as in the results shown in this study. The conditions
881
for the onset of siginficant shear heating may be affected as well.
882
Finally, I note that additional weakening mechanisms that are not in-
883
cluded in this study could be important. In particular, weakening beyond
884
what is provided by shear heating and damage alone may be important dur-
885
ing the initial stages of shear zone formation, and thus help speed up the
886
development of significant grain size reduction or heating. As deformation
887
likely starts outside the diffusion creep regime, viscous anisotropy that re-
888
sults from LPO formation may be one such mechanism. Furthermore, creep
889
mechanisms such as low temperature plasticity and dislocation accommo-
890
dated grain boundary sliding, neglected in this study, may have an impor-
891
tant effect as well (e.g. Hansen et al., 2012c; Thielmann, 2017). Ultimately,
53
892
understanding the deformation mechanisms prevalent early in shear zone for-
893
mation may be critical to understanding how weak shear zones develop at
894
all.
895
6. Conclusions
896
Models of ductile lithospheric shear zones constrain the co-evolution of
897
grain size and temperature for both constant velocity and constant stress
898
boundary conditions. The constant velocity models demonstrate that at low
899
strain rates, shear heating is negligible and grain size shrinks with increasing
900
strain rate. At high strain rates shear heating becomes significant, and high
901
temperatures cause grain size to increase with increasing strain rate due to
902
rapid grain growth rates. For all but very high strain rates, where high
903
temperatures lead to large grain sizes and thus push deformation into the
904
dislocation creep regime, grain size reduction diminishes the amount of shear
905
heating, by lowering the shear zone viscosity and hence the deformational
906
work available to drive heating. As a result, ignoring grain size evolution
907
leads one to significantly overestimate the amount of shear heating during
908
simple shear deformation at geologically relevant conditions. For example,
909
without considering grain size evolution, hundreds of degrees of shear heating
910
can occur at strain rates of 10−14 s−1 at mid- to lower-lithosphere background
911
temperatures, while strain rates must be at least 10-100 times larger for
912
the same level of heating when grain size evolution is considered. Thus
913
long lived, stable shear zones that can be described by constant velocity
914
boundary conditions, are less likely to experience significant shear heating
915
than previously thought, at least for the range of background temperatures 54
916
considered here.
917
With a constant stress boundary condition, grain size reduction is sig-
918
nificant and shear heating minor at low stresses, where stable steady-states
919
can be reached. At higher stresses thermal runaway instability occurs and
920
the shear zone experiences very high temperatures. Consistent with previ-
921
ous studies, the results indicate that grain size reduction lowers the stress
922
threshold where thermal runaway takes place. In fact, the stress needed to
923
trigger thermal runaway is found to be even lower than previous studies that
924
also looked at the combined effects of shear heating and grain size reduc-
925
tion, because grain boundary pinning by secondary phases, which this study
926
incorporates, enhances grain size reduction. However, the timescale needed
927
for thermal runaway to occur in the mid-lithosphere at stresses lower than
928
100 MPa is found to be greater than 100 Myrs, so high stresses may still be
929
necessary for thermal runaway to take place in practice.
930
Finally, rather than acting against each other in terms of providing rhe-
931
ological weakening, combining shear heating and grain size reduction allows
932
weak shear zones to form over a wider range of conditions than would be
933
possible from either mechanism acting in isolation. At low strain rates shear
934
heating is negligible and thus an ineffective weakening mechanism, but grain
935
size reduction is extensive. At high strain rates shear heating is significant
936
while grain size reduction is muted by the high temperatures, and shear heat-
937
ing provides substantial weakening. Thus whichever mechanism allows for
938
the least resistance to flow will control the rheology. Even when both shear
939
heating and grain size reduction are non-negligible, their combined effects
940
provide more weakening than either mechanism acting alone.
55
941
7. Acknowledgements
942
I thank Elvira Mulyukova and Dave Bercovici for discussions on experi-
943
mental constraints on damage theory, and Elvira Mulyukova for sharing some
944
of her codes related to this topic. John D. Platt also helped in the initial
945
stages of this work and showed me how to implement the method of lines for
946
shear zone modeling. I also thank Lars Hansen, Laurent Montesi, and two
947
anonymous reviewers for their comments that helped to significantly improve
948
the manuscript.
949
Appendix A. Constraining f ∗ from laboratory experiments
950
The damage partitioning fraction for damage to the grains directly, f ∗ ,
951
is constrained by laboratory data following Mulyukova & Bercovici (2017).
952
Specifically, Mulyukova & Bercovici (2017) fit experimental results to an
953
equation of the form R = Clab τ −klab ,
(A.1)
954
where R is the grain size after dynamic recrystallization and Clab is a function
955
of temperature as (their equation B.6) Clab =
3λ2 γh∗ 4λ3 f ∗ Bdis
31 .
(A.2)
956
Here λ2 ≈ 3.6 is the second moment of the lognormal grain size distribution
957
and λ3 ≈ 17.8 is the third moment. Rearranging for f ∗ ∗
f =
−3 Clab
3λ2 γh∗n 4λ3 Bdis0
56
exp
Edis − Eh . RT
(A.3)
Damage Partitioning Fraction, f *
f * = exp(-9.41(T/Tref)-0.008) f * inferred from experiments
10-3
10-4
10-5 800 1000 1200 1400 1600 1800 2000 Temperature [K]
Figure A.14: Inferred values of f ∗ from laboratory experiments assuming Eh = 430 kJ mol−1 and calculating h∗n from (12) (circles), and fit to the data using (A.4) (solid line).
958
Inferred values of f ∗ can then be determined from the experimental results
959
shown in Mulyukova & Bercovici (2017) using (A.3) with a given Eh . From
960
inspecting the inferred values of f ∗ for various choices of Eh , it was found that
961
Eh ≈ 430 kJ mol−1 results in f ∗ values that are independent of temperature
962
(Figure A.14). Fitting these f ∗ values with a function of the form f ∗ = exp −Cf
T Tref
kf ! (A.4)
963
where Tref = 1000 K is a reference temperature, yields Cf ≈ 9.41 and kf ≈
964
−0.008, and a value of f ∗ ≈ 10−4 .
965
Appendix B. Composite rheology scaling laws
966
Here I show the scaling laws for steady-state temperature and fineness at
967
the shear zone core with a composite rheology. As in the main text, the peak
57
968
temperature is given by Tp = Tl +
969
(1 − f )Ψ(Tp )L2 8k
and fineness by As =
970
(B.1)
f Ψ(Tp ) γh(Tp )
p1 .
(B.2)
With a constant stress, Ψ is easily given as 2 Ψ = 2(Bdis (Tp )τln+1 + Bdif (Tp )Am s τl ),
(B.3)
971
and combining (B.1) - (B.3) with (6), (7), and (10) gives the scaling laws for
972
Tp and As . With a constant boundary velocity, the rheology law, (5), must
973
be inverted to give stress as a function of strain rate in order to write the
974
scaling laws for Tp and As . With n = 3, τ=
975
η (2/3)1/3 Bdif (Tp )Am s − 181/3 Bdis (Tp ) η
(B.4)
where 1
1
2 3 η = [9Bdis (Tp )2 ε˙0 + {3(27Bdis (Tp )4 ε˙20 + 4Bdis (Tp )3 Bdif (Tp )3 A3m s )} ] . (B.5)
976
Using Ψ = 2τ ε˙0 , with τ given by (B.4)-(B.5), (B.1)-(B.2), and (6), (7),
977
and (10), the scaling laws for Tp and As with constant velocity boundary
978
conditions and a composite rheology can be found.
58
979
Behr, W. M., & Platt, J. P. (2011). A naturally constrained stress profile
980
through the middle crust in an extensional terrane. Earth Planet. Sci.
981
Lett., 303 , 181–192.
982
Bercovici, D., & Karato, S. (2003). Theoretical analysis of shear localization
983
in the lithosphere. In S. Karato, & H. Wenk (Eds.), Reviews in Mineralogy
984
and Geochemistry: Plastic Deformation of Minerals and Rocks chapter 13.
985
(pp. 387–420). Washington, DC: Min. Soc. Am. volume 51.
986
Bercovici, D., & Ricard, Y. (2012). Mechanisms for the generation of plate
987
tectonics by two-phase grain-damage and pinning. Phys. Earth Planet.
988
Inter., 202 , 27–55.
989
Bercovici, D., & Ricard, Y. (2013). Generation of plate tectonics with two-
990
phase grain-damage and pinning: Source-sink model and toroidal flow.
991
Earth Planet. Sci. Lett., 365 , 275–288.
992
993
994
995
996
997
Bercovici, D., & Ricard, Y. (2014). Plate tectonics, damage, and inheritance. Nature, 508 , 513–516. Bercovici, D., & Ricard, Y. (2016). Grain-damage hysteresis and plate tectonic states. Phys. Earth Planet. Inter., 253 , 31–47. Bercovici, D., & Skemer, P. (2017). Grain damage, phase mixing and plateboundary formation. Journal of Geodynamics, 108 , 40–55.
998
Bercovici, D., Tackley, P., & Ricard, Y. (2015). The generation of plate
999
tectonics from mantle dynamics. In D. Bercovici, & G. Schubert (chief
1000
editor) (Eds.), Treatise on Geophysics (pp. 271–318). New York: Elsevier
1001
volume 7, Mantle Dynamics. 59
1002
1003
Bird, P. (2003). An updated digital model of plate boundaries. Geochem., Geophys., Geosyst., 4 , 1027.
1004
Braun, J., Chery, J., Poliakov, A., Mainprice, D., Vauchez, A., Tomassi, A.,
1005
& Daignieres, M. (1999). A simple parameterization of strain localization
1006
in the ductile regime due to grain size reduction: A case study for olivine.
1007
J. Geophys. Res., 104 , 25167–25181.
1008
Camacho, A., McDougall, I., Armstrong, R., & Braun, J. (2001). Evidence
1009
for shear heating, Musgrave Block, central Australia. J. Struct. Geol., 23 ,
1010
1007–1013.
1011
Chapman, A., Luffi, P., Saleeby, J., & Petersen, S. (2011). Metamorphic evo-
1012
lution, partial melting and rapid exhumation above an ancient flat slab:
1013
insights from the san emigdio schist, southern california. Journal of Meta-
1014
morphic Geology, 29 , 601–626.
1015
1016
Crameri, F. (2018a). Geodynamic diagnostics, scientific visualisation and StagLab 3.0. Geosci. Model Dev. Discuss., .
1017
Crameri, F. (2018b). Scientific colour maps (Version 3.0.0). Zenodo, .
1018
Cross, A. J., & Skemer, P. (2017). Ultramylonite generation via phase mixing
1019
in high-strain experiments. J. Geophys. Res. Solid Earth, 122 , 1744–1759.
1020
England, P., & Molnar, P. (1993). The interpretation of inverted metamor-
1021
1022
phic isograds using simple physical calculations. Tectonics, 12 , 145–157. Evans, B., Renner, J., & Hirth, G. (2001). A few remarks on the kinetics of
60
1023
static grain growth in rocks. Int. J. Earth Sciences (Geol. Rundsch.), 90 ,
1024
88–103.
1025
Farla, R. J. M., Karato, S.-i., & Cai, Z. (2013). Role of orthopyroxene
1026
in rheological weakening of the lithosphere via dynamic recrystallization.
1027
Proceedings of the National Academy of Science, 110 , 16355–16360.
1028
1029
Fleitout, L., & Froidevaux, C. (1980). Thermal and mechanical evolution of shear zones. J. Struct. Geol., 2 , 159–164.
1030
Foley, B. J., & Becker, T. W. (2009). Generation of plate-like behavior
1031
and mantle heterogeneity from a spherical, visco-plastic convection model.
1032
Geochem., Geophys., Geosyst., 10 . Q08001.
1033
Foley, B. J., & Bercovici, D. (2014).
Scaling laws for convection with
1034
temperature-dependent viscosity and grain-damage. Geophys. J. Int., 199 ,
1035
580–603.
1036
Foley, B. J., Bercovici, D., & Elkins-Tanton, L. T. (2014). Initiation of plate
1037
tectonics from post-magma ocean thermochemical convection. J. Geophys.
1038
Res. Solid Earth, 119 , 8538–8561.
1039
Foley, B. J., Bercovici, D., & Landuyt, W. (2012). The conditions for plate
1040
tectonics on super-earths: Inferences from convection models with damage.
1041
Earth Planet. Sci. Lett., 331332 , 281 – 290.
1042
Foley, B. J., & Driscoll, P. E. (2016). Whole planet coupling between cli-
1043
mate, mantle, and core: Implications for rocky planet evolution. Geochem.,
1044
Geophys., Geosyst., 17 . 61
1045
Gilley, L. D., Harrison, T. M., Leloup, P. H., Ryerson, F. J., Lovera, O. M.,
1046
& Wang, J.-H. (2003). Direct dating of left-lateral deformation along the
1047
Red River shear zone, China and Vietnam. J. Geophys. Res. Solid Earth,
1048
108 , 2127.
1049
Graham, C. M., & England, P. C. (1976). Thermal regimes and regional
1050
metamorphism in the vicinity of overthrust faults: an example of shear
1051
heating and inverted metamorphic zonation from southern California.
1052
Earth Planet. Sci. Lett., 31 , 142–152.
1053
Hanmer, S. (1988). Great Slave Lake Shear Zone, Canadian Shield: recon-
1054
structed vertical profile of a crustal-scale fault zone. Tectonophysics, 149 ,
1055
245–264.
1056
Hansen, L. N., Zimmerman, M. E., Dillman, A. M., & Kohlstedt, D. L.
1057
(2012a). Strain localization in olivine aggregates at high temperature: A
1058
laboratory comparison of constant-strain-rate and constant-stress bound-
1059
ary conditions. Earth Planet. Sci. Lett., 333 , 134–145.
1060
Hansen, L. N., Zimmerman, M. E., & Kohlstedt, D. L. (2012b). Laboratory
1061
measurements of the viscous anisotropy of olivine aggregates. Nature, 492 ,
1062
415–418.
1063
Hansen, L. N., Zimmerman, M. E., & Kohlstedt, D. L. (2012c). The influence
1064
of microstructure on deformation of olivine in the grain-boundary sliding
1065
regime. J. Geophys. Res. Solid Earth, 117 , B09201.
1066
Hiraga, T., Tachibana, C., Ohashi, N., & Sano, S. (2010). Grain growth
62
1067
systematics for forsterite enstatite aggregates: Effect of lithology on grain
1068
size in the upper mantle. Earth Planet. Sci. Lett., 291 , 10 – 20.
1069
Hirschmann, M. M. (2000). Mantle solidus: Experimental constraints and
1070
the effects of peridotite composition. Geochem., Geophys., Geosyst., 1 .
1071
Hirth, G., & Kohlstedt, D. (2003). Rheology of the upper mantle and the
1072
mantle wedge: a view from the experimentalists. In J. Eiler (Ed.), Sub-
1073
duction Factory Mongraph (pp. 83–105). Washington, DC: Am. Geophys.
1074
Union volume 138.
1075
Jin, D., Karato, S., & Obata, M. (1998). Mechanisms of shear localization
1076
in the continental lithosphere: Inference from the deformation microstruc-
1077
tures of peridotites from the Ivrea zone, northwestern Italy. J. Struct.
1078
Geol., 20 , 195–209.
1079
Kameyama, M., Yuen, D., & Fujimoto, H. (1997). The interaction of viscous
1080
heating with grain-size dependent rheology in the formation of localized
1081
slip zones. Geophys. Res. Lett., 24 , 2523–2526.
1082
Kameyama, M., Yuen, D. A., & Karato, S.-I. (1999). Thermal-mechanical
1083
effects of low-temperature plasticity (the Peierls mechanism) on the defor-
1084
mation of a viscoelastic shear zone. Earth Planet. Sci. Lett., 168 , 159–172.
1085
Karato, S. (1989). Grain growth kinetics in olivine aggregates. Tectono-
1086
1087
1088
physics, 168 , 255–273. Karato, S. (2008). Deformation of Earth Materials: An Introduction to the Rheology of the Solid Earth. New York: Cambridge University Press. 63
1089
1090
1091
Karato, S., & Wu, P. (1993). Rheology of the Upper Mantle: A Synthesis. Science, 260 , 771–778. Kaus, B. J. P., & Podladchikov, Y. Y. (2006).
Initiation of localized
1092
shear zones in viscoelastoplastic rocks. J. Geophy. Res. Solid Earth, 111 ,
1093
B04412.
1094
1095
Kelemen, P. B., & Hirth, G. (2007). A periodic shear-heating mechanism for intermediate-depth earthquakes in the mantle. Nature, 446 , 787–790.
1096
Kidder, S., & Ducea, M. N. (2006). High temperatures and inverted meta-
1097
morphism in the schist of Sierra de Salinas, California. Earth Planet. Sci.
1098
Lett., 241 , 422–437.
1099
Kohlstedt, D., Evans, B., & Mackwell, S. (1995). Strength of the lithosphere:
1100
Constraints imposed by laboratory experiments. J. Geophys. Res., 100 ,
1101
17,587–17,602.
1102
1103
Kohn, M. J. (2014). Himalayan Metamorphism and Its Tectonic Implications. Annual Review of Earth and Planetary Sciences, 42 , 381–419.
1104
Landuyt, W., & Bercovici, D. (2009). Variations in planetary convection via
1105
the effect of climate on damage. Earth and Planetary Science Letters, 277 ,
1106
29–37.
1107
Landuyt, W., Bercovici, D., & Ricard, Y. (2008). Plate generation and two-
1108
phase damage theory in a model of mantle convection. Geophys. J. Int.,
1109
174 , 1065–1080.
64
1110
Leloup, P. H., Arnaud, N., Lacassin, R., Kienast, J. R., Harrison, T. M.,
1111
Trong, T. T. P., Replumaz, A., & Tapponnier, P. (2001). New constraints
1112
on the structure, thermochronology, and timing of the Ailao Shan-Red
1113
River shear zone, SE Asia. J. Geophys. Res., 106 , 6683–6732.
1114
Leloup, P. H., & Kienast, J.-R. (1993). High-temperature metamorphism in a
1115
major strike-slip shear zone: the Ailao Shan-Red River, People’s Republic
1116
of China. Earth Planet. Sci. Lett., 118 , 213–234.
1117
Linckens, J., Herwegh, M., & M¨ untener, O. (2015). Small quantity but large
1118
effecthow minor phases control strain localization in upper mantle shear
1119
zones. Tectonophysics, 643 , 26–43.
1120
Linckens, J., Herwegh, M., M¨ untener, O., & Mercolli, I. (2011). Evolution
1121
of a polymineralic mantle shear zone and the role of second phases in the
1122
localization of deformation. J. Geophys. Res. Solid Earth, 116 , B06210.
1123
Mont´esi, L., & Hirth, G. (2003). Grain size evolution and the rheology of
1124
ductile shear zones: From laboratory experiments to postseismic creep.
1125
Earth Planet. Sci. Lett., 211 , 97–110.
1126
Mont´esi, L., & Zuber, M. (2002).
A unified description of localiza-
1127
tion for application to largescale tectonics.
1128
10.1029/2001JB000465.
J. Geophys. Res., 107 ,
1129
Mont´esi, L. G. J. (2007). A constitutive model for layer development in shear
1130
zones near the brittle-ductile transition. Geophys. Res. Lett., 34 , L08307.
1131
Mont´esi, L. G. J. (2013). Fabric development as the key for forming ductile 65
1132
shear zones and enabling plate tectonics. Journal of Structural Geology,
1133
50 , 254–266.
1134
Moresi, L., & Solomatov, V. (1998). Mantle convection with a brittle litho-
1135
sphere: Thoughts on the global tectonic style of the earth and venus.
1136
Geophys. J. Int., 133 , 669–682.
1137
Mulyukova, E., & Bercovici, D. (2017). Formation of lithospheric shear zones:
1138
Effect of temperature on two-phase grain damage. Phys. Earth Planet.
1139
Inter., 270 , 195 – 212.
1140
Mulyukova, E., & Bercovici, D. (2018). Collapse of passive margins by litho-
1141
spheric damage and plunging grain size. Earth Planet. Sci. Lett., 484 ,
1142
341–352.
1143
1144
1145
1146
Platt, J. P. (2015). Rheology of two-phase systems: a microphysical and observational approach. J. Struct. Geol., 77 , 213–227. Poirier, J. (1980). Shear localization and shear instability in materials in the ductile field. J. Struct. Geol., 2 , 135–142.
1147
Regenauer-Lieb, K., & Yuen, D. A. (1998). Rapid conversion of elastic en-
1148
ergy into plastic shear heating during incipient necking of the lithosphere.
1149
Geophys. Res. Lett., 25 , 2737–2740.
1150
Rozel, A., Ricard, Y., & Bercovici, D. (2011). A thermodynamically self-
1151
consistent damage equation for grain size evolution during dynamic re-
1152
crystallization. Geophys. J. Intl., 184 , 719–728.
66
1153
1154
1155
1156
1157
1158
Schiesser, W. E. (2012). The numerical method of lines: Integration of partial differential equations. Elsevier. Schubert, G., & Yuen, D. A. (1978). Shear heating instability in the earth’s upper mantle. Tectonophysics, 50 , 197–205. Skemer, P., Warren, J. M., & Kelemen, P. B. (2010). Microstructural and rheological evolution of a mantle shear zone. J. Petrol., 51 , 43–53.
1159
Tackley, P. (2000a). The quest for self-consistent generation of plate tec-
1160
tonics in mantle convection models. In M. A. Richards, R. Gordon, &
1161
R. van der Hilst (Eds.), History and Dynamics of Global Plate Motions,
1162
Geophys. Monogr. Ser. (pp. 47–72). Washington, DC: Am. Geophys. Union
1163
volume 121.
1164
Tackley, P. (2000b). Self-consistent generation of tectonic plates in time-
1165
dependent, three-dimensional mantle convection simulations, 1. Pseudo-
1166
plastic yielding. Geochem. Geophys. Geosyst. (G3 ), 1 .
1167
Thielmann, M. (2017). Grain size assisted thermal runaway as a nucleation
1168
mechanism for continental mantle earthquakes: Impact of complex rheolo-
1169
gies. Tectonophysics, .
1170
Thielmann, M., & Kaus, B. J. P. (2012). Shear heating induced lithospheric-
1171
scale localization: Does it result in subduction? Earth Planet. Sci. Lett.,
1172
359 , 1–13.
1173
Thielmann, M., Rozel, A., Kaus, B. J. P., & Ricard, Y. (2015). Intermediate-
1174
depth earthquake generation and shear zone formation caused by grain size
1175
reduction and shear heating. Geology, 43 , 791–794. 67
1176
Tommasi, A., Knoll, M., Vauchez, A., Signorelli, J. W., Thoraval, C., &
1177
Log´e, R. (2009). Structural reactivation in plate tectonics controlled by
1178
olivine crystal anisotropy. Nature Geoscience, 2 , 423–427.
1179
1180
1181
1182
Turcotte, D., & Schubert, G. (2002). Geodynamics. (2nd ed.). New York: Cambridge Univ. Press. van Heck, H., & Tackley, P. (2008). Planforms of self-consistently generated plates in 3d spherical geometry. Geophys. Res. Lett., 35 , L19312.
1183
Warren, J. M., & Hirth, G. (2006). Grain size sensitive deformation mech-
1184
anisms in naturally deformed peridotites. Earth Planet. Sci. Letts., 248 ,
1185
438 – 450.
1186
1187
White, S., Burrows, S., Carreras, J., Shaw, N., & Humphreys, F. (1980). On mylonites in ductile shear zones. J. Struct. Geol., 2 , 175–187.
1188
Yuen, D. A., Fleitout, L., Schubert, G., & Froidevaux, C. (1978). Shear defor-
1189
mation zones along major transform faults and subducting slabs. Geophys.
1190
J. Int., 54 , 93–119.
1191
Yuen, D. A., & Schubert, G. (1979). On the stability of frictionally heated
1192
shear flows in the asthenosphere. Geophysical Journal International , 57 ,
1193
189–207.
68
S0031-9201(18)30026-8 https://doi.org/10.1016/j.pepi.2018.07.008 PEPI 6174
To appear in:
Physics of the Earth and Planetary Interiors
Received Date: Accepted Date:
17 January 2018 13 July 2018
Please cite this article as: Foley, B.J., On the dynamics of coupled grain size evolution and shear heating in lithospheric shear zones, Physics of the Earth and Planetary Interiors (2018), doi: https://doi.org/10.1016/j.pepi. 2018.07.008
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
On the dynamics of coupled grain size evolution and shear heating in lithospheric shear zones Bradford J. Foley∗ Department of Geosciences, Pennsylvania State University, University Park PA
Abstract The formation of weak, localized shear zones in the lithosphere via tectonic forces is a requirement for the operation of plate tectonics on Earth. However, how shear zones form is not well understood. Shear heating and grain size reduction are two commonly cited mechanisms, and both processes are likely to be active in lithospheric shear zones. However, shear heating and grain size reduction can potentially act against each other, as elevated temperatures increase grain growth rates, and fast grain growth can eliminate weak shear zones formed via grain size reduction. Here, an effectively one-dimensional simple shear model that considers shear heating and grain size evolution is used to determine how grain size reduction and shear heating interact. The results indicate that for constant velocity boundary conditions at typical tectonic strain rates and mid- to lower-lithosphere temperatures, grain size reduction suppresses shear heating by weakening the shear zone. As a result shear heating at these conditions is small (∼ 10 K or less), unless grain size reduction via damage is significantly less effective than the models assume. ∗
Corresponding author Email address: [email protected] (Bradford J. Foley)
Preprint submitted to Physics of the Earth and Planetary Interiors
July 17, 2018
For constant stress boundary conditions, grain size reduction enhances shear heating and lowers the stress threshold required to trigger thermal runaway instabilities, as previous studies have found. Finally, combining shear heating and grain size reduction allows weak shear zones to form over a wider range of conditions than would result from either mechanism acting in isolation; coupling between shear heating and grain size reduction does not impede weak shear zone formation. Keywords:
1
1. Introduction
2
Why Earth is alone among the terrestrial planets of our solar system in
3
operating under a plate tectonic, rather than a stagnant lid, mode of man-
4
tle convection has long been an open question. In order for plate tectonics
5
to operate, mantle convective forces must be capable of forming weak, lo-
6
calized shear zones in the mid- to lower-lithosphere that effectively act as
7
plate boundaries (e.g. Tackley, 2000a; Bercovici et al., 2015). Without such
8
shear zones, the lithosphere is too strong to sink under its own weight and
9
form subduction zones. In particular, lithospheric strength is highest in the
10
mid- to lower-lithosphere, where deformation is predominantly ductile (e.g.
11
Kohlstedt et al., 1995); here the high viscosity of the lithosphere, due to the
12
temperature-dependence of rock rheology, produces a strong viscous resis-
13
tance to plates sliding past each other, or subducting beneath each other.
14
Forming low viscosity, ductile shear zones in this region is thus essential for
15
plate tectonics to operate. Moreover, localization of deformation into narrow
16
shear zones, that leave the interiors of plates stable and non-deforming, is
2
17
a defining characteristic of plate tectonics. Indeed, localized, ductile shear
18
zones are commonly seen in exhumed sections of the lower crust or litho-
19
sphere on Earth (e.g. Poirier, 1980; White et al., 1980; Hanmer, 1988; Jin
20
et al., 1998; Warren & Hirth, 2006; Skemer et al., 2010; Linckens et al., 2011).
21
However, how such shear zones form is not well understood.
22
Weakness and localization are linked, and both are necessary for the op-
23
eration of plate tectonics, but there are different requirements for shear zone
24
weakening and shear localization. Weak ductile shear zones can form via any
25
mechanism that allows deformation to induce a reduction in rock viscosity,
26
while localization additionally requires that stress decreases with increas-
27
ing strain rate (Mont´esi & Zuber, 2002; Bercovici & Karato, 2003). This
28
study focuses primarily on the processes that control shear zone strength.
29
Many mechanisms have been proposed for forming weak shear zones, in-
30
cluding shear heating (e.g. Yuen et al., 1978; Fleitout & Froidevaux, 1980;
31
Kaus & Podladchikov, 2006; Thielmann & Kaus, 2012), the development
32
of aligned mineral fabrics (e.g. Tommasi et al., 2009; Mont´esi, 2013), and
33
grain size reduction (e.g. Braun et al., 1999; Bercovici & Karato, 2003; War-
34
ren & Hirth, 2006; Landuyt & Bercovici, 2009; Bercovici & Ricard, 2012;
35
Platt, 2015; Bercovici & Ricard, 2016; Mulyukova & Bercovici, 2017). In
36
the above mechanisms deformation either causes dissipative heating, which
37
lowers viscosity through the temperature-dependence of mantle and crustal
38
rheology (e.g. Karato & Wu, 1993; Hirth & Kohlstedt, 2003), grain size re-
39
duction, which leads to weakening when deformation occurs in a grain size
40
sensitive creep regime such as diffusion creep or grain boundary sliding (e.g.
41
Hirth & Kohlstedt, 2003; Mont´esi & Hirth, 2003), or formation of mineral
3
42
alignment fabrics or foliations that cause viscous anisotropy, where viscos-
43
ity is lower when shearing in the direction minerals are aligned along (e.g.
44
Mont´esi, 2007; Hansen et al., 2012b).
45
Ultimately, in order to understand the origin and evolution of plate tec-
46
tonics on Earth, these shear zone weakening mechanisms must be incorpo-
47
rated into mantle convection models. While previous models have typically
48
invoked failure of the lithosphere as a result of pseudoplastic yielding for
49
weak shear zone formation (e.g. Moresi & Solomatov, 1998; Tackley, 2000b;
50
van Heck & Tackley, 2008; Foley & Becker, 2009), more sophisticated mecha-
51
nisms such as shear heating (Thielmann & Kaus, 2012) and grain size reduc-
52
tion (e.g. Landuyt et al., 2008; Foley et al., 2012; Bercovici & Ricard, 2014)
53
have been recently considered. However, most previous work has focused
54
on individual weakening mechanisms operating in isolation, when in reality
55
multiple mechanisms are likely to be active simultaneously, and thus may
56
interact in potentially unexpected ways.
57
In particular, this study focuses on coupling between grain size reduction
58
and shear heating, because these processes can potentially offset each other,
59
and thus impede the formation of weak shear zones (e.g. Kameyama et al.,
60
1997). Grain size reduction is a well supported mechanism for plate bound-
61
ary formation, as it is commonly seen in exhumed shear zones in the form
62
of mylonites and ultra-mylonites (e.g. White et al., 1980; Warren & Hirth,
63
2006). However, grains grow faster at warmer temperatures than cooler tem-
64
peratures (e.g. Evans et al., 2001), meaning that any shear heating induced
65
by flow in the shear zone may lead to faster grain growth, larger grains, and
66
viscous strengthening of the shear zone, if it deforms via grain size sensitive
4
67
creep. Thus, determining how temperature and grain size co-evolve in a shear
68
zone is essential for determining the conditions under which mylonites form,
69
and more generally what processes control the strength of plate boundaries
70
on Earth.
71
The combined effects of shear heating and grain size evolution have been
72
studied previously. Kameyama et al. (1997) shows that shear heating largely
73
eliminates grain size reduction during simple shear, unless grain growth in
74
the lithosphere is slower than experiments on pure olivine samples indicate.
75
However, grain growth in the lithosphere likely is slower than what is observed
76
in monominerallic olivine experiments, due to grain boundary pinning effects
77
from secondary phases (e.g. Evans et al., 2001; Hiraga et al., 2010). Moreover,
78
secondary phases not only impede grain growth, but also enhance grain size
79
reduction (Bercovici & Ricard, 2012; Farla et al., 2013; Linckens et al., 2015;
80
Cross & Skemer, 2017). Thielmann et al. (2015) & Thielmann (2017) also
81
studied coupled shear heating and grain size reduction in the context of deep
82
earthquake formation, and included some of the effects of secondary phases
83
on grain size evolution. However, as Thielmann et al. (2015) & Thielmann
84
(2017) focus on deep earthquakes induced by runaway shear heating instabil-
85
ities, they do not study how shear heating and grain size reduction interact
86
in stably flowing, ductile shear zones, a major goal of this paper. Moreover,
87
this paper more explicitly includes the effects of secondary phases, using a
88
grain size evolution theory based on Bercovici & Ricard (2012), and uses dy-
89
namic recrystallization experiments to constrain the key model parameters
90
in the theory, as in Mulyukova & Bercovici (2017).
91
In particular, I use an effectively one-dimensional model of a shear zone
5
92
undergoing simple shear to determine the conditions under which shear heat-
93
ing is negligible, and those where it is significant, in terms of the applied stress
94
or velocity driving deformation in the shear zone, the ambient lithospheric
95
temperature, shear zone width, and the parameters governing grain size evo-
96
lution. When shear heating is negligible, grain size reduction is unaffected
97
by heating and controls the shear zone strength, while when shear heating is
98
significant it impacts the grain size in the shear zone and the overall strength.
99
Scaling laws for the steady-state grain size and temperature are developed
100
for both constant stress and constant velocity boundary conditions, and used
101
to determine the conditions where shear heating becomes significant. The
102
models are also used to determine how coupled shear heating and grain size
103
evolution influences shear zone strength, in comparison to the cases where
104
grain size reduction or shear heating are considered in isolation.
105
The paper is organized as follows: §2 gives the theoretical formulation of
106
the shear zone model, §3 gives steady-state model results and scaling anal-
107
yses for both constant velocity and constant stress boundary conditions, §4
108
shows model results and scaling analysis for the time evolution of shear zone
109
development, §5 discusses the implications of the results for plate boundary
110
formation and the possible effects of uncertainties in the model formulation,
111
and conclusions are summarized in §6.
112
2. Theory
113
2.1. Governing equations
114
A two-dimensional shear zone is modeled (Figure 1). The center lies at
115
x = 0 and boundary at x = L/2, where L is the total width of the shear 6
116
zone, L/2 is the half width, and x measures distance from the shear zone
117
center (or core). Only half of the shear zone is modeled due to symmetry.
118
The shear zone boundary has a fixed temperature, equal to the background
119
temperature in the lithosphere of Tl , and is held at either a fixed velocity, U ,
120
or stress, τl . In reality, different thermal boundary conditions could apply, as
121
the point where temperature converges to the background temperature could
122
lie beyond the edge of the deforming shear zone. As discussed in §5.3, the
123
model setup may thus underestimate the degree of shear heating. The shear
124
zone is assumed to have no pressure gradients and no flow in the x-direction,
125
and to be infinitely long in the y-direction. As a result, the velocity profile is
126
also assumed to be uniform in y, and thus the shear zone model is effectively
127
one-dimensional (i.e. there is no variation in temperature, grain size, or
128
velocity the y-direction). With these assumptions, as well as assuming the
129
viscosity is isotropic (see §2.2), conservation of momentum reduces to ∂τxy = 0, ∂x
(1)
130
where τxy is the deviatoric shear stress in the shear zone. Equation (1)
131
implies that shear stress is constant in x in the shear zone. Using the viscous
132
constitutive law τxy = 2µε˙xy = µ
∂w , ∂x
(2)
133
where µ is viscosity, ε˙xy is the strain rate, and w is the velocity in the shear
134
zone in the y-direction. Combining (1) and (2) gives ∂ ∂x
∂w µ = 0. ∂x 7
(3)
y"
At"x"="0:" w"="0" dT/dx"="0"
0"
At"x"="L/2:"" T"="Tl" w"="U"or" τ"="τl"
L/2"
x"
Figure 1: Sketch of the shear zone model setup, including a schematic velocity profile.
135
Both temperature and grain size vary within the shear zone. Evolution of
136
temperature, T , is governed by energy conservation, where T is assumed to
137
vary only in the x-direction, resulting in no advective heat flow in the shear
138
zone as a result of the assumptions for velocity outlined above, and with
139
no heat sources aside from shear heating. In reality grain growth converts
140
surface energy into heating, and acts as a heat source, but this effect is found
141
to be small for grains larger than 10−8 m (Thielmann, 2017), and is neglected
142
here. With these assumptions temperature evolves as (1 − f )Ψ ∂T ∂ 2T =κ 2 + ∂t ∂x ρCp
(4)
143
where t is time, κ is thermal diffusivity, Ψ is the rate of deformational work,
144
f is the fraction of deformational work that goes into grain size reduction
145
(and thus 1 − f is the fraction available to drive shear heating), ρ is density,
146
and Cp is the specific heat (see Table 1 for a list of all variables and Table 2 8
147
for a list of parameters and their assumed values). The deformational work
148
rate is defined as Ψ = τij ε˙ij . For a one-dimensional shear zone undergoing
149
simple shear with isotropic viscosity, as modeled here, the only non-zero
150
components of τij are τxy = τyx , and for ε˙ij are ε˙xy = ε˙yx . As a result
151
Ψ = 2τxy ε˙xy , or Ψ = 2τ ε, ˙ where τ and ε˙ are the second invariants of the
153
deviatoric stress and strain-rate tensors, respectively; these quantities are p p defined as τ = (1/2)τij τij and ε˙ = (1/2)ε˙ij ε˙ij .
154
2.2. Rheology
152
155
A purely viscous, composite rheology that includes both dislocation and
156
diffusion creep is assumed. The purely viscous rheology is chosen as this study
157
focuses on ductile shear zone formation in the lower lithosphere. The model
158
therefore ignores elasticity, which can lead to thermal runaway instabilities
159
during shear zone formation, even with constant velocity boundary condi-
160
tions, as built up elastic stress is released (Regenauer-Lieb & Yuen, 1998;
161
Kelemen & Hirth, 2007; Thielmann et al., 2015). However, at temperatures
162
greater than ∼ 1000 K, as are expected in the mid- to lower-lithosphere, and
163
at typical geologic strain rates, elasticity is negligible (Thielmann, 2017), and
164
thus the purely viscous rheology used in this study is justifiable. Note also
165
that by including models with constant stress boundaries, thermal runaway
166
instabilities are still captured in the purely viscous models presented here,
167
and thus some first order comparisons to viscoelastic models can be made
168
(see §3.1.2 & 5.2).
169
The composite rheology used in this study is given by: ε˙ = Bdis τ n + Bdif Am τ 9
(5)
170
where Bdis is the compliance for dislocation creep, Bdif the compliance for
171
diffusion creep, A is the fineness, or inverse grain size, and n and m are
172
constants. I use n = 3, rather than n = 3.5 as given in Hirth & Kohlstedt
173
(2003), as an analytical expression for stress as a function of strain rate can
174
not be obtained with n = 3.5, and such an expression is necessary for the
175
scaling analysis performed in §3.2 and Appendix B. Using n = 3 instead
176
of n = 3.5 does not significantly change the strain rate obtained from the
177
dislocation creep flow law at a stress of 1 MPa, which is in the middle of the
178
range of stresses covered in this study. I also use m = 3 as appropriate for
179
Coble creep (e.g. Hirth & Kohlstedt, 2003). The dislocation and diffusion
180
creep compliances are temperature dependent, following Bdis
181
Edis = Bdis0 exp − RT
(6)
and
Bdif
Edif = Bdif 0 exp − RT
(7)
182
where Edis and Edif are the activation energies for dislocation and diffusion
183
creep, respectively, and R is the universal gas constant. I use Edis = 530 kJ
184
mol−1 and Edif = 375 kJ mol−1 (Hirth & Kohlstedt, 2003). The composite
185
rheology will be controlled by whichever mechanism, diffusion or dislocation
186
creep, leads to the larger strain rate. The field boundary between deformation
187
controlled by diffusion creep versus that controlled by dislocation creep can
188
be found by equating the strain rates from the two mechanisms. Writing this
10
189
boundary in terms of fineness gives the field boundary fineness, Af ield , as Af ield =
Bdis0 τ n−1 Bdif 0
m1
exp
Edif − Edis . mRT
(8)
190
When fineness is larger than Af ield deformation is in the diffusion creep
191
regime, and in the dislocation creep regime when fineness is lower than Af ield .
192
2.3. Damage theory
193
The evolution of fineness, A, is modeled using a simplified version of
194
the two-phase damage theory developed in Bercovici & Ricard (2012) and
195
extended in Bercovici & Ricard (2016). Bercovici & Ricard (2012) shows
196
that deformation of a mixture of phases, predominately olivine and pyroxene
197
in mantle peridotite, quickly enters a “pinned” state, where the grain sizes
198
of both phases are controlled by the curvature of the interface between the
199
phases, or the interface roughness. In the pinned state the grain size is
200
equal to the interface roughness multiplied by the factor π/2; the factor
201
(π/2)m is thus incorporated into Bdif 0 in the diffusion creep compliance in
202
this study. With the above assumptions, the fineness, A, evolves as (e.g.
203
Foley & Bercovici, 2014; Foley et al., 2014) dA fΨ = − hAp dt γ
(9)
204
where γ is the surface tension, h is the grain growth, or healing, rate, and p
205
is a constant (p = 5 is used, as estimated for the pinned state by Bercovici
206
& Ricard (2012)). Note that p in the above formulation is analogous to q in
207
the two-phase theory of Bercovici & Ricard (2012); p and q are related as 11
208
p = q +1. The healing rate is temperature-dependent, following (e.g. Karato,
209
1989; Evans et al., 2001) Eh h = hn exp − , RT
(10)
210
where hn is a pre-exponential constant and Eh is the activation energy for
211
grain growth. The first term on the right hand side of (9) represents grain
212
size reduction (fineness growth) as a result of partitioning a fraction, f , of the
213
deformational work into surface energy (f is referred to here as the damage
214
partitioning fraction). The second term represents grain growth (fineness
215
reduction) which acts to reduce the surface energy in the rock. A steady-state
216
fineness exists where grain size reduction and grain growth are in balance.
217
Note that in this formulation, where grain size is assumed to be pinned by
218
the interface between the two mineral phases that make up the bulk rock,
219
both dislocation and diffusion creep contribute to the deformational work
220
that drives grain size reduction. While dynamic recrystallization can only
221
occur via dislocation creep (e.g. Karato, 2008), damage to, or warping of,
222
the interface between phases can occur via any solid-state creep mechanism
223
(Bercovici & Ricard, 2012).
224
The bulk grain growth, or healing, rate in a rock in the pinned state is
225
calibrated to grain growth experiments. The pre-exponential constant, hn ,
226
is calculated in terms of a standard grain growth pre-exponential constant,
227
h∗n , valid for grain growth without pinning, from the following relationship
228
(Bercovici & Ricard, 2013): h∗n (p − 1) ∗ p−3 hn = r . 500 12
(11)
Table 1: Model Variables
Variable τxy x µ ε˙xy w T t Ψ τ ε˙ Bdis Bdif A h
Meaning Shear stress Distance from shear zone center Viscosity Strain rate Velocity Temperature Time Deformational work Second invariant of stress tensor Second invariant of strain rate tensor Dislocation creep compliance Diffusion creep compliance Fineness (inverse grain size) Healing rate
Units Pa m Pa s s−1 m s−1 K s W m−3 Pa s−1 −n Pa s−1 mm Pa−1 s−1 m−1 p−1 m s−1
Equation (1) (1) (2) (2) (2) (4) (4) (4) below (4) below (4) (6) (7) (9) (10)
229
Here r∗ is the interface roughness where the pinned state is reached; for a
230
mixture of olivine and pyroxene typical of peridotite r∗ = 10−6 m (Bercovici
231
& Ricard, 2013). As in Foley & Bercovici (2014) and Foley et al. (2014), h∗n
232
is calculated for a given grain growth activation energy from the standard
233
grain growth rate of 3 × 10−15 m2 s−1 at T = 1630 K given in Bercovici &
234
Ricard (2012), yielding h∗n =
3 × 10−15 . −Eh exp R1630
(12)
235
In the pinned state, f represents the fraction of deformational work
236
that partitions into damage along the interface between phases, while Eh
237
is, strictly speaking, the activation energy for coarsening of the interface be-
238
tween phases; as pining then controls grain size, these quantities effectively
239
govern grain size evolution. However, values of these parameters for dam-
240
age to the grains themselves, when pinning is absent, could be different. I
241
thus use f ∗ , h∗ , and Eh∗ to represent damage partitioning fraction, healing 13
242
rate, and grain growth activation energy for grains in the absence of pinning.
243
Values for f and Eh are not well known, but dynamic recrystallization ex-
244
periments can be used to constrain f ∗ or Eh∗ (Rozel et al., 2011; Mulyukova
245
& Bercovici, 2017). In particular, laboratory experiments show a trend in
246
recrystallized grain size with temperature, which constrains the temperature
247
dependence of the ratio f ∗ /h∗ , but does not constrain f ∗ and Eh∗ uniquely
248
(see Appendix A). With a chosen Eh∗ , f ∗ can be inferred from experimental
249
results and is generally a function of temperature, such that the temperature
250
dependence of f ∗ /h∗ matches the experiments (Rozel et al., 2011; Mulyukova
251
& Bercovici, 2017). In this study I chose a value of Eh∗ that results in f ∗ in-
252
dependent of temperature, in order to simplify the model. This is found at
253
Eh∗ ≈ 430 kJ mol−1 , consistent with the value used in Mulyukova & Bercovici
254
(2017). With this Eh∗ and h∗n as given in (12), f ∗ ≈ 10−4 (see Appendix A).
255
The simplest assumption is then that f = f ∗ and Eh = Eh∗ , and I make this
256
assumption for Eh in this study. Mulyukova & Bercovici (2017) suggests f
257
could be lower than f ∗ , finding that f ∼ 10−5 f ∗ can fit field observations of
258
exhumed mylonites. I thus test a range of f = 10−4 − 10−9 in this study.
259
The chosen value of Eh gives hn ≈ 1.4 × 10−15 m4 s−1 .
260
2.4. Numerical solution scheme
261
Equations (4) and (9) are solved to obtain the one-dimensional profiles of
262
temperature and fineness through the shear zone as functions of time. The
263
shear zone is discretized using 401 grid points in the x−direction. Solution
264
of (9) is straightforward because fineness grid points are independent of each
265
other in the grain size evolution formulation employed in this study, where
266
fineness evolution is assumed to be only a function of the local strain rate, 14
Parameter L Tl U τl κ f ρ Cp n m Bdis0 Edis Bdif 0 Edis R γ p hn ε˙0
Table 2: Model Parameters and Assumed Values Meaning Assumed value Total shear zone width 2-200 km, reference value: 20 km Background lithosphere temperature 900-1300 K, reference value: 1000 K Imposed shear zone boundary velocity Calculated from ε˙0 (m s−1 ) Imposed shear zone boundary stress 104 − 109 Pa −6 Thermal diffusivity 10 m s−1 (ref. 4) −9 −4 Damage partitioning fraction 10 − 10 , reference value: 10−4 (ref. 2) Lithosphere density 2800 kg m−3 (ref. 4) Heat capacity 1000 J kg−1 K−1 (ref. 4) Stress exponent 3 (ref. 3) grain size sensitivity exponent 3 (ref. 1) Dislocation creep compliance constant 1.1 × 10−13 Pa−n s−1 (ref. 1) Dislocation creep activation energy 530 kJ mol−1 (ref. 1) Diffusion creep compliance constant 1.5 × 10−15 mm Pa−1 s−1 (ref. 1) Diffusion creep activation energy 375 kJ mol−1 (ref. 1) Universal gas constant 8.314 J mol−1 (ref. 4) Surface tension 1 J m−2 (ref. 5) grain growth exponent 5 (ref. 2) grain growth pre-exponential constant 1.4 × 10−15 m4 s−1 (ref. 2) Initial strain rate 10−18 − 10−9 s−1
References for chosen parameter values: ref. 1) Hirth & Kohlstedt (2003); ref. 2) derived in this paper; ref. 3) n = 3.5 in Hirth & Kohlstedt (2003), but n = 3 greatly simplifies model analysis and is chosen here. Using n = 3 instead of n = 3.5 does not significantly change the values of Bdis0 and Edis (Mulyukova & Bercovici, 2017); ref. 4) Turcotte & Schubert (2002); ref. 5) Mulyukova & Bercovici (2017)
15
Equation above (1) above (1) above (1) above (1) (4) (4),(9) (4) (4) (5) (5) (6) (6) (7) (7) (7) (9) (9) (10) below (13)
267
temperature, and stress. Thus (9) is solved at each grid point using MAT-
268
LAB’s ordinary differential equation (ODE) solver “ode15s.” The equation
269
for temperature is more difficult, as it contains a diffusion term. To solve
270
(4) the method of lines is used (Schiesser, 2012), and the resulting system
271
of ODEs is also solved with ode15s in MATLAB. As introduced in §2.1, the
272
boundary at x = 0 corresponds to the shear zone center, and at x = L/2 to
273
the shear zone edge. Due to symmetry in temperature at the shear zone cen-
274
ter, ∂T /∂x = 0 at x = 0. The shear zone edge is held at a fixed temperature
275
Tl , and thus T = Tl at x = L/2.
276
277
For a constant stress boundary condition at the shear zone edge, Ψ can be written in a straightforward manner as Ψ = 2τl ε˙ = 2Bdis τln+1 + 2Bdif Am τl2 .
(13)
278
However, for a constant velocity boundary condition, the stress varies with
279
time as temperature and fineness evolve. For constant velocity boundary
280
conditions, the initial strain rate in the shear zone, ε˙0 , is specified, and the
281
constant edge velocity U = Lε˙0 /2. For a given stress, strain rate as a function
282
of x can be obtained for a known temperature and fineness profile in the shear
283
zone from (5). Integrating ε(x) ˙ from x = 0 to x = L/2 gives the velocity at
284
the edge of the shear zone for a given stress. Thus, stress can be determined
285
for the constant velocity shear zones by finding the stress that satisfies Z
x=L/2
ε(x)dx ˙ − U = 0.
(14)
x=0
286
The MATLAB root finding function “fzero” is used to find τ that satisfies 16
287
(14) at each timestep, and Ψ is then determined from (13). The shear zone
288
is assumed to initially have a uniform temperature equal to Tl and uniform
289
fineness equal to A0 , the reference fineness; A0 = 100, corresponding to a
290
reference grain size of 1 cm, is used.
291
When high temperatures are reached as a result of shear heating, melting
292
can occur. Melting is not explicitly included in the models, which assume a
293
solid state rheology. A reasonable estimate for the melting temperature in
294
the lower lithosphere is ≈ 1600 K, which is the solidus temperature of olivine
295
at ≈ 50 km depth based on the melting curve of Hirschmann (2000). In
296
the figures shown below, models that yield temperatures above this nominal
297
melting temperature are highlighted.
298
3. Results: Steady-state fineness and temperature
299
3.1. Numerical model results
300
3.1.1. Constant velocity shear zones
301
I first consider results from shear zones with constant velocity boundary
302
conditions. Figure 2 shows typical fineness and temperature profiles through
303
shear zones at steady-state. These steady-state profiles are found to be stable
304
to perturbations in fineness (perturbations up to 100 times the steady-state
305
fineness over regions extending up to 0.1L from the shear zone center were
306
tested). At low initial strain rates, fineness is constant spatially across the
307
shear zone, and the temperature rise due to shear heating is negligible (Fig-
308
ures 2A & B). As a result, strain rate and deformational work rate are also
309
constant spatially, and there is no localization in the shear zone (Figures 2C
310
& D). At steady-state and when diffusion creep is the dominant deforma17
311
tion mechanism (as Figures 2 & 3 show is the case for the majority of the
312
model results), fineness is a monotonic function of stress (e.g. Foley et al.,
313
2012; Foley & Bercovici, 2014; Foley & Driscoll, 2016). As stress is constant
314
spatially as required by (1), fineness will also be constant and localization
315
will not occur, when shear heating is negligible. Note that although local-
316
ization does not occur for the simple damage formulation used in this study,
317
more sophisticated grain size reduction theories do allow for localization (see
318
§5.3). Higher initial strain rates then increase the deformational work rate,
319
and hence lead to higher fineness in the shear zone (Figure 2B).
320
Eventually, initial strain rates are high enough for shear heating to be-
321
come significant. In this situation the peak temperature (Tp ) occurs at the
322
shear zone center, away from the cooler fixed temperature boundaries, while
323
the opposite is true for fineness (Figures 2A & B). High temperatures en-
324
hance grain growth and decrease fineness; thus the smallest fineness occurs
325
at the shear zone center, while fineness grows towards the shear zone edge as
326
temperatures cool. The deformational work rate is also larger at the shear
327
zone core than at the margins (Figure 2D), but the effect of temperature is
328
stronger and results in the documented grain coarsening at the shear zone
329
center. When shear heating is significant, localization can occur when the
330
temperature rise is on the order of the background temperature (e.g. Bercovici
331
& Karato, 2003), as seen in the case with an initial strain rate of 10−10 s−1
332
(Figures 2C & D). Furthermore, the fineness throughout the shear zone is
333
smaller at ε˙0 = 10−10 s−1 than at ε˙0 = 10−12 s−1 , as a result of shear heating
334
inducing grain growth.
18
1200 1000 800
Strain Rate [s-1]
Fineness [m-1]
1400
A
105 104 103 102 101 100
C
10-10 10-12 10-14 10-16 10-18
0 1 2 3 4 5 6 7 8 9 10
Distance [km]
Def. Work Rate [Pa s-1]
Temperature [K]
106 1600
B
D
10-4 10-6 10-8 10-10 10-12 10-14
0 1 2 3 4 5 6 7 8 9 10
Distance [km] Strain Rate = 10-18 s-1 Strain Rate = 10-15 s-1 Strain Rate = 10-12 s-1 Strain Rate = 10-10 s-1
Figure 2: Temperature (A), fineness (B), strain rate (C), and deformational work rate (D) as a function of distance from the shear zone center (x). L/2 = 10 km, Tl = 1000 K, f = 10−4 , and other parameters as given in Table 2. Four different initial strain rates, ε˙0 , are used: ε˙0 = 10−18 s−1 (solid line), ε˙0 = 10−15 s−1 (dashed line), ε˙0 = 10−12 s−1 (dot-dashed line), ε˙0 = 10−10 s−1 (dotted line). Blue lines give the field boundary fineness for each initial strain rate, where fineness values larger than the field boundary result in deformation in the diffusion creep regime.
19
Strain Rate [s-1]
Strain Rate [s-1]
10-16 10-14 10-12 10-10
10-16 10-14 10-12 10-10
107
10-6 10-8
100 10-2 10-4 L/2 = 1 km L/2 = 10 km L/2 = 100 km
10-6 10-8 104
102
107
100 10-2 10-4 -4
10-8
105 103
f = 10 f = 10-6 f = 10-9
10-16 10-14 10-12 10-10 -1
Strain Rate [s ]
106
L/2 = 1 km L/2 = 10 km L/2 = 100 km
105 104 C
f = 10-4 f = 10-6 f = 10-9
106 105 104
103
G
108
J
104
L/2 = 1 km L/2 = 10 km L/2 = 100 km
107
103
105
102 101
106
108
F
106 104
Tl = 1100 K Tl = 1000 K Tl = 900 K
107
103
101
I
102
10-6
102
107
E
102
Tl = 1100 K Tl = 1000 K Tl = 900 K
103 101
Fineness [m-1]
104
Stress [Pa]
Tl = 1100 K Tl = 1000 K Tl = 900 K
105 104
Stress [Pa]
10-4
10
108
B
6
Stress [Pa]
100 10-2
Fineness [m-1]
A
102
107 106 105 104 103
10-16 10-14 10-12 10-10
-4
K
f = 10 f = 10-6 f = 10-9
10-16 10-14 10-12 10-10
Strain Rate [s-1]
-1
Strain Rate [s ]
Def. Work Rate [Pa s-1] Def. Work Rate [Pa s-1] Def. Work Rate [Pa s-1]
Strain Rate [s-1] 10-16 10-14 10-12 10-10
Fineness [m-1]
Temperature Rise [K] Temperature Rise [K] Temperature Rise [K]
104
Strain Rate [s-1] 10-16 10-14 10-12 10-10
Tl = 1100 K Tl = 1000 K Tl = 900 K
10-2 10-4 10-6 10-8 10-10 10-12
D
L/2 = 1 km L/2 = 10 km L/2 = 100 km
10-2 10-4 10-6 10-8 10-10 10-12
H
10-2 L 10-4 10-6 10-8 10-10 10-12
f = 10-4 f = 10-6 f = 10-9
10-16 10-14 10-12 10-10
Strain Rate [s-1]
Figure 3: Steady-state temperature rise (Tp − Tl ), fineness, stress, and deformational work rate at the shear zone core from the full numerical solution (circles) and scaling laws for the steady-state solutions, developed in §3.2.1, (lines). Background temperature is varied for a constant f = 10−4 and L/2 = 10 km (A, B, C, & D), shear zone width varied for constant f = 10−4 and Tl = 1000 K (E, F, G, & H), and f is varied for constant L/2 = 10 km and Tl = 1000 K (I, J, K, & L). Symbols representing numerical model results and lines representing the scaling laws are colored to represent Tl , L, or f as indicated in the legend for each figure. Open symbols and dotted lines denote model results where shear heating leads to peak temperature above the nominal melting temperature. The dashed lines in the figures for fineness as a function of strain rate mark the field boundary between dislocation and diffusion creep, where fineness larger than the field boundary means deformation takes place in the diffusion creep regime.
20
335
Thus, the temperature rise (Tp − Tl ) and fineness at the shear zone cen-
336
ter in steady-state show two clear trends with initial strain rate; one at low
337
strain rates where shear heating is negligible, and one at higher strain rates
338
where shear heating is significant (Figure 3). Fineness increases with strain
339
rate until a critical threshold is reached, where the high temperatures caused
340
by shear heating counteract damage, and steady-state fineness begins to de-
341
crease with increasing strain rate. Eventually temperatures become so high
342
that melting would occur, meaning further grain growth and temperature
343
increases beyond this point would not take place in nature. Temperature
344
rise increases with increasing strain rate for all strain rates, but the slope
345
becomes smaller at the initial strain rate threshold where fineness begins
346
to decrease. Here, significant shear heating causes the viscosity to decrease
347
sharply with increasing strain rate, such that stress also decreases and the
348
deformational work rate becomes only a weakly increasing function of strain
349
rate (Figure 3). As the temperature rise is controlled by Ψ, temperature rise
350
also becomes a weakly increasing function of strain rate when shear heating
351
is significant. That stress decreases with increasing strain rate when shear
352
heating is significant also explains why localization can happen in this case,
353
as Mont´esi & Zuber (2002) shows that stress being inversely proportional to
354
strain rate is a key requirement for shear localization.
355
Changing the shear zone background temperature has a minor effect on
356
temperature rise, but a large influence on steady-state fineness when shear
357
heating is negligible (Figure 3A & B); when shear heating is significant the
358
background temperature has no strong influence on both temperature rise
359
and shear zone core fineness. When shear heating is negligible, lower back-
21
360
ground temperatures suppress healing and therefore allow smaller steady-
361
state grain sizes (larger fineness) to be reached. However, when shear heat-
362
ing is significant the temperature in the shear zone is controlled by shear
363
heating rather than the background temperature, and thus grain sizes ap-
364
proach the same value for different background temperatures. Varying shear
365
zone width has a strong influence on the temperature rise, and, when shear
366
heating is significant, on the shear zone core fineness as well (Figure 3E &
367
F). The larger L, the larger the temperature rise because the temperature
368
gradient across the shear zone is smaller. As a result, when shear heating is
369
significant, fineness is lower for wider shear zones, because the higher tem-
370
peratures enhance grain growth. When shear heating is negligible, fineness
371
is unaffected by shear zone width because in these models the shear zone
372
boundary velocity is varied with L to produce a given initial strain rate (i.e.
373
varying L does not change the initial strain rate). Increasing the amount
374
of deformational work that goes into grain size reduction, f , leads to larger
375
fineness in the shear zone. As a result, the effective viscosity of the shear
376
zone is decreased. As Ψ = 2τ ε˙0 and τ = 2µε˙0 , where µ is the effective vis-
377
cosity, decreasing µ decreases Ψ. Increasing f therefore also decreases the
378
temperature rise (Figures 3I-L). Grain size reduction thus suppresses shear
379
heating by weakening the shear zone.
380
3.1.2. Constant stress shear zones
381
With a constant stress boundary condition, a runaway shear heating in-
382
stability is possible above a critical driving stress, where heating leads to
383
higher strain rates in the shear zone, which induce further heating and in-
384
crease strain rates even more; the instability continues until the shear zone 22
385
experiences melting (e.g. Schubert & Yuen, 1978; Yuen & Schubert, 1979;
386
Kameyama et al., 1999; Turcotte & Schubert, 2002; Kaus & Podladchikov,
387
2006). As shown by Thielmann et al. (2015), grain size reduction lowers
388
the critical stress where runaway heating occurs, because grain size reduc-
389
tion weakens the shear zone, thereby enhancing deformational work rate and
390
shear heating. The results of this study confirm that grain size reduction low-
391
ers the threshold stress where thermal runaway occurs, as shown in §3.3.2.
392
In this section model results that reach a steady-state are shown; any mod-
393
els that undergo thermal runaway have no steady-state solution and are not
394
shown. As with the constant strain rate shear zones, tests of the stability of
395
the steady-states reached were performed. Steady-states are stable as long as
396
perturbations are not so large that they induce thermal runaway by weaken-
397
ing the shear zone. Any perturbations that do not trigger thermal runaway
398
decay and the unperturbed steady-state is recovered.
399
Consistent with the constant velocity case, shear localization does not
400
occur when shear heating is negligible (Figure 4). In fact, significant local-
401
ization is not seen in any steady-state solution, because large degrees of shear
402
heating are only found in models that experience thermal runaway, and thus
403
do not reach a steady-state. When thermal runaway is ongoing very high val-
404
ues of temperature and fineness in the shear zone core do lead to significant
405
shear localization.
406
Results for steady-state temperature rise and fineness at the shear zone
407
core (Figure 5) are similar to those from the constant velocity models at low
408
initial strain rates. Increasing stress leads to higher fineness and temperature
409
at the shear zone core, while lower background temperatures and higher
23
A
1010 1000
106
Fineness [m-1]
Temperature [K]
1020
C
10-12 10-14 10-16 -18
0 1 2 3 4 5 6 7 8 9 10
Distance [km]
Def. Work Rate [Pa s-1]
Strain Rate [s-1]
10-10
10-20
104 103
B
102 101 100
990
10
105
10-4
D
10-6 10-8 10-10 10-12 10-14 10-16
0 1 2 3 4 5 6 7 8 9 10
Distance [km] Stess = 104 Pa Stress = 105 Pa Stress = 4 × 105 Pa Stress = 7 × 105 Pa
Figure 4: Profiles of temperature (A), fineness (B), strain rate (C), and deformational work rate (D) across a shear zone with L/2 = 10 km, Tl = 1000 K, f = 10−4 , and other parameters as given in Table 2. Four different stresses, τl , are used: τl = 104 Pa (solid line), τl = 105 Pa (dashed line), τl = 4 × 105 Pa (dot-dashed line), τl = 7 × 105 Pa (dotted line). Blue lines give the field boundary fineness for each stress.
24
10-8
Tl = 1100 K Tl = 1000 K Tl = 900 K
10-10
100
E
10-2 10-4 10-6 10-8
L/2 = 1 km L/2 = 10 km L/2 = 100 km
10-10 102 100 10-2 10-6 10-8
f = 10-4 f = 10-6 f = 10-9
10-10 104
105
106
107
Stress [Pa]
105 Tl = 1100 K Tl = 1000 K Tl = 900 K
103 102
106
F
10-10 G 10-12 10-14 10-16 10-18 10-20 10-22 10-24
105 104
L/2 = 1 km L/2 = 10 km L/2 = 100 km
103 102 101
I
10-4
B
104
104 10-10 C 10-12 -14 10 10-16 10-18 10-20 10-22 10-24
107
101
Fineness [m-1]
102
106
106
f = 10-4 f = 10-6 f = 10-9
J
105 104 103 102 101
104
105
106
107
Stress [Pa]
10-10 K 10-12 10-14 10-16 10-18 10-20 10-22 10-24 104
105
106
107
Tl = 1100 K Tl = 1000 K Tl = 900 K
L/2 = 1 km L/2 = 10 km L/2 = 100 km
f = 10-4 f = 10-6 f = 10-9
105
106
107
Stress [Pa]
Stress [Pa] Def. Work Rate [Pa s-1] Def. Work Rate [Pa s-1] Def. Work Rate [Pa s-1]
10-6
106
105
Strain Rate [s-1]
10-4
104
Strain Rate [s-1]
10-2
Stress [Pa]
Stress [Pa]
107
Strain Rate [s-1]
106
Fineness [m-1]
105
Fineness [m-1]
Temperature Rise [K]
Temperature Rise [K]
Temperature Rise [K]
Stress [Pa] 104 102 A 100
104 10-4 D 10-6 10-8 10-10 10-12 10-14 10-16 10-18 10-20
105
107
Tl = 1100 K Tl = 1000 K Tl = 900 K
10-4 H 10-6 10-8 10-10 10-12 10-14 10-16 10-18 10-20 10-4 L 10-6 10-8 10-10 10-12 10-14 10-16 10-18 10-20 104
106
L/2 = 1 km L/2 = 10 km L/2 = 100 km
f = 10-4 f = 10-6 f = 10-9
105
106
107
Stress [Pa]
Figure 5: Same as Figure 3, but for constant stress boundary conditions rather than constant velocity. Scaling laws for the steady-state shear zone center temperature and fineness are developed in §3.2.2, and strain rate (C,G, & K) is strain rate at the shear zone core.
25
410
damage partitioning fractions lead to larger fineness. The effect of shear zone
411
width on temperature rise and fineness is also consistent with the constant
412
velocity models. However, one major difference is that with a constant stress
413
boundary condition, a lower shear zone viscosity leads to higher rates of
414
deformational work, rather than lower. At fixed τl , lower viscosity causes a
415
higher strain rate, and thus the product of stress and strain rate is higher. As
416
a result, increasing background temperature or f causes the temperature rise
417
to increase, rather than decrease as in the constant velocity models (Figure
418
5A, D, I & L). Thus, with constant stress boundary conditions, damage
419
enhances, rather than hinders, shear heating.
420
3.2. Scaling analysis
421
3.2.1. Constant velocity shear zones
422
The model results can be understood by a simple scaling analysis. The
423
full numerical models demonstrate that shear localization does not occur,
424
and thus deformational work rate, Ψ, is constant in x, when shear heating
425
is negligible (Figure 2). I thus assume Ψ is constant in space for this scal-
426
ing analysis. Localization can occur when shear heating is significant, but
427
the scaling laws derived here are found to provide a good approximation of
428
the numerical model results even when shear heating is important. With a
429
constant Ψ, (4) can be integrated to give the peak temperature at the shear
430
zone core, Tp : Tp = Tl +
(1 − f )ΨL2 . 8k
(15)
431
Deformational work is then defined at the peak temperature and steady-state
432
fineness in the shear zone core. The model results nearly all lie within the 26
433
diffusion creep regime as a result of extensive grain size reduction. I thus
434
neglect dislocation creep in the following derivation; scaling laws derived
435
using the full composite rheology are given in Appendix B. For constant
436
velocity boundary conditions, Ψ is written in terms of strain rate as Ψ = 2τ ε˙0 =
2ε˙20 . Bdif Am
(16)
437
The steady-state fineness at the shear zone center can then be found from
438
(16), and (7) - (10) as, As =
439
2f ε˙20 γBdif 0 hn
1 p+m
exp
Eh + Edif . RTp (p + m)
(17)
The equation for peak temperature, combining (15) - (17), and (7) is (1 − f )L2 Tp − Tl = 4k
2f γhn
m − p+m
ε˙20 Bdif 0
p p+m
exp
pEdif − mEh . (18) RTp (p + m)
440
(18) is a transcendental equation for Tp that can be solved numerically (in
441
this work the fzero function in MATLAB is used). With Tp found, the shear
442
zone core steady-state fineness is directly computed from (17).
443
The scaling laws for Tp and As are shown in Figure 3, and match the
444
model results well. In order to match the full range of results obtained from
445
the numerical shear zone models, the scaling laws derived in Appendix B with
446
a full composite rheology are shown in Figure 3. However, scaling laws with
447
a composite rheology only differ from those that neglect dislocation creep,
448
(17) & (18), for the limited range of conditions where deformation is in the
449
dislocation creep regime (e.g. at ε˙0 > 10−12 s−1 for the f = 10−9 models 27
450
in Figure 3J). At low strain rates, where shear heating is negligible, the
451
scaling laws and numerical model results match exactly. At high strain rates
452
the scaling laws underestimate Tp and As , because modest shear localization
453
in the shear zone core causes higher rates of deformational work in this
454
region than the scaling analysis predicts. Nevertheless, (17)-(18) provide a
455
good first order estimate of the shear zone core temperature and fineness at
456
steady-state.
457
When shear heating is negligible Tp ≈ Tl , and thus Tp can be replaced
458
with Tl in the right hand sides of (17) - (18), and Tp − Tl and As can be
459
computed directly. The trends in Tp and As with imposed strain rate can
460
then be easily understood from the scaling laws. Fineness scales with strain
461
rate as As ∼ ε˙0
462
rise scales as Tp −Tl ∼ ε˙0
463
strain rate than fineness, increasing strain rate will eventually cause shear
464
heating to become significant and impede grain size reduction, as the model
465
results show.
2/(p+m)
when shear heating is negligible, while the temperature 2p/(p+m)
. As temperature increases more sharply with
466
When Tp − Tl ∼ 10, the shear zone compliance, and hence effective
467
viscosity, begins to be significantly affected by shear heating, and Tp is no
468
longer approximately equal to Tl . The weak increase in Tp − Tl seen in
469
the numerical models at these high strain rate conditions is captured by the
470
scaling law for Tp . From (18), there is a competition between increasing strain
471
rate, which acts to increase shear zone temperature as described above, and
472
the resulting increase in temperature which acts to dampen shear heating by
473
decreasing viscosity and in turn stress. For parameters used in this study,
474
the result of this competition is the weak increase in temperature rise with
28
475
strain rate seen from the numerical models. Likewise the decrease in fineness
476
with strain rate when shear heating is significant is captured by (17) through
477
a similar competition; increasing strain rate acts to increase fineness, but
478
higher temperatures act to decrease fineness as a result of both faster grain
479
growth and a lower viscosity. In the end the temperature effect wins out,
480
resulting in the decreasing fineness trend seen in the models.
481
The scaling laws also demonstrate how other model parameters and shear
482
zone characteristics control the steady-state temperature and fineness. For
483
example, Tp − Tl scales as L2 , because a wider shear zone results in slower
484
diffusive heat transport, and thus a higher temperature. Fineness is not
485
directly influenced by L, which does not enter (17), but will be indirectly
486
influenced by variations in shear zone temperature caused by variations in L.
487
Shear zone center fineness increases with the damage to healing ratio, f /hn ,
488
as As ∼ (f /hn )1/(p+m) , because a higher damage to healing ratio means
489
more deformational work goes into reducing grain size and/or grain growth
490
is slower. Damage to healing ratio also enters into the scaling law for tem-
491
perature rise, as increasing fineness acts to dampen shear heating; (18) shows
492
that Tp − Tl ∼ (f /hn )−m/(p+m) .
493
3.2.2. Constant stress shear zones
494
Scaling laws for constant stress boundary conditions can be derived in a
495
similar fashion to those in §3.2.1. Again assuming that Ψ is constant in x
496
(Figure 4) and diffusion creep dominates, Ψ = 2Bdif Am τl2 .
29
(19)
497
Defining Ψ at the shear zone core temperature and fineness, As =
498
2f Bdif 0 τl2 γhn
1 p−m
exp
Eh − Edif RTp (p − m)
(20)
and (1 − f )L2 Tp = Tl + 4k
2f γhn
m p−m
τl2 Bdif 0
p p−m
exp
mEh − pEdif . RTp (p − m)
(21)
499
For most conditions where a steady-state solution exists, Tp ≈ Tl . Only at
500
stresses just before thermal runaway takes place does significant shear heating
501
of 10-100 K occur. Thus for the discussion of steady-state shear zones at low
502
stress in this section, Tp = Tl in the right hand side of (20) - (21) is assumed.
503
The scaling laws show that As ∼ τl
504
the increasing trends in As and Tp with stress seen in the numerical models.
505
Steady-state shear zone core fineness is an increasing function of f /hn , as in
506
the constant velocity case, scaling as As ∼ (f /hn )1/(p−m) . Furthermore, peak
507
temperature increases with increasing damage to healing ratio, scaling as
508
Tp ∼ (f /hn )m/(p−m) , because lowering shear zone viscosity via more extensive
509
damage increases the rate of deformational work (see §3.1.2). The preceding
510
reasoning also explains why grain size reduction lowers the threshold stress
511
where thermal instability occurs, as shown in Thielmann et al. (2015). The
512
dependences of fineness and peak temperature on background temperature
513
are the result of two competing effects. For fineness, increasing temperature
514
acts to decrease viscosity, and thus increase Ψ, but also increases the grain
515
growth rate. Thus the difference between the grain growth and diffusion
516
creep viscosity activation energies controls which of these two competing
2/(p−m)
30
2p/(p−m)
and Tp ∼ τl
, explaining
517
affects dominates, as (20) demonstrates. In the case of peak temperature,
518
the same competing effects of varying Tl are at play, and the sign of the term
519
mEh − pEdif in (21) determines which competing effect wins out. For the
520
baseline parameter values chosen here, mEh − pEdif is negative, and thus
521
increasing Tl increases Tp .
522
3.3. Steady-state solutions and the onset of significant shear heating
523
3.3.1. Constant velocity shear zones
524
Using the scaling laws for temperature and fineness at the shear zone
525
core, (17) and (18), steady-state fineness and temperature for constant ve-
526
locity shear zones are determined for a range of strain rates, background
527
temperatures, and two different damage partitioning fractions (Figure 6).
528
The results show the same trends seen in the numerical models, and criti-
529
cally, map out the conditions under which shear heating becomes significant.
530
Shear heating significantly affects shear zone rheology when fineness begins
531
to decrease with increasing strain rate, and thus the onset of significant shear
532
heating is represented by the peak in fineness with increasing strain rate, for
533
a given background temperature, seen in Figure 6. The fineness peak oc-
534
curs at lower values of ε˙0 as Tl decreases, though the influence of Tl is weak.
535
Lowering f shifts the fineness peak to lower ε˙0 , as smaller strain rates are
536
needed to cause significant shear heating when damage is less effective. The
537
fineness peak is found to occur at a constant temperature rise of ≈ 20 − 30
538
K. The onset of significant shear heating can thus be described by defining
539
a critical temperature rise, ∆Tcrit , above which shear heating is important.
540
Setting Tp − Tl = ∆Tcrit in (18) one can solve for the critical strain rate, ε˙crit ,
31
-14 -16
-10 2
-12
1 0 -1 -2 -3 -4 -5 -6
4
-18 900 1000 1100 1200 1300 Background Temperature [K]
B
3
-14 -16
f = 10-9 -10
C
-12
2
1 0 -1 -2 -3 -4
-18 900 1000 1100 1200 1300 Background Temperature [K]
-7 -6 -5 -4 -3 -2 -1 0 1 2 3 log10(Temperature Rise [K])
-14 4
log10(Strain Rate [s-1])
-12
-16
-18 900 1000 1100 1200 1300 Background Temperature [K]
f = 10-9 -10
-14
log10(Strain Rate [s-1])
-12
A
-16
3
D
-18 900 1000 1100 1200 1300 Background Temperature [K]
2
3
4
5
log10(Strain Rate [s-1])
-10
f = 10-4
5
log10(Strain Rate [s-1])
f = 10-4
6
log10(Fineness [m-1])
Figure 6: Steady-state temperature rise (A) and fineness (B) at the shear zone core for f = 10−4 , and for f = 10−9 (C & D). Dashed line marks the onset of significant shear heating, as given by (22), and the dot-dashed line in C & D shows where temperatures reach the nominal melting temperature (i.e. melting would occur above this line). Other model parameters are given in Table 2. The perceptually-uniform color map “roma” is used in this study to prevent visual distortion of the data (Crameri, 2018b,a).
32
Half Width = 104 m
-14
4.5
-12
3.5
3.0
-11
A -9
-8
-7
-6
-5
-4
-3
B
1200 -12
-13
4.0
1300
1100
1000
900
-13
log10(Half Width [m])
-15
-11
5.0
Background Temperature [K]
Background Temp. = 1000 K
-9
-8
-7
log10(f)
-6
-5
-4
-3
log10(f)
-15
-14
-13
-12
-11
-10
log10(Strain Rate [s-1]) Figure 7: Critical strain rate marking the onset of significant shear heating, ε˙crit , as a function of f and shear zone half width, L/2 (A), and f and background temperature, Tl (B).
541
marking the fineness peak as ε˙crit =
542
2∆Tcrit k (1 − f )L2
p+m 2p
4f γhn
m 2p
1 2
(2Bdif 0 ) exp
mEh − pEdif . (22) 2pR(Tl + ∆Tcrit )
Equation (22) is plotted in Figure 6 assuming ∆Tcrit = 20 K.
543
Tradeoffs between how damage and shear zone parameters control whether
544
significant shear heating occurs can be assessed from (22). For example,
545
ε˙crit ∼ L−(p+m)/p , demonstrating that a narrow shear zone requires a higher
546
driving strain rate in order for shear heating to be important, due the influ-
547
ence of shear zone width on heat conduction (e.g. see Figure 3F). Likewise
548
ε˙crit ∼ (f /hn )m/(2p) , as the higher the damage to healing ratio, the more shear 33
549
heating is suppressed by grain size reduction, and thus the larger strain rate
550
must be for significant shear heating.
551
Figure 7A shows ε˙crit for variable f and L. For the baseline value of f =
552
10−4 , the critical strain rate for shear heating to become significant is larger
553
than typical tectonic strain rates. For a 5 cm yr−1 boundary velocity, tectonic
554
strain rates range from 1.5×10−12 −1.5×10−14 s−1 for L/2 = 1−100 km. Thus
555
tectonic strain rates are insufficient to induce significant shear heating when
556
f = 10−4 in the purely viscous models presented here. With the lower bound
557
of f = 10−9 , shear zones of 1 km half width or smaller still require strain
558
rates higher than typical tectonic strain rates for significant shear heating,
559
though wider shear zones can experience shear heating at tectonic strain
560
rates; however, even in this case, the temperature rises are typically low, on
561
the order of 10 K. Only shear zones on the order of 100 km wide, as wide as
562
an entire plate boundary (e.g. Bird, 2003), would experience shear heating
563
of 100 K or more as a result of lower lithosphere ductile flow based on these
564
models. The background temperature also influences the onset of significant
565
shear heating, though with a lesser influence than shear zone width. As
566
Tl decreases ε˙crit also decreases, and thus shear zones are more likely to
567
experience appreciable shear heating at lower background temperatures, as
568
also seen in Figures 3 & 6.
569
The constant velocity shear zone results clearly demonstrate that grain
570
size reduction dampens shear heating, by lowering the effective shear zone
571
viscosity, and hence lowering the deformational work rate. As a result, an
572
important implication is that ignoring grain size reduction in shear zone
573
models leads one to significantly over predict the magnitude of shear heating.
34
-10 -12 -14 -16
A
Damage, f = 10-9 3
Melt
2
1 0 -1 -2 -3 -4 -5 -6
No Damage 3
B
Melt
-10 -12
2
2
1
-14
1
0
0
-1
-2
-3 -4
-18 900 1000 1100 1200 1300 Background Temperature [K]
C
-5
900 1000 1100 1200 1300 Background Temperature [K]
-1 -2 -3 -4
-16
-18 900 1000 1100 1200 1300 Background Temperature [K]
log10(Strain Rate [s-1])
log10(Strain Rate [s-1])
Damage, f = 10-4
-7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 log10(Temperature Rise [K])
Figure 8: Temperature rise caused by shear heating when grain size evolution is considered with f = 10−4 (A) and with f = 10−9 (B), and for a constant grain size of 1 cm (i.e. with no damage considered) (C). Shear zone half width is L/2 = 10 km, and other parameters are given in Table 2. Above the dot-dashed line melting would occur in the shear zone.
574
Figure 8 compares the steady-state temperature rise for shear zones where
575
grain size reduction is taken into account, and those where grain size is
576
assumed to be constant. Shear heating causes a temperature rise of 100 K at
577
strain rates of ε˙ ∼ 10−12 −10−14 s−1 for background temperatures of 1300-900
578
K, respectively, when grain size is fixed, while with grain size reduction shear
579
heating at such conditions is suppressed. A temperature rise of 100 K is not
580
seen until strain rates reach ∼ 10−11 s−1 with f = 10−4 or ∼ 10−12 − 10−13
581
with f = 10−9 , due to the effects of grain size reduction.
582
3.3.2. Constant stress shear zones
583
For constant stress shear zones, steady-state fineness and temperature
584
rise are calculated from (20) and (21) and shown in Figure 9. The region
585
where thermal runaway occurs, and no steady-state solution exists, is also
586
shown. As in the constant velocity case, the trends in Tp and As with Tl , 35
f = 10-4
A
7
8
B
Thermal Runaway
Thermal Runaway
7
6
6 2 0
5
5
-2 -4 -6 -8
5
log10(Stress [Pa])
log10(Stress [Pa])
8
f = 10-4
4
4 4 900 1000 1100 1200 1300 900 1000 1100 1200 1300 Background Temperature [K] Background Temperature [K]
log10(Stress [Pa])
8
Thermal Runaway
7 6 5
f = 10-9 Thermal Runaway
C
2 0
8
D
7
4
-2 -4 -6 -8 -10 -12 -14
6
3
5
2
log10(Stress [Pa])
f = 10-9
4 4 900 1000 1100 1200 1300 900 1000 1100 1200 1300 Background Temperature [K] Background Temperature [K]
-14-12-10 -8 -6 -4 -2 0 2 log10(Temperature Rise [K])
1
2
3
4
5
6
log10(Fineness [m-1])
Figure 9: Steady-state temperature rise (A) and fineness (B) at the shear zone core for constant stress boundary conditions and f = 10−4 . Where thermal runaway occurs (dark red region) no steady-state solution exists. Steady-state temperature rise (C) and fineness (D) are also shown for f = 10−9 . Dashed line denotes where thermal runaway occurs. Shear zone half width is L/2 = 10 km, and other parameters are given in Table 2.
36
587
τl , and f are the same as previously discussed. A salient feature of Figure 9
588
is that lowering f significantly increases the stress required to drive thermal
589
runaway. The critical stress needed to induce thermal runaway can be found
590
as follows. First, (21) is rewritten as (1 − f )L2 F 2 = Tl + 4k
2f γhn
m p−m
τl2 Bdif 0
p p−m
exp
mEh − pEdif RTp (p − m)
− Tp , (23)
591
where F2 = 0 when a steady-state solution for Tp is found. When steady-
592
state solutions exist, there are two solutions, one at low Tp that matches the
593
numerical model results, and one at higher Tp . Between these two steady-
594
state solutions lies a minimum in F2 , that is less than 0. With increasing
595
τl , the minimum of F2 increases, and eventually reaches 0; at this point
596
both steady-state solutions occur at the same Tp . For τl any higher than
597
this critical point, no steady-state solutions exist and the system goes into
598
thermal runaway. Assume that Tpcrit is a solution to (21) (i.e. it satisfies
599
F2 = 0) at the critical stress for the onset of thermal runaway, τcrit . Tpcrit
600
will also be a minimum of F2 , and thus τcrit can be found from dF2 /dTp = 0,
601
evaluated at Tp = Tpcrit :
τcrit =
−1 Bdif20
2 k 4R(p − m)Tpcrit (pEdif − mEh )(1 − f )L2
p−m 2p
2f hn γ
m − 2p
exp
pEdif − mEh . 2pRTpcrit (24)
602
Tpcrit is then found by plugging (24) into (21), at Tp = Tpcrit and τl = τcrit ,
603
and solving for Tpcrit : Tpcrit =
pEdif − mEh 2R(p − m)
4RTl (p − m) 1− 1− pEdif − mEh 37
21 ! .
(25)
604
Together, (24) and (25) give τcrit as a function of Tl , L, and parameters
605
for diffusion creep rheology and grain size evolution. Equations (24) and
606
(25) demonstrate how grain size evolution alters the critical stress threshold
607
where runaway thermal instability occurs. As previously explained, decreas-
608
ing grain size softens the shear zone, thereby enhancing the rate of deforma-
609
tional work and causing thermal runaway to occur at lower stresses than it
610
would without grain size reduction (see also Thielmann et al., 2015). Equa-
611
tion (24) captures this effect, showing that τcrit ∼ (f /hn )−m/(2p) , such that
612
with increasing damage to healing ratio thermal runaway occurs at a lower
613
driving stress. Likewise increasing Bdif 0 , which would lower the bulk shear
614
zone viscosity, also lowers τcrit . The critical stress is also inversely propor-
615
tional to L, as a wider shear zone allows for more heating and thus can more
616
easily experience thermal runaway.
617
3.4. Implications for shear zone strength and the operation of plate tectonics
618
A critical motivating question of this work is whether coupling between
619
shear heating and grain size reduction impedes weak shear zone formation
620
when both effects are significant. Figure 10 shows the shear zone core viscos-
621
ity and stress (assuming constant velocity boundary conditions) as a function
622
of strain rate, for a shear zone with only grain size reduction, one with only
623
shear heating, and one where shear heating and grain size reduction are com-
624
bined. At low strain rates, grain size reduction produces lower stresses and
625
viscosities, i.e. a weaker shear zone, than shear heating. Likewise at high
626
strain rates the opposite is true, and shear heating leads to greater weakening
627
than grain size reduction. With combined shear heating and grain size reduc-
628
tion, shear zone strength is controlled by whichever mechanism provides the
38
100 10-2 10-4 10-6
Stress [Pa]
109
A C
108 107 106 105 104 10-18 10-16 10-14 10-12 10-10 Strain Rate [s-1]
Fineness [m-1]
102
Strain Rate [s-1] 10-18 10-16 10-14 10-12 10-10 107
B
106 105 104 103
Viscosity [Pa s]
Temperature Rise [K]
Strain Rate [s-1] 10-18 10-16 10-14 10-12 10-10 104
1026 D 1024 22 10 1020 1018 1016 1014 1012 10-18 10-16 10-14 10-12 10-10 Strain Rate [s-1]
Shear heating and damage scaling law Damage only scaling law Shear heating only scaling law Shear heating and damage model Damage only model Shear heating only model
Figure 10: Temperature rise caused by shear heating (A), shear zone core fineness (B), stress (C), and viscosity (D) as a function of strain rate. Numerical shear zone model results are shown as symbols, and scaling laws as lines. A case where only shear heating is considered (diamonds and dotted line), where only grain size reduction (damage) is considered (squares and dashed line), and one where damage and shear heating are combined (cirles and solid line). Boundary conditions are constant velocity, f = 10−4 , L/2 = 10 km and Tl = 1000 K.
39
629
greatest degree of weakening, and hence the least resistance to deformation.
630
At low strain rates a shear zone where shear heating and grain size reduction
631
is combined has identical strength to a shear zone with only grain size reduc-
632
tion, and at high strain rates it has identical strength to one where only shear
633
heating is considered. At the transition point where shear heating becomes
634
significant, and thus both grain size evolution and shear heating materially
635
impact the rheology, the combined model is significantly weaker than either
636
end-member case, as grain size reduction and shear heating both contribute
637
to weakening the rock. Thus, rather than impeding each other and limiting
638
weakening, combined shear heating and grain size reduction allow weak plate
639
boundaries to form over a wider range of conditions than either mechanism
640
acting in isolation. Specifically at low strain rates grain size reduction leads
641
to significant weakening, well beyond that which would be provided by shear
642
heating, while at high strain rates shear heating induces significantly more
643
weakening than grain size reduction alone can provide.
644
4. Results: Time evolution
645
4.1. Constant stress shear zones
646
In addition to the final steady-state a shear zone will reach, I also con-
647
strain the time to steady-state to assess whether the timescales for shear zone
648
temperature and fineness evolution are geologically reasonable. I first con-
649
sider constant stress boundary conditions (§3.1.2), where the timescale for
650
reaching steady-state, or thermal runaway, can be well described by a sim-
651
ple scaling law. Time evolution results for a range of conditions are shown
652
in Figure 11; models that experience thermal runaway are stopped as shear
653
zone temperature reaches the melting temperature. 40
Temperature Rise [K]
102 100 10-2 10-4 10-6 10-8 10-10 10-12
C
102 100 10-2 10-4 10-6 10-8 10-10 10-12
E
Time [Myrs] 10-5 100 105 1010 1015 106 Fineness [m-1]
Temperature Rise [K]
A
B Tl = 900 K Tl = 1000 K Tl = 1100 K Tl = 1200 K Numerical Model Scaling Law
105 104 103 102
Fineness [m-1]
106
D τl = 105 Pa τl = 106 Pa τl = 107 Pa τl = 108 Pa Numerical Model Scaling Law
105 104 103 102
102 G 100 10-2 10-4 10-6 10-8 10-10 10-12 10-5 100 105 1010 1015 Time [Myrs]
Fineness [m-1]
106
F
105
L/2 = 1 km L/2 = 10 km L/2 = 100 km Numerical Model Scaling Law
104 103 102 106
Fineness [m-1]
Temperature Rise [K]
102 100 10-2 10-4 10-6 10-8 10-10 10-12
Temperature Rise [K]
Time [Myrs] 10-5 100 105 1010 1015
H
105
f = 10-4 f = 10-6 f = 10-9 Numerical Model Scaling Law
104 103 102 10-5 100 105 1010 1015 Time [Myrs]
Figure 11: Time evolution of temperature rise at the shear zone center, Tp − Tl , and fineness at the shear zone center. In panels A & B, Tl is varied for a constant τl = 107 Pa, f = 10−4 , and L/2 = 10 km. In panels C & D, τl is varied for a constant Tl = 1000 K, f = 10−4 , and L/2 = 10 km. Panels E & F vary L/2 with a constant τl = 107 Pa, f = 10−4 , and Tl = 1000 K, and G & H vary f for constant τl = 107 Pa, L/2 = 10 km, and Tl = 1000 K. Blue lines give the time to reach steady-state from the scaling laws developed in §4.1.
41
654
With a constant stress, fineness and temperature show a rapid rise to
655
their steady-state values, or into a thermal runaway when no steady-state
656
exists, after a long period of slow change.
657
and fineness occur at the same time, as grain size reduction enhances shear
658
heating, and vice-versa, for constant stress boundary conditions. The timing
659
of this rapid rise can be understood from the equation for fineness (9) ne-
660
glecting the grain growth term, as during the rapid rise in fineness grain size
661
reduction is the dominant process:
The rapid rises in temperature
dA fΨ = . dt γ
(26)
662
As the initial fineness is small, deformation typically begins in the dislocation
663
creep regime. I first determine the time to reach the field boundary fineness
664
by integrating (26) from A0 to Af ield for a pure dislocation creep rheology: tf ield
(Af ield − A0 )γ = exp 2f τln+1 Bdis0
Edis . RTl
(27)
665
If A0 > Af ield , meaning deformation begins in the diffusion creep regime,
666
then tf ield is set to 0. In the diffusion creep regime, (26) can be integrated
667
to give A=
668
1 1−m (1 − m)2f τl2 Bdif 0 t Edif exp − + A1−m , f ield γ RTl
(28)
which has a singularity at γA1−m f ield ts = exp (m − 1)2f τl2 Bdif 0 42
Edif . RTl
(29)
669
Runaway grain size reduction occurs as this singularity is approached, and
670
continues until grain growth becomes important and steady-state is reached,
671
or until the shear zone melts during thermal runaway. The total time for
672
rapid grain size reduction to begin is then given by ttot = tf ield + ts , which
673
matches the numerical model results well (Figure 11). Specifically, the scal-
674
ing laws accurately capture the increase in time needed for the onset of rapid
675
grain size reduction (fineness growth) caused by decreasing Tl , τl , or f . Low-
676
ering Tl , τl , or f all slow grain size reduction by either lowering the rate of
677
deformational work (Tl and τl ), or the amount of deformational work that
678
drives grain size reduction (f ). Changing L has no discernible effect on the
679
time to rapid grain size reduction, because deformational work rate is not
680
dependent on L when the driving stress is fixed.
681
Using (27) and (29), the timescale for the fineness rise is calculated for
682
a range of driving stresses and background temperatures (Figure 12). Only
683
at high stress, high temperature conditions (from above 1000 K at τl ≈ 109
684
Pa to above 1300 K at τl ≈ 107 Pa) can significant grain size reduction and
685
heating occur on geologically reasonable timescales of ∼ 100 Myrs or less,
686
for the constant stress boundary condition models employed here. Lower
687
stress and temperature conditions produce improbably long timescales for
688
significant grain size reduction to occur, or for thermal runaway to develop
689
when there is no steady-state solution. Thus, although thermal runaway is
690
predicted to occur at stresses as low as 1 MPa in this study, it may take
691
too long for the instability to develop to be seen in natural settings at these
692
conditions. However, note that additional deformation mechanisms, such as
693
geometric softening due to viscous anisotropy (e.g. Tommasi et al., 2009;
43
1200 2 2
4
1100
4
6
1000
1000 6
10
12
5
8
900
1100
6 7 8 log10(Stress [Pa])
9
-2 0 2 4 6 8 10 12 14 log10(Time [Myrs])
-20 -18 -16 -14 -12 -10 log10(Strainrate [s-1])
0
2
4
6
8
900
10 12
log10(Time [Myrs])
Figure 12: Time for rapid grain size reduction to begin as a function of stress, τl , and background temperature, Tl , as calculated from the scaling laws in §4.1. Shear zone width of L/2 = 10 km and f = 10−4 are assumed.
44
Background Temperature [K]
0
1200
8
Background Temperature [K]
1300
-2
1300
694
Hansen et al., 2012b), could be important during the early stages of shear
695
zone formation, and may speed up the development of fine grained mylonite
696
zones, or the onset of thermal runaway.
697
Dynamic recrystallization experiments constrain the damage to healing
698
ratio, but not f uniquely, so there is significant uncertainty in the exact
699
value of f . As the timing of rapid grain size reduction scales with 1/f , this
700
uncertainty maps into the timescale analysis presented in this section. For
701
example, using f = 0.1 would result in significant grain size reduction devel-
702
oping 1000 times faster than shown in Figure 12, and allow weak shear zones
703
to form in 100 Myrs or less for a broader range of stresses and background
704
temperatures. The influence of f on the timing of rapid grain size reduction
705
likely plays a major role in why Mulyukova & Bercovici (2018) find shorter
706
timescales for significant grain size reduction to occur (generally on the order
707
of 10-100 Myrs) then what is presented in Figure 12; they use f = 10−2 −10−1
708
at lithospheric temperatures and stresses on the order of 108 . Note also that
709
Mulyukova & Bercovici (2018) involves a driving stress that increases with
710
time, as they model passive margin evolution, instead of the constant stress
711
used in this study.
712
4.2. Constant velocity shear zones
713
With a constant velocity boundary condition, timescales to reach steady
714
state are similar to the constant stress case. However, unlike the constant
715
stress case, there is no grain size reduction singularity, and grain size re-
716
duction takes place more gradually from the initial grain size to the final,
717
steady-state grain size (Figure 13). As a result, a simple scaling law for the
718
time when steady-state grain size is reached is more difficult to derive, and 45
100 10-2 10-4
106 Fineness [m-1]
A
Time [Myrs] 10-5 100 105 1010 1015
102
C
100 10-2 10-4 10-6 10-8
E
102 100 10-2 10-4
D ε˙ 0 = 10-17 s-1 ε˙ 0 = 10-15 s-1 ε˙ 0 = 10-13 s-1 ε˙ 0 = 10-11 s-1
105 104 103
F
105
L/2 = 1 km L/2 = 10 km L/2 = 100 km
104 103 102
G
100 10-2 10-4 10-6
103
106
10-6 102
104
102 Fineness [m-1]
104
Tl = 900 K Tl = 1000 K Tl = 1100 K Tl = 1200 K
105
106 Fineness [m-1]
104
B
102
10-6
10-5 100 105 1010 1015 Time [Myrs]
106 Fineness [m-1]
Temperature Rise [K] Temperature Rise [K] Temperature Rise [K] Temperature Rise [K]
102
Time [Myrs] 10-5 100 105 1010 1015
H
105
f = 10-4 f = 10-6 f = 10-9
104 103 102
10-5 100 105 1010 1015 Time [Myrs]
Figure 13: Time evolution of temperature rise at the shear zone center, Tp − Tl , and fineness at the shear zone center. In panels A & B, Tl is varied for a constant ε˙0 = 10−15 s−1 , f = 10−4 , and L/2 = 10 km. In panels C & D, ε˙0 is varied for a constant Tl = 1000 K, f = 10−4 , and L/2 = 10 km. Panels E & F vary L/2 with a constant ε˙0 = 10−15 s−1 , f = 10−4 , and Tl = 1000 K, and G & H vary f for constant ε˙0 = 10−15 s−1 , L/2 = 10 km, and Tl = 1000 K.
46
719
hence not attempted here. However, the numerical results give guidance on
720
the typical timescales for grain size reduction and shear heating. The models
721
show that there are two important time scales: one where significant grain
722
size reduction occurs and there is a peak in fineness, and a longer timescale
723
to reach the final steady-state. Unlike the constant stress case, background
724
temperature has only a minor influence on the time for significant grain size
725
reduction to occur (Figure 13B), because of two competing effects: Lowering
726
the background temperature increases the deformational work rate and hence
727
rate of grain size reduction, but at the same time smaller grain sizes must be
728
reached before grain growth becomes important, at least when shear heating
729
is minor. The time to reach the final steady-state does increase modestly with
730
decreasing background temperature. Increasing the imposed strain rate, ε˙0 ,
731
has a strong effect on the time for significant grain size reduction to develop,
732
and, intuitively, lowers this timescale (Figure 13D). However, the time to
733
reach steady-state is not strongly sensitive to ε˙0 . Shear zone width has very
734
little influence on the time to reach steady-state, while damage fraction, f ,
735
does; increasing f leads to faster grain size reduction and thus shortens the
736
time for both significant grain size reduction to occur and for steady-state to
737
be reached (Figures 13F & H).
738
Overall, timescales are unreasonably long, greater than 1 Gyrs, for most
739
conditions modeled. Only at high imposed strain rates, ε˙0 ∼ 10−13 − 10−11
740
s−1 , can shear zones form on 1 − 10 Myr timescales. Although such strain
741
rates are not unreasonable for 1-10 km wide shear zones, driving a shear
742
zone at ε˙0 ∼ 10−13 − 10−11 s−1 before significant grain size reduction and/or
743
shear heating has developed requires high stresses. Specifically stresses of
47
744
τ ∼ 1 − 10 GPa are required to drive a shear zone at strain rates of ε˙0 ∼
745
10−13 − 10−11 s−1 , with a background temperature of 1000 K; such stresses
746
are higher than typical tectonic stresses (e.g. Behr & Platt, 2011). Thus,
747
as discussed in §4.1, additional weakening mechanisms that are active early
748
in shear zone formation may be necessary for mylonite zones to form on
749
geologically reasonable timescales. Note that localization might also hasten
750
shear zone development, by increasing strain rates in the localized zone.
751
Therefore more sophisticated damage theories that allow for localization, as
752
discussed in §5.3, may also be important for shear zone formation.
753
5. Discussion
754
5.1. Implications for observations of natural shear zones
755
Overall, the results demonstrate that in stably flowing shear zones that
756
can be described with constant velocity boundary conditions, like plate bound-
757
aries, and where elasticity is negligible, shear heating is diminshed by grain
758
size reduction for most of the conditions explored. Grain size reduction
759
weakens the shear zone and thus lowers the deformational work that drives
760
shear heating. At the mid- to lower-lithosphere temperatures explored in this
761
study either weakly effective damage or very high strain rates are needed for
762
shear heating of ∼ 10 K, which is sufficient to have an impact on the shear
763
zone grain size and rheology. Shear heating of 100 K or more is only seen
764
for an even narrower range of conditions. Mont´esi (2013) argues that shear
765
heating of greater than 100 K is needed for significant localization, so other
766
mechanisms, such as fabric development or feedbacks between grain size re-
767
duction and phase mixing (see §5.3), are likely to be more important for 48
768
shear localization in lower lithosphere shear zones. However, with a constant
769
stress boundary condition the results indicate that grain size reduction sig-
770
nificantly enhances the likelihood of thermal runaway instabilities, consistent
771
with previous studies (see §5.2), by lowering the stress threshold necessary
772
for triggering thermal runaway. Thus transient pulses of significant shear
773
heating at shallower depths, where elasticity is important, may be common,
774
so long as instability can develop on geologically reasonable timescales.
775
Whether significant shear heating is observed in natural ductile shear
776
zones, however, has been debated. Metamorphic assemblages in exhumed
777
fault zones and shear zones have been interpreted as showing shear heating
778
in some studies (Graham & England, 1976; England & Molnar, 1993; Leloup
779
& Kienast, 1993; Camacho et al., 2001), though alternative explanations
780
that do not require shear heating have also been proposed for many of these
781
localities (Leloup et al., 2001; Gilley et al., 2003; Kidder & Ducea, 2006;
782
Chapman et al., 2011; Kohn, 2014). If shear heating in the deeper lithosphere
783
is generally muted, then suppression of shear heating via grain size reduction
784
provides a possible explanation. This study also predicts that for constant
785
strain rate shear zones when shear heating is non-negligible, the temperature
786
distribution in the shear zone leads to grain growth in the shear zone core.
787
As a result, grain sizes are largest in the shear zone center and become
788
progressively smaller towards the shear zone edge, where temperatures are
789
cooler. Such a grain size distribution within a shear zone has not been
790
observed in nature to my knowledge, but observing such a trend would be
791
an indicator that non-negligible levels of shear heating took place.
49
792
5.2. Comparison to previous studies
793
This work reaffirms the results of Thielmann et al. (2015) & Thielmann
794
(2017), that grain size reduction lowers the critical stress needed to trigger
795
thermal runaway. In this study, grain size reduction is more extensive than
796
in Thielmann et al. (2015) & Thielmann (2017), as the pinned state formu-
797
lation for grain damage both suppresses grain growth and allows grain size
798
reduction to continue in the diffusion creep regime. As a result, I find the
799
stress threshold for thermal runaway can be lowered even further than what
800
is shown in Thielmann et al. (2015) & Thielmann (2017). However, the long
801
time scales needed for instability to develop mean high stresses of 100 MPa
802
or more may still be necessary for thermal runaway to be seen in natural
803
shear zones, as explained in §4.1.
804
Hansen et al. (2012a) compared rock deformation experiments performed
805
at constant strain rate and constant stress conditions, and found that the
806
constant stress conditions enhance shear localization and lead to unstable
807
deformation, where strain accelerates during deformation. This is similar
808
to what is found here, where constant velocity boundary conditions due not
809
allow for a runaway instability, while constant stress conditions do. However,
810
there are important differences to note as well. In Hansen et al. (2012a),
811
localization occurs due to the combined effects of grain size reduction and
812
fabric development, while here fabric development is not modeled (viscosity
813
is assumed to be isotropic) and grain size reduction alone does not lead to
814
localization, as explained in §3.1.1. Localization seen in Hansen et al. (2012a)
815
could be a transient effect that would be erased at higher strain, or a result
816
of effects not included in the grain size formulation used in this study, such
50
817
as fabric development or the influence grain size itself has on the effectiveness
818
of grain size reduction (discussed further in §5.3).
819
5.3. Model limitations
820
As the goal of this study is to constrain the first order effects of coupled
821
shear heating and grain size reduction on ductile shear zone strength, a sim-
822
ple model is used that can be easily understood and analyzed. Namely an
823
effectively one-dimensional shear zone model with a constant temperature at
824
the shear zone edge, which allows steady-state solutions to exist, is employed.
825
Other possibilities are a no heat flux boundary at the shear zone edge or far
826
from the shear zone edge (i.e. as x goes to infinity), neither of which allow
827
for steady-state solutions. Ultimately vertical heat transport to the cold sur-
828
face and to the base of the lithosphere will limit the amount of shear heating
829
that can take place, and allow a steady-state thermal profile to be reached,
830
in a way that is analogous to the simple one-dimensional model used in this
831
study. Nevertheless, it is likely that the models shown here underestimate
832
the peak temperature reached in the shear zone core and skew the onset of
833
significant shear heating to higher driving stress or strain rate, as a result
834
of the constant temperature condition applied to the shear zone edge. How-
835
ever, the overall results, that shear heating is insignificant and grain size
836
reduction is extensive at low stress or strain rate, and shear heating is sig-
837
nificant while grain size reduction is hampered at higher τl or ε˙0 , is unlikely
838
to be significantly affected by the temperature boundary condition. Likewise
839
the result illustrated in §3.4, that combined grain size reduction and shear
840
heating allows weak shear zones to form over a wider range of conditions
841
than either mechanism acting in isolation, should also hold true regardless of 51
842
the temperature boundary condition, as it is caused by the tendency of rock
843
deformation to occur via whichever mechanism allows for the highest strain
844
rate.
845
The effectively one-dimensional, simple shear model setup used also means
846
any possible vertical flow in the shear zone is not captured. If vertical flow
847
is significant it would change the deformational work rate, and hence grain
848
size and temperature evolution, in the shear zone, but ultimately this issue is
849
beyond the scope of this study as it involves a different shear zone geometry.
850
The rheology assumed is also purely viscous, when in reality elasticity can
851
be important at low temperatures. In particular, elasticity can drive thermal
852
runaway even with constant velocity boundary conditions, by releasing stress
853
that builds up during shearing. The models presented here therefore cannot
854
capture these thermal runaway scenarios, which would result in lower grain
855
sizes and higher temperatures in shear zones than shown in this study. How-
856
ever, as explained in §2.2, elasticity is negligible for most of the conditions
857
explored in this paper.
858
The models also use a constant damage partitioning fraction, f . While the
859
choice of Eh in this study eliminates any possible temperature dependence
860
for f , at least as inferred from dynamic recrystallization experiments (see
861
Appendix A), f may also be a function of grain size (Bercovici & Ricard,
862
2016; Mulyukova & Bercovici, 2017). Damage is most efficient when pyroxene
863
and olivine grains are well mixed, allowing Zener pinning to both suppress
864
grain growth and enhance grain size reduction (Bercovici & Ricard, 2012).
865
As mixing is thought to be more efficient at smaller grain sizes than larger
866
ones (Bercovici & Skemer, 2017; Cross & Skemer, 2017), the pinned state is
52
867
not reached until grains are reduced to small enough sizes for the phases to
868
easily mix. Thus there are two limits to the effective value of f ; a smaller
869
value of f characteristic of rock that is not in the pinned state when grain
870
size is large, and a larger value, characteristic of the pinned state, when grain
871
sizes are small. In this study, the pinned state is assumed to prevail at all
872
conditions, meaning that the value of f used is the value for the pinned state.
873
Allowing f to vary with grain size would allow shear localization, forming
874
a fine-grained region at the shear zone core where damage is operating in
875
the pinned state, with the rest of the shear zone seeing coarser grains and
876
low strain rates (Bercovici & Ricard, 2016). How shear heating would be
877
affected with a grain size dependent f is not known without running full
878
models, but it is likely that when significant grain size reduction and shear
879
localization occurs, shear heating would be suppressed for constant velocity
880
boundary conditions, as in the results shown in this study. The conditions
881
for the onset of siginficant shear heating may be affected as well.
882
Finally, I note that additional weakening mechanisms that are not in-
883
cluded in this study could be important. In particular, weakening beyond
884
what is provided by shear heating and damage alone may be important dur-
885
ing the initial stages of shear zone formation, and thus help speed up the
886
development of significant grain size reduction or heating. As deformation
887
likely starts outside the diffusion creep regime, viscous anisotropy that re-
888
sults from LPO formation may be one such mechanism. Furthermore, creep
889
mechanisms such as low temperature plasticity and dislocation accommo-
890
dated grain boundary sliding, neglected in this study, may have an impor-
891
tant effect as well (e.g. Hansen et al., 2012c; Thielmann, 2017). Ultimately,
53
892
understanding the deformation mechanisms prevalent early in shear zone for-
893
mation may be critical to understanding how weak shear zones develop at
894
all.
895
6. Conclusions
896
Models of ductile lithospheric shear zones constrain the co-evolution of
897
grain size and temperature for both constant velocity and constant stress
898
boundary conditions. The constant velocity models demonstrate that at low
899
strain rates, shear heating is negligible and grain size shrinks with increasing
900
strain rate. At high strain rates shear heating becomes significant, and high
901
temperatures cause grain size to increase with increasing strain rate due to
902
rapid grain growth rates. For all but very high strain rates, where high
903
temperatures lead to large grain sizes and thus push deformation into the
904
dislocation creep regime, grain size reduction diminishes the amount of shear
905
heating, by lowering the shear zone viscosity and hence the deformational
906
work available to drive heating. As a result, ignoring grain size evolution
907
leads one to significantly overestimate the amount of shear heating during
908
simple shear deformation at geologically relevant conditions. For example,
909
without considering grain size evolution, hundreds of degrees of shear heating
910
can occur at strain rates of 10−14 s−1 at mid- to lower-lithosphere background
911
temperatures, while strain rates must be at least 10-100 times larger for
912
the same level of heating when grain size evolution is considered. Thus
913
long lived, stable shear zones that can be described by constant velocity
914
boundary conditions, are less likely to experience significant shear heating
915
than previously thought, at least for the range of background temperatures 54
916
considered here.
917
With a constant stress boundary condition, grain size reduction is sig-
918
nificant and shear heating minor at low stresses, where stable steady-states
919
can be reached. At higher stresses thermal runaway instability occurs and
920
the shear zone experiences very high temperatures. Consistent with previ-
921
ous studies, the results indicate that grain size reduction lowers the stress
922
threshold where thermal runaway takes place. In fact, the stress needed to
923
trigger thermal runaway is found to be even lower than previous studies that
924
also looked at the combined effects of shear heating and grain size reduc-
925
tion, because grain boundary pinning by secondary phases, which this study
926
incorporates, enhances grain size reduction. However, the timescale needed
927
for thermal runaway to occur in the mid-lithosphere at stresses lower than
928
100 MPa is found to be greater than 100 Myrs, so high stresses may still be
929
necessary for thermal runaway to take place in practice.
930
Finally, rather than acting against each other in terms of providing rhe-
931
ological weakening, combining shear heating and grain size reduction allows
932
weak shear zones to form over a wider range of conditions than would be
933
possible from either mechanism acting in isolation. At low strain rates shear
934
heating is negligible and thus an ineffective weakening mechanism, but grain
935
size reduction is extensive. At high strain rates shear heating is significant
936
while grain size reduction is muted by the high temperatures, and shear heat-
937
ing provides substantial weakening. Thus whichever mechanism allows for
938
the least resistance to flow will control the rheology. Even when both shear
939
heating and grain size reduction are non-negligible, their combined effects
940
provide more weakening than either mechanism acting alone.
55
941
7. Acknowledgements
942
I thank Elvira Mulyukova and Dave Bercovici for discussions on experi-
943
mental constraints on damage theory, and Elvira Mulyukova for sharing some
944
of her codes related to this topic. John D. Platt also helped in the initial
945
stages of this work and showed me how to implement the method of lines for
946
shear zone modeling. I also thank Lars Hansen, Laurent Montesi, and two
947
anonymous reviewers for their comments that helped to significantly improve
948
the manuscript.
949
Appendix A. Constraining f ∗ from laboratory experiments
950
The damage partitioning fraction for damage to the grains directly, f ∗ ,
951
is constrained by laboratory data following Mulyukova & Bercovici (2017).
952
Specifically, Mulyukova & Bercovici (2017) fit experimental results to an
953
equation of the form R = Clab τ −klab ,
(A.1)
954
where R is the grain size after dynamic recrystallization and Clab is a function
955
of temperature as (their equation B.6) Clab =
3λ2 γh∗ 4λ3 f ∗ Bdis
31 .
(A.2)
956
Here λ2 ≈ 3.6 is the second moment of the lognormal grain size distribution
957
and λ3 ≈ 17.8 is the third moment. Rearranging for f ∗ ∗
f =
−3 Clab
3λ2 γh∗n 4λ3 Bdis0
56
exp
Edis − Eh . RT
(A.3)
Damage Partitioning Fraction, f *
f * = exp(-9.41(T/Tref)-0.008) f * inferred from experiments
10-3
10-4
10-5 800 1000 1200 1400 1600 1800 2000 Temperature [K]
Figure A.14: Inferred values of f ∗ from laboratory experiments assuming Eh = 430 kJ mol−1 and calculating h∗n from (12) (circles), and fit to the data using (A.4) (solid line).
958
Inferred values of f ∗ can then be determined from the experimental results
959
shown in Mulyukova & Bercovici (2017) using (A.3) with a given Eh . From
960
inspecting the inferred values of f ∗ for various choices of Eh , it was found that
961
Eh ≈ 430 kJ mol−1 results in f ∗ values that are independent of temperature
962
(Figure A.14). Fitting these f ∗ values with a function of the form f ∗ = exp −Cf
T Tref
kf ! (A.4)
963
where Tref = 1000 K is a reference temperature, yields Cf ≈ 9.41 and kf ≈
964
−0.008, and a value of f ∗ ≈ 10−4 .
965
Appendix B. Composite rheology scaling laws
966
Here I show the scaling laws for steady-state temperature and fineness at
967
the shear zone core with a composite rheology. As in the main text, the peak
57
968
temperature is given by Tp = Tl +
969
(1 − f )Ψ(Tp )L2 8k
and fineness by As =
970
(B.1)
f Ψ(Tp ) γh(Tp )
p1 .
(B.2)
With a constant stress, Ψ is easily given as 2 Ψ = 2(Bdis (Tp )τln+1 + Bdif (Tp )Am s τl ),
(B.3)
971
and combining (B.1) - (B.3) with (6), (7), and (10) gives the scaling laws for
972
Tp and As . With a constant boundary velocity, the rheology law, (5), must
973
be inverted to give stress as a function of strain rate in order to write the
974
scaling laws for Tp and As . With n = 3, τ=
975
η (2/3)1/3 Bdif (Tp )Am s − 181/3 Bdis (Tp ) η
(B.4)
where 1
1
2 3 η = [9Bdis (Tp )2 ε˙0 + {3(27Bdis (Tp )4 ε˙20 + 4Bdis (Tp )3 Bdif (Tp )3 A3m s )} ] . (B.5)
976
Using Ψ = 2τ ε˙0 , with τ given by (B.4)-(B.5), (B.1)-(B.2), and (6), (7),
977
and (10), the scaling laws for Tp and As with constant velocity boundary
978
conditions and a composite rheology can be found.
58
979
Behr, W. M., & Platt, J. P. (2011). A naturally constrained stress profile
980
through the middle crust in an extensional terrane. Earth Planet. Sci.
981
Lett., 303 , 181–192.
982
Bercovici, D., & Karato, S. (2003). Theoretical analysis of shear localization
983
in the lithosphere. In S. Karato, & H. Wenk (Eds.), Reviews in Mineralogy
984
and Geochemistry: Plastic Deformation of Minerals and Rocks chapter 13.
985
(pp. 387–420). Washington, DC: Min. Soc. Am. volume 51.
986
Bercovici, D., & Ricard, Y. (2012). Mechanisms for the generation of plate
987
tectonics by two-phase grain-damage and pinning. Phys. Earth Planet.
988
Inter., 202 , 27–55.
989
Bercovici, D., & Ricard, Y. (2013). Generation of plate tectonics with two-
990
phase grain-damage and pinning: Source-sink model and toroidal flow.
991
Earth Planet. Sci. Lett., 365 , 275–288.
992
993
994
995
996
997
Bercovici, D., & Ricard, Y. (2014). Plate tectonics, damage, and inheritance. Nature, 508 , 513–516. Bercovici, D., & Ricard, Y. (2016). Grain-damage hysteresis and plate tectonic states. Phys. Earth Planet. Inter., 253 , 31–47. Bercovici, D., & Skemer, P. (2017). Grain damage, phase mixing and plateboundary formation. Journal of Geodynamics, 108 , 40–55.
998
Bercovici, D., Tackley, P., & Ricard, Y. (2015). The generation of plate
999
tectonics from mantle dynamics. In D. Bercovici, & G. Schubert (chief
1000
editor) (Eds.), Treatise on Geophysics (pp. 271–318). New York: Elsevier
1001
volume 7, Mantle Dynamics. 59
1002
1003
Bird, P. (2003). An updated digital model of plate boundaries. Geochem., Geophys., Geosyst., 4 , 1027.
1004
Braun, J., Chery, J., Poliakov, A., Mainprice, D., Vauchez, A., Tomassi, A.,
1005
& Daignieres, M. (1999). A simple parameterization of strain localization
1006
in the ductile regime due to grain size reduction: A case study for olivine.
1007
J. Geophys. Res., 104 , 25167–25181.
1008
Camacho, A., McDougall, I., Armstrong, R., & Braun, J. (2001). Evidence
1009
for shear heating, Musgrave Block, central Australia. J. Struct. Geol., 23 ,
1010
1007–1013.
1011
Chapman, A., Luffi, P., Saleeby, J., & Petersen, S. (2011). Metamorphic evo-
1012
lution, partial melting and rapid exhumation above an ancient flat slab:
1013
insights from the san emigdio schist, southern california. Journal of Meta-
1014
morphic Geology, 29 , 601–626.
1015
1016
Crameri, F. (2018a). Geodynamic diagnostics, scientific visualisation and StagLab 3.0. Geosci. Model Dev. Discuss., .
1017
Crameri, F. (2018b). Scientific colour maps (Version 3.0.0). Zenodo, .
1018
Cross, A. J., & Skemer, P. (2017). Ultramylonite generation via phase mixing
1019
in high-strain experiments. J. Geophys. Res. Solid Earth, 122 , 1744–1759.
1020
England, P., & Molnar, P. (1993). The interpretation of inverted metamor-
1021
1022
phic isograds using simple physical calculations. Tectonics, 12 , 145–157. Evans, B., Renner, J., & Hirth, G. (2001). A few remarks on the kinetics of
60
1023
static grain growth in rocks. Int. J. Earth Sciences (Geol. Rundsch.), 90 ,
1024
88–103.
1025
Farla, R. J. M., Karato, S.-i., & Cai, Z. (2013). Role of orthopyroxene
1026
in rheological weakening of the lithosphere via dynamic recrystallization.
1027
Proceedings of the National Academy of Science, 110 , 16355–16360.
1028
1029
Fleitout, L., & Froidevaux, C. (1980). Thermal and mechanical evolution of shear zones. J. Struct. Geol., 2 , 159–164.
1030
Foley, B. J., & Becker, T. W. (2009). Generation of plate-like behavior
1031
and mantle heterogeneity from a spherical, visco-plastic convection model.
1032
Geochem., Geophys., Geosyst., 10 . Q08001.
1033
Foley, B. J., & Bercovici, D. (2014).
Scaling laws for convection with
1034
temperature-dependent viscosity and grain-damage. Geophys. J. Int., 199 ,
1035
580–603.
1036
Foley, B. J., Bercovici, D., & Elkins-Tanton, L. T. (2014). Initiation of plate
1037
tectonics from post-magma ocean thermochemical convection. J. Geophys.
1038
Res. Solid Earth, 119 , 8538–8561.
1039
Foley, B. J., Bercovici, D., & Landuyt, W. (2012). The conditions for plate
1040
tectonics on super-earths: Inferences from convection models with damage.
1041
Earth Planet. Sci. Lett., 331332 , 281 – 290.
1042
Foley, B. J., & Driscoll, P. E. (2016). Whole planet coupling between cli-
1043
mate, mantle, and core: Implications for rocky planet evolution. Geochem.,
1044
Geophys., Geosyst., 17 . 61
1045
Gilley, L. D., Harrison, T. M., Leloup, P. H., Ryerson, F. J., Lovera, O. M.,
1046
& Wang, J.-H. (2003). Direct dating of left-lateral deformation along the
1047
Red River shear zone, China and Vietnam. J. Geophys. Res. Solid Earth,
1048
108 , 2127.
1049
Graham, C. M., & England, P. C. (1976). Thermal regimes and regional
1050
metamorphism in the vicinity of overthrust faults: an example of shear
1051
heating and inverted metamorphic zonation from southern California.
1052
Earth Planet. Sci. Lett., 31 , 142–152.
1053
Hanmer, S. (1988). Great Slave Lake Shear Zone, Canadian Shield: recon-
1054
structed vertical profile of a crustal-scale fault zone. Tectonophysics, 149 ,
1055
245–264.
1056
Hansen, L. N., Zimmerman, M. E., Dillman, A. M., & Kohlstedt, D. L.
1057
(2012a). Strain localization in olivine aggregates at high temperature: A
1058
laboratory comparison of constant-strain-rate and constant-stress bound-
1059
ary conditions. Earth Planet. Sci. Lett., 333 , 134–145.
1060
Hansen, L. N., Zimmerman, M. E., & Kohlstedt, D. L. (2012b). Laboratory
1061
measurements of the viscous anisotropy of olivine aggregates. Nature, 492 ,
1062
415–418.
1063
Hansen, L. N., Zimmerman, M. E., & Kohlstedt, D. L. (2012c). The influence
1064
of microstructure on deformation of olivine in the grain-boundary sliding
1065
regime. J. Geophys. Res. Solid Earth, 117 , B09201.
1066
Hiraga, T., Tachibana, C., Ohashi, N., & Sano, S. (2010). Grain growth
62
1067
systematics for forsterite enstatite aggregates: Effect of lithology on grain
1068
size in the upper mantle. Earth Planet. Sci. Lett., 291 , 10 – 20.
1069
Hirschmann, M. M. (2000). Mantle solidus: Experimental constraints and
1070
the effects of peridotite composition. Geochem., Geophys., Geosyst., 1 .
1071
Hirth, G., & Kohlstedt, D. (2003). Rheology of the upper mantle and the
1072
mantle wedge: a view from the experimentalists. In J. Eiler (Ed.), Sub-
1073
duction Factory Mongraph (pp. 83–105). Washington, DC: Am. Geophys.
1074
Union volume 138.
1075
Jin, D., Karato, S., & Obata, M. (1998). Mechanisms of shear localization
1076
in the continental lithosphere: Inference from the deformation microstruc-
1077
tures of peridotites from the Ivrea zone, northwestern Italy. J. Struct.
1078
Geol., 20 , 195–209.
1079
Kameyama, M., Yuen, D., & Fujimoto, H. (1997). The interaction of viscous
1080
heating with grain-size dependent rheology in the formation of localized
1081
slip zones. Geophys. Res. Lett., 24 , 2523–2526.
1082
Kameyama, M., Yuen, D. A., & Karato, S.-I. (1999). Thermal-mechanical
1083
effects of low-temperature plasticity (the Peierls mechanism) on the defor-
1084
mation of a viscoelastic shear zone. Earth Planet. Sci. Lett., 168 , 159–172.
1085
Karato, S. (1989). Grain growth kinetics in olivine aggregates. Tectono-
1086
1087
1088
physics, 168 , 255–273. Karato, S. (2008). Deformation of Earth Materials: An Introduction to the Rheology of the Solid Earth. New York: Cambridge University Press. 63
1089
1090
1091
Karato, S., & Wu, P. (1993). Rheology of the Upper Mantle: A Synthesis. Science, 260 , 771–778. Kaus, B. J. P., & Podladchikov, Y. Y. (2006).
Initiation of localized
1092
shear zones in viscoelastoplastic rocks. J. Geophy. Res. Solid Earth, 111 ,
1093
B04412.
1094
1095
Kelemen, P. B., & Hirth, G. (2007). A periodic shear-heating mechanism for intermediate-depth earthquakes in the mantle. Nature, 446 , 787–790.
1096
Kidder, S., & Ducea, M. N. (2006). High temperatures and inverted meta-
1097
morphism in the schist of Sierra de Salinas, California. Earth Planet. Sci.
1098
Lett., 241 , 422–437.
1099
Kohlstedt, D., Evans, B., & Mackwell, S. (1995). Strength of the lithosphere:
1100
Constraints imposed by laboratory experiments. J. Geophys. Res., 100 ,
1101
17,587–17,602.
1102
1103
Kohn, M. J. (2014). Himalayan Metamorphism and Its Tectonic Implications. Annual Review of Earth and Planetary Sciences, 42 , 381–419.
1104
Landuyt, W., & Bercovici, D. (2009). Variations in planetary convection via
1105
the effect of climate on damage. Earth and Planetary Science Letters, 277 ,
1106
29–37.
1107
Landuyt, W., Bercovici, D., & Ricard, Y. (2008). Plate generation and two-
1108
phase damage theory in a model of mantle convection. Geophys. J. Int.,
1109
174 , 1065–1080.
64
1110
Leloup, P. H., Arnaud, N., Lacassin, R., Kienast, J. R., Harrison, T. M.,
1111
Trong, T. T. P., Replumaz, A., & Tapponnier, P. (2001). New constraints
1112
on the structure, thermochronology, and timing of the Ailao Shan-Red
1113
River shear zone, SE Asia. J. Geophys. Res., 106 , 6683–6732.
1114
Leloup, P. H., & Kienast, J.-R. (1993). High-temperature metamorphism in a
1115
major strike-slip shear zone: the Ailao Shan-Red River, People’s Republic
1116
of China. Earth Planet. Sci. Lett., 118 , 213–234.
1117
Linckens, J., Herwegh, M., & M¨ untener, O. (2015). Small quantity but large
1118
effecthow minor phases control strain localization in upper mantle shear
1119
zones. Tectonophysics, 643 , 26–43.
1120
Linckens, J., Herwegh, M., M¨ untener, O., & Mercolli, I. (2011). Evolution
1121
of a polymineralic mantle shear zone and the role of second phases in the
1122
localization of deformation. J. Geophys. Res. Solid Earth, 116 , B06210.
1123
Mont´esi, L., & Hirth, G. (2003). Grain size evolution and the rheology of
1124
ductile shear zones: From laboratory experiments to postseismic creep.
1125
Earth Planet. Sci. Lett., 211 , 97–110.
1126
Mont´esi, L., & Zuber, M. (2002).
A unified description of localiza-
1127
tion for application to largescale tectonics.
1128
10.1029/2001JB000465.
J. Geophys. Res., 107 ,
1129
Mont´esi, L. G. J. (2007). A constitutive model for layer development in shear
1130
zones near the brittle-ductile transition. Geophys. Res. Lett., 34 , L08307.
1131
Mont´esi, L. G. J. (2013). Fabric development as the key for forming ductile 65
1132
shear zones and enabling plate tectonics. Journal of Structural Geology,
1133
50 , 254–266.
1134
Moresi, L., & Solomatov, V. (1998). Mantle convection with a brittle litho-
1135
sphere: Thoughts on the global tectonic style of the earth and venus.
1136
Geophys. J. Int., 133 , 669–682.
1137
Mulyukova, E., & Bercovici, D. (2017). Formation of lithospheric shear zones:
1138
Effect of temperature on two-phase grain damage. Phys. Earth Planet.
1139
Inter., 270 , 195 – 212.
1140
Mulyukova, E., & Bercovici, D. (2018). Collapse of passive margins by litho-
1141
spheric damage and plunging grain size. Earth Planet. Sci. Lett., 484 ,
1142
341–352.
1143
1144
1145
1146
Platt, J. P. (2015). Rheology of two-phase systems: a microphysical and observational approach. J. Struct. Geol., 77 , 213–227. Poirier, J. (1980). Shear localization and shear instability in materials in the ductile field. J. Struct. Geol., 2 , 135–142.
1147
Regenauer-Lieb, K., & Yuen, D. A. (1998). Rapid conversion of elastic en-
1148
ergy into plastic shear heating during incipient necking of the lithosphere.
1149
Geophys. Res. Lett., 25 , 2737–2740.
1150
Rozel, A., Ricard, Y., & Bercovici, D. (2011). A thermodynamically self-
1151
consistent damage equation for grain size evolution during dynamic re-
1152
crystallization. Geophys. J. Intl., 184 , 719–728.
66
1153
1154
1155
1156
1157
1158
Schiesser, W. E. (2012). The numerical method of lines: Integration of partial differential equations. Elsevier. Schubert, G., & Yuen, D. A. (1978). Shear heating instability in the earth’s upper mantle. Tectonophysics, 50 , 197–205. Skemer, P., Warren, J. M., & Kelemen, P. B. (2010). Microstructural and rheological evolution of a mantle shear zone. J. Petrol., 51 , 43–53.
1159
Tackley, P. (2000a). The quest for self-consistent generation of plate tec-
1160
tonics in mantle convection models. In M. A. Richards, R. Gordon, &
1161
R. van der Hilst (Eds.), History and Dynamics of Global Plate Motions,
1162
Geophys. Monogr. Ser. (pp. 47–72). Washington, DC: Am. Geophys. Union
1163
volume 121.
1164
Tackley, P. (2000b). Self-consistent generation of tectonic plates in time-
1165
dependent, three-dimensional mantle convection simulations, 1. Pseudo-
1166
plastic yielding. Geochem. Geophys. Geosyst. (G3 ), 1 .
1167
Thielmann, M. (2017). Grain size assisted thermal runaway as a nucleation
1168
mechanism for continental mantle earthquakes: Impact of complex rheolo-
1169
gies. Tectonophysics, .
1170
Thielmann, M., & Kaus, B. J. P. (2012). Shear heating induced lithospheric-
1171
scale localization: Does it result in subduction? Earth Planet. Sci. Lett.,
1172
359 , 1–13.
1173
Thielmann, M., Rozel, A., Kaus, B. J. P., & Ricard, Y. (2015). Intermediate-
1174
depth earthquake generation and shear zone formation caused by grain size
1175
reduction and shear heating. Geology, 43 , 791–794. 67
1176
Tommasi, A., Knoll, M., Vauchez, A., Signorelli, J. W., Thoraval, C., &
1177
Log´e, R. (2009). Structural reactivation in plate tectonics controlled by
1178
olivine crystal anisotropy. Nature Geoscience, 2 , 423–427.
1179
1180
1181
1182
Turcotte, D., & Schubert, G. (2002). Geodynamics. (2nd ed.). New York: Cambridge Univ. Press. van Heck, H., & Tackley, P. (2008). Planforms of self-consistently generated plates in 3d spherical geometry. Geophys. Res. Lett., 35 , L19312.
1183
Warren, J. M., & Hirth, G. (2006). Grain size sensitive deformation mech-
1184
anisms in naturally deformed peridotites. Earth Planet. Sci. Letts., 248 ,
1185
438 – 450.
1186
1187
White, S., Burrows, S., Carreras, J., Shaw, N., & Humphreys, F. (1980). On mylonites in ductile shear zones. J. Struct. Geol., 2 , 175–187.
1188
Yuen, D. A., Fleitout, L., Schubert, G., & Froidevaux, C. (1978). Shear defor-
1189
mation zones along major transform faults and subducting slabs. Geophys.
1190
J. Int., 54 , 93–119.
1191
Yuen, D. A., & Schubert, G. (1979). On the stability of frictionally heated
1192
shear flows in the asthenosphere. Geophysical Journal International , 57 ,
1193
189–207.
68