A mantle conveyor belt beneath the Tethyan collisional belt

Thorsten W. Becker ∗ Department of Earth Sciences, University of Southern California, Los Angeles, CA

Claudio Faccenna Dipartimento Scienze Geologiche, Universit`a Roma TRE and IGAG, CNR Rome

Abstract

Collisional belts are generated by the arrival of continental lithosphere into a subduction zone, leading to stacking of crustal slices during indentation. The Tethyan suture from the Bitlis to the Himalayas is a prime example where the Arabian and Indian plates collided with Eurasia during the Cenozoic, generating the highest mountain belts on Earth (Argand, 1924). While the kinematics of this process are well established, its dynamics are more uncertain. India and Arabia intriguingly keep advancing in spite of large collisional resisting forces. We perform global mantle circulation computations to test the role of deep mantle flow as a driving force for the kinematics of the Tethyan collisional belt, evaluating different boundary conditions and mantle density distributions as inferred from seismic tomography or slab models. Our results show that mantle drag exerted on the base of the lithosphere by a large-scale upwelling is likely the main cause for the ongoing indentation of the Indian and Arabian plates into Eurasia. Key words: continental collision; mantle upwelling; plate motions; plate driving forces

Preprint submitted to Elsevier

26 April 2011

1

1 Introduction

2

Collisional orogens are among the most impressive manifestations of plate tecton-

3

ics, generating high relief by piling up of buoyant crustal slivers scraped from the

4

downgoing plate while trenches advance into the upper continental plate. The clas-

5

sical theory for orogeny suggests that crustal thickening persists after the entrance

6

of continental lithosphere at the trench and ceases when the driving pull of the

7

subducting oceanic lithosphere vanishes, e.g. after detachment of the subducting

8

lithosphere (Cloos, 1993). However, along the Himalayas, collision was sustained

9

over most of the Cenozoic, pushing the Indian plate into Asia for more than 600 km,

10

forming the Tibetan plateau and the Tien Shan (Yin and Harrison, 2000; Johnson,

11

2002; Guillot et al., 2003). On a smaller scale and at a reduced rate, Arabia also ad-

12

vanced into Eurasia, creating the Bitlis-Zagros collisional belt, during the spreading

13

of the Red Sea-Gulf of Aden ocean (McQuarrie et al., 2003; Hatzfeld and Molnar,

14

2010).

15

A longstanding, but not entirely resolved question then concerns the driving forces

16

for the northward motion of Arabia and India, which has to be sufficient to over-

17

come the large resistance of the collisional system for a long time span.

∗ Address: Department of Earth Sciences, University of Southern California, MC 0740, 3651 Trousdale Pkwy, Los Angeles, CA 90089–0740, USA. Phone: ++1 (213) 740 8365, Fax: ++1 (213) 740 8801 Email address: [email protected] (Thorsten W. Becker).

2

18

1.1

Plate kinematics and dynamic models

19

Despite the differences in size and shape, the present kinematic patterns associated

20

with the motion of the India and Arabia plates with respect to fixed Eurasia shows

21

striking similarities (Fig. 1A; Zhang et al., 2004; Gan et al., 2007; ArRajehi et al.,

22

2010). Separated via ridge systems from the large and stable African continent, the

23

two plates show northward indentation into Eurasia. On both plates, the geodeti-

24

cally inferred velocities decrease gradually inside the collision zone, progressively

25

turning laterally outward, toward active subduction zones: the Hellenic and Java-

26

Sumatra trenches for Arabia and India plate, respectively.

27

The geologic history of separation of the two plates also shares common fea-

28

tures. Both India and Arabia split from a large, continental plate after the arrival

29

of mantle plumes, as indicated by the emplacement of Large Igneous Provinces

30

(LIPs), followed by seafloor spreading. India separated from Antarctica-Australia

31

at ∼125 Ma, associated with the emplacement of the Kerguelen LIP (Gaina et al.,

32

2007), from Madagascar after the emplacement of the Morondova LIP (91-84 Ma;

33

Torsvik et al., 2000), and from the Seychelles micro continent just after the em-

34

placement of the Deccan LIP (∼65 Ma), related to the Reunion plume. During this

35

period, India’s plate motion increased up to the highest reliably recorded velocities

36

(∼16-18 cm/yr; e.g. Patriat and Achache, 1984; Copley et al., 2010). It is typical-

37

ly inferred that the continental lithosphere of India collided with Asia at ∼50 Ma

38

(Molnar and Tapponnier, 1975; Guillot et al., 2003; Hatzfeld and Molnar, 2010),

39

but an alternative view is that collision initiated later in the Tertiary (Aitchison

40

et al., 2007). After collision, India’s plate velocity decreased, to ∼4-5 cm/yr at the

41

present day (Zhang et al., 2004; Copley et al., 2010). Approximately half of the

42

plate velocity during the last ∼25 Ma is accommodated by overriding plate thick3

43

ening and extrusion, with the corresponding subduction velocity being reduced to

44

. 2 cm/yr. This process is associated with shallow dipping under-thrusting of In-

45

dia’s continental lithosphere below Tibet (Capitanio et al., 2010; Replumaz et al.,

46

2010).

47

Arabia broke apart from Africa at ∼35 Ma, after the arrival of the Afar plume and

48

spreading of the Aden Gulf-Red Sea (Ebinger and Sleep, 1998; McQuarrie et al.,

49

2003). Continental under-thrusting probably initiated at ∼30 Ma while final colli-

50

sion occurred ∼12 Ma with a convergence rate of ∼ 2.5 cm/yr (Jolivet and Faccen-

51

na, 2000; Allen and Armstrong, 2008; Hatzfeld and Molnar, 2010). For the case of

52

Arabia, it appears that its translation occurred during continental under-thrusting

53

(McQuarrie et al., 2003; Hatzfeld and Molnar, 2010). Similar to the case of In-

54

dia, about half of the present-day shortening is accommodated within the Zagros

55

belt for the Arabia-Eurasia collision (McClusky et al., 2003; ArRajehi et al., 2010).

56

Hence, the subduction velocity corrected for under-thrusting, is likely . 1.5 cm/yr.

57

We interpret the geologic record for the collisional system such that both conver-

58

gence and subduction velocities decreased rapidly after continental collision, but

59

then remained fairly constant over the last 20-30 Myrs. By inference, the present-

60

day kinematics, for example as recorded by geodesy, may then be representative of

61

the long-term force equilibrium (Meijer and Wortel, 1999).

62

Different classes of models have been proposed to understand the dynamics of the

63

Tethyan collisional system. The first, and perhaps most popular, class of models

64

proposes that the negative buoyancy force exerted by subduction of oceanic litho-

65

sphere (slab pull) may propel India and Arabia against Eurasia. This mechanism

66

might also work for the case of continental subduction, provided that the buoy-

67

ant, upper-crustal layers are scraped off from the lithospheric mantle (Cloos, 1993;

68

Capitanio et al., 2010). Laboratory and numerical models have reproduced the re4

69

duction of the subducting and convergence velocity by ratios similar to what is

70

observed for India (Bellahsen et al., 2003; Capitanio et al., 2010). Slab pull was

71

most likely a major driving force during India’s fast drift, modulating the decrease

72

in the convergence motion (Capitanio et al., 2010; van Hinsbergen et al., 2011).

73

However, at present slab pull is expected to be reduced to its minimum level, giv-

74

en the inferred repeated episodes of slab break off (Schott and Schmeling, 1998;

75

Chemenda et al., 2000; Keskin, 2003; Faccenna et al., 2006; Replumaz et al., 2010).

76

Another important contribution is represented by the lithospheric thickening (ridge

77

push) force of the oceanic lithosphere. Together with the pull force of the Tethyan

78

subducting slab, ridge push has been considered as driving the indentation (Copley

79

et al., 2010; Capitanio et al., 2010). Recent estimates of ridge push forces based on

80

plate structure and on force equilibrium give values between 1 and 2 · 1012 N/m,

81

however, too low to explain the deformation associated with the collisional system

82

(Ghosh et al., 2006). Intraplate forcing transmitted from the surrounding plates,

83

especially Indochina and Australia, could also contribute to explain the motion of

84

the colliding plates (Cloetingh and Wortel, 1986; Coblentz et al., 1998; Meijer and

85

Wortel, 1999; Flesch et al., 2005; Li et al., 2008b) given that the boundary between

86

Australia and India is not well defined.

87

Lastly, a potential contribution is represented by the drag exerted by large scale

88

mantle flow, a “continental undertow” (Alvarez, 1990). Mantle drag associated with

89

regional, plume-like upwellings has been suggested as an efficient mechanism for

90

rapid drift of continental plates, such as Laurentia or Baltica, when velocities are

91

in excess of 10 cm/yr (Gurnis and Torsvik, 1994; van Hinsbergen et al., 2011). The

92

sustained convergence of the Tethyan belt may then be associated with a long-term

93

drag exerted by larger-scale mantle flow (Alvarez, 2010; van Hinsbergen et al.,

94

2011). 5

95

Here, we use global mantle flow computations to test the influence of different

96

models of density anomalies within the mantle and/or the subducting lithosphere

97

on present-day plate motions. Using different boundary conditions, we are able to

98

investigate the mutual role of mantle drag, plate interactions, and the subducting

99

slab on the kinematics of the India and Arabia collisional system.

100

2 Methods

101

To examine the mantle drag contributions due to different buoyancy force distri-

102

butions, we use instantaneous, 3-D spherical mantle flow models (Hager and O’-

103

Connell, 1981). To model mantle flow, we solve the infinite Prandtl number, Stokes

104

equation for incompressible, instantaneous fluid flow (Boussinesq approximation)

105

in a global, spherical shell using the global, finite-element code CitcomS (Zhong

106

et al., 2000), with our own modifications as shared through CIG (geodynamics.org).

107

We employ a numerical resolution of ∼20 km laterally in the upper mantle which

108

we tested was sufficient to resolve flow behavior of the lithosphere on the 100 km

109

scales we are interested in (Faccenna and Becker, 2010). Our parameter choic-

110

es are motivated by previous work on fitting plate motions and the geoid (Ghosh

111

et al., 2010). The rheology of the mantle is Newtonian viscous and most models

112

only have radial viscosity variations besides weak zones which are confined to the

113

lithosphere (down to 100 km depth), have 100 km width, and a viscosity reduction

114

to 0.01 the ambient viscosity. The reference, radial viscosity structure is 5 · 1022 Pas

115

in the lithosphere, 1021 Pas from 100-660 km, and 5 · 1022 Pas in the lower mantle

116

(see supplementary material).

117

The solution for the instantaneous mantle flow field relies on assuming that densi-

118

ties are known within the computational domain. Density models may be construct6

119

ed by scaling the velocity anomalies from seismic tomography to temperature (e.g.

120

Hager et al., 1985), by inferring density for slabs alone, based on seismicity in

121

Wadati-Benioff zones and/or inferred past subduction (e.g. Hager, 1984; Ricard

122

et al., 1993; Lithgow-Bertelloni and Richards, 1998), or by mixed models that ex-

123

plore the respective role of slabs versus tomography and edge forces (e.g. Becker

124

and O’Connell, 2001; Conrad and Lithgow-Bertelloni, 2002; Ghosh et al., 2010;

125

Stadler et al., 2010). The density structure in our reference model is based on the

126

composite S wave tomography model SMEAN (Becker and Boschi, 2002). We re-

127

move velocity anomalies underneath cratonic keels (from 3SMAC by Nataf and Ri-

128

card, 1996) down to 250 km to avoid being affected by presumably compositional

129

anomalies there (cf. discussion in Forte, 2007). At a Rayleigh number of 3.4 · 108

130

(definition of Zhong et al., 2000), non-dimensional temperatures are scaled such

131

that density, ρ, anomalies go as

132

133

134

d ln ρ = 0.2 d lnVS

(1)

for S wave tomography, and d ln ρ = 0.4 d lnVP

(2)

135

for P wave models, consistent with previous global circulation modeling (Becker

136

and O’Connell, 2001). We adopt such a scaling for simplicity, realizing that com-

137

positional anomalies and depth-dependent mineral physics relationships will com-

138

plicate the mapping between velocity anomalies and density (e.g. Steinberger and

139

Calderwood, 2006; Simmons et al., 2010). We assume that, by taking into account

140

cratonic roots, all remaining anomalies in the upper mantle are of thermal nature.

141

In addition to tomography-only density, we consider upper mantle models that are

142

entirely or partially based on slab structure inferred from the Wadati-Benioff zone 7

143

geometry, based on the RUM model (Gudmundsson and Sambridge, 1998), as in

144

Ghosh et al. (2010).

145

Our approach to model plate motions follows Ricard and Vigny (1989) and Zhong

146

et al. (2000) and is described in detail elsewhere (Becker, 2006; Faccenna and Beck-

147

er, 2010). It is geared toward exploring complex collisional systems by allowing

148

for different sets of surface boundary conditions, and by exploring the role of litho-

149

spheric weak zones in models that allow for lateral variations in viscosity (cf. King

150

et al., 1992). Surface kinematic boundary conditions for most models are shear-

151

stress free (i.e. the surface is free slip, allowing for dynamically-consistent plate

152

motions), with weak zones following plate boundaries. The intraplate weak zones

153

in Asia (Fig. 2A) focus deformation in India and a southeast corner of the Eurasian

154

plate. Results are presented showing surface motions around the Tethyan with re-

155

spect to a best-fit, fixed Eurasian plate reference frame (white vectors in Fig. 2B,

156

C), allowing for significant intraplate deformation. We also run models to test the

157

influence of single or overall plate motion on the kinematics of the colliding system

158

by prescribing velocities on part of the surface, while leaving the rest shear-stress

159

free, as in Faccenna and Becker (2010).

160

The normal stresses generated by viscous flow at the surface are used to infer the

161

equivalent “dynamic” topography (Ricard et al., 1984; Richards and Hager, 1984),

162

which can be compared to observed topography once the effect of isostatic ad-

163

justment is removed to arrive at “residual topography” (Fig. 1B). Model results

164

are then compared with the present-day geodetic velocity fields and with geomor-

165

phological information concerning the vertical motion (e.g. Gurnis et al., 2000;

166

Daradich et al., 2003). 8

167

3 Results

168

Fig. 2C shows surface velocities for our reference model. The motion of the Indi-

169

an plate is fairly well reproduced compared to geodetic velocities (orange arrows,

170

averaged from Fig. 1A). The direction of the motion of Arabia is also matched,

171

although its rate is overestimated. The Pacific and Philippine Sea plate are fit by

172

the model as well, with an under-predicted rate. The motion of the Australian plate,

173

conversely, is not correctly reproduced, but as shown next, this is not significantly

174

affecting the kinematics of the colliding system. We therefore refrain from optimiz-

175

ing our models for global motion fit (cf. Forte, 2007), for simplicity.

176

In detail, the motion of India shows a progressive decrease in velocity moving

177

from the Carlsberg ridge to the collisional zone (Fig. 2C). This is seen in sections

178

across the Indian and Arabia plates, showing an upwelling in the upper and low-

179

er mantle (Figure 2D), pushing plates toward the collisional zone, where a high

180

velocity anomaly –a trace of an older subduction zone– (Figs. 2E and F) rep-

181

resents the return flow into the upper and lower mantle. Overall, the upwelling-

182

downwelling circuit form a convection “cell” (“conveyor belt”, realizing that there

183

is three-dimensionality to the flow), encompassing the motion of the India plate,

184

with northward flow concentrated in the upper mantle and the south-directed re-

185

turn flow in the lower mantle. The low velocity anomalies associated with this flow

186

are found beneath Ethiopia-Afar, spreading northward towards Arabia and beneath

187

Reunion spreading northward beneath the Carlsberg ridge.

188

The convective pattern associated with the conveyor belt produces a multi-scale

189

wavelength topographic signal that, to first order, visually matches many features

190

of the residual topography (cf. Figs. 1B and 2C). This indicates that at least part of 9

191

the topography of the Tethyan and surrounding regions may be dynamically sup-

192

ported by the mantle. A previously recognized example is the prominent positive

193

topography over east Africa-eastern Arabia. This produces an overall uplift and tilt-

194

ing of the Arabian plate (Daradich et al., 2003), already attributed to the large scale

195

low-velocity mantle anomaly mapped under Africa (Lithgow-Bertelloni and Silver,

196

1998; Gurnis et al., 2000). This large-scale topographic signal is rooted toward the

197

south, on what is presumably a hotspot acting over the last 30 Ma, first over east

198

Africa and presently beneath Afar. The model matches the residual topography es-

199

timate (Fig. 1B), showing a broad upwelling extending north toward the Middle

200

East (cf. Boschi et al., 2010).

201

Positive dynamic topography also extends over the presently active hotspot beneath

202

Reunion, the origin of the Deccan trap eruption at the end of the Cretaceous. This

203

feature extends northward along the Carlsberg and Central Indian Ridges and al-

204

so connects to the south-Africa lower mantle, low velocity anomaly (Gurnis et al.,

205

2000; Steinberger et al., 2001). In our model, this upwelling is responsible for a

206

large part of the driving force of the India motion. Other high dynamic topography

207

features are found in the north China Sea and in the Baikal region. The Hainan vol-

208

cano lies on top of the broad upwelling zone in south China which has been related

209

to plumes forming on top of a prominent low-velocity anomaly in the mantle (Lei

210

et al., 2009). Doming, volcanism and uplift beneath the Baikal-Mongolia plateau

211

has been associated with lithosphere plume interactions (Windley and Allen, 1993).

212

The Tethyan belt is only locally (Caucasus-Anatolia-Iran) marked by a positive to-

213

pography signal. Comparison with the residual topography map shows that most

214

of the elevation is compensated by crustal thickening, predicting small, if any, dy-

215

namic uplift in the Tibetan plateau for this particular model.

216

To isolate the contribution of mantle flow on the micro plate motion, we set to ze10

217

ro the motion of all plates outside the Eurasian/Arabian domain (Fig. 3A). This

218

model shows that the influence of the surrounding plate motions is to turn the ve-

219

locity field northerly for India and Arabia by ∼ 10◦ -20◦ , but does not affect the

220

velocity amplitude significantly. Prescribing the plate motion of Australia or Africa

221

(Figs. 3B and C) to conform to NUVEL-1A (DeMets et al., 1994) does also not

222

produce a significant deviations from our reference model (compare Figs. 3C and

223

D with Fig. 2C). This implies that plate interactions produce only a moderate in-

224

fluence on the orientation of the colliding motion, indicating the primary role of

225

driving mantle tractions associated with the conveyor belt.

226

The individual role of the upwelling and downwelling (subduction zone) on the

227

plate motions can be inferred by comparing our reference model, which is based on

228

SMEAN, with one where density is inferred from a recent P-wave model (Li et al.,

229

2008a) (MIT08). By virtue of the construction with its reliance on body waves and

230

the specific datasets, the MIT P model has superior resolution close to slabs, and

231

underneath parts of Asia (Li et al., 2008a). The high velocity zone in the shallow

232

layer in the SMEAN model covers a large part of the collisional zone, whereas in

233

MIT08 is restricted to the subduction zones (Supplementary Fig. S1). More impor-

234

tantly, the MIT08 model does not show the low velocity anomaly beneath Reunion

235

and Carlsberg ridge, presumably because of lack of coverage (SMEAN includes in-

236

formation from surface waves). As a result, India is not moving correctly (Figs. 3B

237

and S2).

238

The role of subduction zones in driving plate motions can also be tested by assign-

239

ing density anomalies to regions with Wadati-Benioff seismicity. In our study re-

240

gion, deep seismicity is distributed along the Pacific margins, along Java-Sumatra-

241

Burma, beneath the Hindu Kush cluster, and, far to the west, along the Hellenic

242

trench (Fig. 1A). The flow predictions for slabs only (Fig. 3E) indicate that sub11

243

duction does drive Australia northward, although at a significantly reduce rate. In-

244

dia and Arabia move slowly, if at all, with Arabia-Anatolia mainly pulled to the

245

west from the Hellenic slab. The northward motion of India and Arabia increase

246

by adding lower mantle anomalies and, even more, by adding upper mantle slow

247

velocity anomalies (Fig. 3F). In addition, this latter model shows a good fit with

248

crustal velocities. We take this to indicate the role of active subduction zones on

249

the side of the colliding system, which leads to an important regional modifica-

250

tion of the overall convergence kinematics that are themselves mainly driven by

251

a large-scale mantle upwelling. This is particularly evident on the Mediterranean

252

side, where the Hellenic slab is pulling Anatolia westward (cf. Boschi et al., 2010).

253

4 Discussion

254

Given our results, we can now infer the most important driving forces for the

255

Tethyan collisional system at present-day by elimination. The pull exerted by sub-

256

ducting lithosphere is here modeled by density anomalies either restricted to the

257

high velocity regions in the upper mantle as imaged by tomography, or to Wadati

258

Benioff zones. Models with well resolved high velocities in the subduction zone,

259

but poorly resolved low velocity anomalies in upper mantle regions, such as under-

260

neath the Carlsberg ridge, produce no motion of India (Figs. 3A and S2). Similarly,

261

given that Wadati-Benioff seismicity is almost absent beneath the colliding area,

262

the restriction of the density anomaly to the seismically active slabs, expectedly,

263

does not lead to any significant motion of the colliding plates.

264

While the lack of a direct link between slab-related driving forces and India and

265

Arabia plate motions in our models is not entirely representative of all aspects of

266

the slab pull mechanism, we infer that the role of subducting slabs on the present12

267

day velocity field cannot be dominant. Geological evidence indicates that the Hi-

268

malayan orogeny has been punctuated by episodes of slab rupture, the first probably

269

occurring soon after collision (Chemenda et al., 2000; Replumaz et al., 2010). This

270

is similar to what has been proposed for Arabia (Keskin, 2003; Faccenna et al.,

271

2006), and would limit the possible action of slab-pull forces. Recent tomographic

272

images illustrate that the slabs beneath India and Arabia (e.g. Keskin, 2003; Li

273

et al., 2008a) are now rather small, perhaps under-thrusting sub-horizontally, and

274

shallow slab segments appear separated from deeper slabs, which are either pond-

275

ing above the 660 km phase transition, or in the lower mantle behind the present-

276

day subduction zone (van der Voo et al., 1999; Hafkenscheid et al., 2006).

277

The potential energy forces related to lithospheric thickening and continental struc-

278

ture are here not taken into account directly. However, including the first 100 km

279

depth velocity anomalies from tomography means including most of the half-space

280

cooling signature of oceanic lithosphere and hence ridge push (Fig. 3F). A com-

281

parison with a model without that layer (not shown) indicates that this contribution

282

is, however, not predominant. Ghosh et al. (2006) estimate the contribution of the

283

ridge push in the Indian ocean and conclude that the gravitational potential ener-

284

gy exerted by the Tibetan plateau (∼ 3 · 1012 N/m) exceeds the ridge push level.

285

Consistent with our findings, Ghosh et al. also point out that an additional force,

286

presumably mantle drag, is required to fit the present-day stress field of the region.

287

For the case of Arabia, in addition, the ridge push contribution cannot be large

288

because the Gulf of Aden-Red Sea is only a narrow, incipient ocean. In terms of

289

plate interactions, the contribution of the northward advancing Australian plate is

290

relevant in determining the direction of motion of India, but appears to not affect

291

the amplitudes of India motion strongly. This is shown by comparing our reference

292

model with the one where Australian motion is prescribed (Figs. 2C and 3C). 13

293

This leads us to suggest that the most important force contribution is related to

294

mantle flow. The reference model (Fig. 2), and a range of tests for different tomo-

295

graphic models and rheology (including Figs. S3 and S4), consistently show that an

296

important contribution for plate motion in the system is derived from the large-scale

297

convection cell beneath the two colliding plates (Fig. 4). In particular, the compari-

298

son between different tomography models shows that upwelling, upper mantle flow

299

beneath the Reunion-Carlsberg ridge is, and probably has been, a fundamental in-

300

gredient for the collisional history of India. Comparing the reference model with

301

Fig. 3F also shows that incorporating the density distribution within upper mantle

302

subduction zones in addition to active upwellings gives a better match of the over-

303

all velocity field. For example, implementing density anomaly within the Hellenic

304

trench, promote the westward motion of the Anatolia-Arabia system (cf. Ghosh

305

et al., 2010; Stadler et al., 2010).

306

All of our main models fail to reproduce the toroidal flow trajectories that are pro-

307

nounced on the eastern side of the India indenter, where motions turn to SE ori-

308

entations with respect to Eurasia. We deduce that this flow component is related

309

to lower crustal, or decoupled lithospheric, flow underneath the Yunnan domain as

310

suggested by Royden et al. (1997). Material appears to escape laterally from the

311

convergence zone under the gravitational potential energy of the high plateau with

312

respect to the surrounding lowlands (cf. Royden, 2008). This motion is perhaps

313

assisted by uppermost mantle upwellings (Fig. 3A) or slab suction (Fig. 3E), but

314

these effects are swamped in our models by the overall convergent flow, and rhe-

315

ologically more complex models with regional decoupling may be required (Clark

316

et al., 2005; Flesch et al., 2005; Sol et al., 2007). 14

317

5 Conclusions

318

We substantiate previous suggestions that mantle drag is highly efficient in driving

319

continental plates (Gurnis and Torsvik, 1994; Alvarez, 2010; van Hinsbergen et al.,

320

2011). Our models indicate that a plume-related, “conveyor belt” represents the pri-

321

mary driving force for India and Arabia to collide, indenting against the Eurasian

322

plate and supporting the thick Tibetan and Iranian mountain ranges. We posit that

323

particularly the upwelling component of the conveyor belt is a significant driving

324

force, breaking continental plates into small pieces and dragging them away against

325

the upper plate. The upwelling extends from East Africa to the Indian Ocean, and

326

represents the shallow, upper mantle expression of the deep south African low-

327

velocity zone. Our models suggest that the gravitational pull induced by subduc-

328

tion represents a secondary mechanism, accommodating convergence and leaving

329

behind drips of high velocity material in the mantle beneath the advancing colli-

330

sional zone.

331

Acknowledgments 332

333

We thank A. Ghosh for assistance with assembling mantle density models. Com-

334

putations were performed on USC’s High Performance Computing Center, and we

335

thank CIG and seismologists sharing their tomographic models in electronic form.

336

This research was supported by NSF grants EAR 0930046 and 0643365.

337

References

338

Aitchison, J. C., Ali, J. R., Davis, A. M., 2007. When and where did India and Asia

339

collide? J. Geophys. Res. 112, B05423. doi:10.1029/2006JB004706. 15

340

Allen, M. B., Armstrong, H. A., 2008. ArabiaEurasia collision and the forcing of

341

mid-Cenozoic global cooling. Palaeogeogr., Palaeoclim., Palaeoeco. 265, 3152–

342

3158.

343

Alvarez, W., 1990. Geologic evidence for the plate-driving mechanism: the conti-

344

nental undertow hypothesis and the Australian-Antarctic Discordance. Tectonics

345

9, 1213–1220.

346

347

348

349

Alvarez, W., 2010. Protracted continental collisions argue for continental plates driven by basal traction. Earth Planet. Sci. Lett. 296, 434–442. Argand, E., 1924. La tectonique de l’asie. In: Compte-rendu du 13e congr`es G´eologique International. Brussels, pp. 171–372.

350

ArRajehi, A., McClusky, S., Reilinger, R., Daoud, M., Alchalbi, A., Ergintav, S.,

351

Gomez, F., Sholan, J., Bou-Rabee, F., Ogubazghi, G., Haileab, B., Fisseha, S.,

352

Asfaw, L., Mahmoud, S., Rayan, A., Bendik, R., Kogan, L., 2010. Geodetic

353

constraints on present day motion of the Arabian Plate: Implications for Red Sea

354

and Gulf of Aden rifting. Tectonics 29, TC3011. doi:10.1029/2009TC002482.

355

356

Bassin, C., Laske, G., Masters, G., 2000. The current limits of resolution for surface wave tomography in North America (abstract). Eos Trans. AGU 81, F897.

357

Becker, T. W., 2006. On the effect of temperature and strain-rate dependent viscos-

358

ity on global mantle flow, net rotation, and plate-driving forces. Geophys. J. Int.

359

167, 943–957.

360

Becker, T. W., Boschi, L., 2002.

361

namic mantle models.

362

2001GC000168.

363

364

A comparison of tomographic and geody-

Geochem., Geophys., Geosys. 3, 1.

doi:10.1029/

Becker, T. W., O’Connell, R. J., 2001. Predicting plate velocities with geodynamic models. Geochem., Geophys., Geosys. 2, 12. doi:10.1029/2001GC000171.

365

Bellahsen, N., Faccenna, C., R., F., J.M., D., Jolivet, L., 2003. Why did Africa

366

separate from Africa? Insight from 3D laboratory experiments. Earth Planet. 16

367

Sci. Lett. 216, 365–381.

368

Boschi, L., Faccenna, C., Becker, T. W., 2010. Mantle structure and dynamic to-

369

pography in the Mediterranean Basin. Geophys. Res. Lett. 37, L20303. doi:

370

10.1029/2010GL045001.

371

Capitanio, F. A., Morra, G., Goes, S., Weinberg, R. F., Moresi, L., 2010. India-Asia

372

convergence driven by the subduction of the Greater Indian continent. Nature

373

Geosciences 3, 136–139.

374

Chemenda, A. I., Burg, J. P., Mattauer, M., 2000. Evolutionary model of the

375

Himalaya-Tibet system: geopoem based on new modelling, geological and geo-

376

physical data. Earth Planet. Sci. Lett. 174, 397–409.

377

Clark, M. K., Bush, J. W. M., Royden, L. H., 2005. Dynamic topography produced

378

by lower crustal flow against rheological strength heterogeneities bordering the

379

Tibetan Plateau. Geophys. J. Int. 162, 575–590.

380

381

Cloetingh, S., Wortel, M. J. R., 1986. Stress in the Indo-Australian plate. Tectonophys. 132, 49–67.

382

Cloos, M. N., 1993. Lithospheric buoyancy and collisional orogenesis: subduc-

383

tion of oceanic plateaus, continental margins, island arcs, spreading ridges and

384

seamounts. Geol. Soc. Am. Bull. 105, 715–737.

385

Coblentz, D. D., Zhou, S., Hillis, R. R., Richardson, R. M., Sandiford, M., 1998.

386

Topography, boundary forces, and the Indo-Australian intraplate stress field. J.

387

Geophys. Res. 103, 919–931.

388

389

Conrad, C. P., Lithgow-Bertelloni, C., 2002. How mantle slabs drive plate tectonics. Science 298, 207–209.

390

Copley, A., Avouac, J.-P., Royer, J.-Y., 2010. India-Asia collision and the Cenozoic

391

slowdown of the Indian plate: Implications for the forces driving plate motions.

392

J. Geophys. Res. 115, B03410. doi:10.1029/2009JB006634.

393

Daradich, A., Mitrovica, J. X., Pysklywec, R. N., Willet, S. D., Forte, A., 2003. 17

394

Mantle flow, dynamic topography, and rift-flank uplift of Arabia. Geology 31,

395

901–904.

396

DeMets, C., Gordon, R. G., Argus, D. F., Stein, S., 1994. Effect of recent revisions

397

to the geomagnetic reversal time scale on estimates of current plate motions.

398

Geophys. Res. Lett. 21, 2191–2194.

399

400

Ebinger, C. J., Sleep, N. H., 1998. Cenozoic magmatism throughout east Africa resulting from impact of a single plume. Nature 395, 788–791.

401

Engdahl, E. R., van der Hilst, R. D., Buland, R., 1998. Global teleseismic earth-

402

quake relocation with improved travel times and procedures for depth determi-

403

nation. Bull. Seismol. Soc. Am. 88, 722–743.

404

405

Faccenna, C., Becker, T. W., 2010. Shaping mobile belts by small-scale convection. Nature 465, 602–605.

406

Faccenna, C., Bellier, O., Martinod, J., Piromallo, C., Regard, V., 2006. Slab de-

407

tachment beneath eastern Anatolia: A possible cause for the formation of the

408

North Anatolian fault. Earth Planet. Sci. Lett. 242, 85–97.

409

Flesch, L. M., Holt, E. E., Silver, P. G., Stephenson, M., Wang, C.-Y., Chan, W. W.,

410

2005. Constraining the extent of crustmantle coupling in central Asia using GPS,

411

geologic, and shear wave splitting data. Earth Planet. Sci. Lett. 238, 248–268.

412

Forte, A. M., 2007. Constraints on seismic models from other disciplines – Impli-

413

cations for mantle dynamics and composition. In: Schubert, G., Bercovici, D.

414

(eds.), Treatise on Geophysics, Elsevier, Amsterdam. pp. 805–858.

415

Gaina, C., M¨uller, R., Brown, B., Ishihara, T., Ivanov, S., 2007. Breakup and early

416

seafloor spreading between India and Antarctica. Geophys. J. Int. 170, 151–170.

417

Gan, W., Zhang, P., Shen, Z., Niu, Z., Wang, M., Wan, Y., Zhou, D., Cheng, J.,

418

2007. Presentday crustal motion within the Tibetan Plateau inferred from GPS

419

measurements. J. Geophys. Res. 113, B08416. doi:10.1029/2005JB004120.

420

Ghosh, A., Becker, T. W., Zhong, S., 2010. Effects of lateral viscosity variations 18

421

on the geoid. Geophys. Res. Lett. 37, L01301. doi:10.1029/2009GL040426.

422

Ghosh, A., Holt, W. E., Flesch, L. M., Haines, A. J., 2006. Gravitational potential

423

energy of the Tibetan Plateau and the forces driving the Indian plate. Geology

424

34, 321–324.

425

426

Gudmundsson, O., Sambridge, M., 1998. A regionalized upper mantle (RUM) seismic model. J. Geophys. Res. 103, 7121–7136.

427

Guillot, S., Garzanti, E., Baratoux, D., Marquer, D., Mah´eo, G., de Sigoyer,

428

J., 2003. Reconstructing the total shortening history of the NW Himalaya.

429

Geochem., Geophys., Geosys. 4, 7. doi:10.1029/2002GC000484.

430

Gurnis, M., Mitrovica, J. X., Ritsema, J., van Heijst, H.-J., 2000. Constraining man-

431

tle density structure using geological evidence of surface uplift rates: The case of

432

the African superplume. Geochem., Geophys., Geosys. 1. doi:1999GC000035.

433

Gurnis, M., Torsvik, T., 1994. Rapid drift of large continents during the Late Pre-

434

cambrian and Paleozoic: Paleomagnetic constraints and dynamics models. Ge-

435

ology 22, 1023–1026.

436

Hafkenscheid, E., Wortel, M. J. R., Spakman, W., 2006. Subduction history of the

437

Tethyan region derived from seismic tomography and tectonic reconstructions.

438

J. Geophys. Res. 111, B08401. doi:10.1029/2005JB003791.

439

440

Hager, B. H., 1984. Subducted slabs and the geoid: constraints on mantle rheology and flow. J. Geophys. Res. 89, 6003–6015.

441

Hager, B. H., Clayton, R. W., Richards, M. A., Comer, R. P., Dziewo´nski, A. M.,

442

1985. Lower mantle heterogeneity, dynamic topography and the geoid. Nature

443

313, 541–545.

444

445

Hager, B. H., O’Connell, R. J., 1981. A simple global model of plate dynamics and mantle convection. J. Geophys. Res. 86, 4843–4867.

446

Hatzfeld, D., Molnar, P., 2010. Comparisons of the kinematics and deep structures

447

of the Zagros and Himalaya and of the Iranian and Tibetan plateaus and geody19

448

449

450

451

452

namic implications. Rev. Geophys. 48, RG2005. doi:10.1029/2009RG000304. Johnson, M. R. W., 2002. Shortening budgets and the role of continental subduction during the India-Asia collision. Earth Sci. Rev. 59, 101–123. Jolivet, L., Faccenna, C., 2000. Mediterranean extension and the Africa-Eurasia collision. Tectonics 6, 1095–1107.

453

Keskin, M., 2003. Magma generation by slab steepening and breakoff beneath

454

a subduction-accretion complex: an alternative model for collision-related vol-

455

canism in eastern Anatolia, Turkey.

456

10.1029/2003GL018019.

Geophys. Res. Lett. 30, 8046.

doi:

457

King, S. D., Gable, C. W., Weinstein, S. A., 1992. Models of convection-driven

458

tectonic plates: a comparison of methods and results. Geophys. J. Int. 109, 481–

459

487.

460

Lei, J., Zhao, D., Steinberger, B., Wu, B., Shen, F., Li, Z., 2009. New seismic con-

461

straints on the upper mantle structure of the Hainan plume. Phys. Earth Planet.

462

Inter. 33–50.

463

Li, C., van der Hilst, R. D., Engdahl, E. R., Burdick, S., 2008a. A new global model

464

for P wave speed variations in Earth’s mantle. Geochem., Geophys., Geosys. 9,

465

Q05018. doi:10.1029/2007GC001806.

466

Li, C., van der Hilst, R. D., Meltzer, A. S., Engdahl, E. R., 2008b. Subduction of

467

the Indian lithosphere beneath the Tibetan Plateau and Burma. Earth Planet. Sci.

468

Lett. 274, 157–168.

469

470

471

472

Lithgow-Bertelloni, C., Richards, M. A., 1998. The dynamics of Cenozoic and Mesozoic plate motions. Rev. Geophys. 36, 27–78. Lithgow-Bertelloni, C., Silver, P. G., 1998. Dynamic topography, plate driving forces and the African superswell. Nature 395, 269–272.

473

McClusky, S., Reilinger, R., Mahmoud, S., Sari, D. B., Tealeb, A., 2003. GPS

474

constraints on Africa (Nubia) and Arabia plate motions. Geophys. J. Int. 155, 20

475

126–138.

476

McQuarrie, N., Stock, J. M., Verdel, C., Wernicke, B. P., 2003. Cenozoic evolution

477

of Neotethys and implications for the causes of plate motions. Geophys. Res.

478

Lett. 30, 2036. doi:10.1029/2003GL017992.

479

480

481

482

Meijer, P. T., Wortel, M. J. R., 1999. Cenozoic dynamics of the African plate with emphasis on the Africa-Eurasia collision. J. Geophys. Res. 104, 7405–7418. Molnar, P., Tapponnier, P., 1975. Cenozoic tectonics of Asia: Effects of a continental collision. Science 189, 419–426.

483

Nataf, H.-C., Ricard, Y., 1996. 3SMAC: An a priori tomographic model of the

484

upper mantle based on geophysical modeling. Phys. Earth Planet. Inter. 95, 101–

485

122.

486

Patriat, P., Achache, J., 1984. India-Eurasia collision chronology has implications

487

for crustal shortening and driving mechanism of plates. Nature 311, 615–621.

488

Replumaz, A., Negredo, A., Guillot, S., Villase˜nor, A., 2010. Multiple episodes

489

of continental subduction during India/Asia convergence: Insight from seismic

490

tomography and tectonic reconstruction. Tectonophys. 483, 125–134.

491

492

Ricard, Y., Fleitout, L., Froidevaux, C., 1984. Geoid heights and lithospheric stresses for a dynamic Earth. Ann. Geophys. 2, 267–286.

493

Ricard, Y., Richards, M. A., Lithgow-Bertelloni, C., Le Stunff, Y., 1993. A geo-

494

dynamic model of mantle density heterogeneity. J. Geophys. Res. 98, 21,895–

495

21,909.

496

497

498

499

500

501

Ricard, Y., Vigny, C., 1989. Mantle dynamics with induced plate tectonics. J. Geophys. Res. 94, 17,543–17,559. Richards, M., Hager, B. H., 1984. Geoid anomalies in a dynamic Earth. J. Geophys. Res. 89, 5987–6002. Royden, L. H., 2008. The geological evolution of the Tibetan Plateau. Science 321, 1054–1058. 21

502

Royden, L. H., Burchfiel, B. C., King, B. C., Wang, E., Chen, Z., Shen, F., Liu, Y.,

503

1997. Surface deformation and lower crustal flow in eastern tibet. Science 276,

504

788–790.

505

506

Schott, B., Schmeling, H., 1998. Delamination and detachment of a lithospheric root. Tectonophys. 296, 225–247.

507

Simmons, N. A., Forte, A. M., Boschi, L., Grand, S. P., 2010. GyPSuM: A joint

508

tomographic model of mantle density and seismic wave speeds. J. Geophys. Res.

509

115, B12310. doi:10.1029/2010JB007631.

510

Sol, S., Meltzer, A., B¨urgmann, R., van der Hilst, R. D., Chen, Z., Koons, P. O.,

511

Lev, E., Liu, Y. P., Zeitler, P. K., Zhang, X., Zhang, J., Zurek, B., 2007. Geody-

512

namics of the southeastern Tibetan plateau from seismic anisotropy and geodesy.

513

Geology 35, 563–566.

514

Stadler, G., Gurnis, M., Burstedde, C., Wilcox, L. C., Alisic, L., Ghattas, O., 2010.

515

The dynamics of plate tectonics and mantle flow: From local to global scales.

516

Science 329, 1033–1038.

517

Steinberger, B., Calderwood, A., 2006. Models of large-scale viscous flow in the

518

earth’s mantle with constraints from mineral physics and surface observations.

519

Geophys. J. Int. 167, 1461–1481.

520

Steinberger, B., Schmeling, H., Marquart, G., 2001. Large-scale lithospheric stress

521

field and topography induced by global mantle circulation. Earth Planet. Sci.

522

Lett. 186, 75–91.

523

Torsvik, T. H., Tucker, R. D., Ashwal, L. D., Carter, L. M., Jamtveit, B., Vidyadha-

524

ran, K. T., Venkataramana, P., 2000. Late Cretaceous India-Madagascar fit and

525

timing of break-up related magmatism. Terra Nova 13, 220–224.

526

527

528

van der Voo, R., Spakman, W., Bijwaard, H., 1999. Tethyan subducted slabs under India. Earth Planet. Sci. Lett. 171, 7–20. van Hinsbergen, D. J., Steinberger, B., Doubrovine, P., Gassm¨oller, R., 2011. 22

529

Acceleration-deceleration cycles of India-Asia convergence: roles of mantle

530

plumes and continental collision.

531

2010JB008051.

532

533

534

535

J. Geophys. Res. in press.

doi:10.1029/

Windley, B. F., Allen, M. B., 1993. Mongolian plateau: Evidence for a late Cenozoic mantle plume under central Asia. Geology 21, 295–298. Yin, A., Harrison, T. M., 2000. Geologic evolution of the Himalayan-Tibetan orogen. Ann. Rev. Earth Planet. Sci. 28, 211–280.

536

Zhang, Z.-K., Wang, M., Gan, W., B¨urgmann, R., Wang, Q., Niu, Z., Sun, J., Wu, J.,

537

Hanrong, S., Xinzhao, Y., 2004. Continuous deformation of the Tibetan Plateau

538

from global positioning system data. Geology 32, 809–812.

539

Zhong, S., Zuber, M. T., Moresi, L., Gurnis, M., 2000. Role of temperature-

540

dependent viscosity and surface plates in spherical shell models of mantle con-

541

vection. J. Geophys. Res. 105, 11,063–11,082.

23

Fig. 1. Topography and deformation indicators for the Tethyan belt. A) Geodetic velocity field (Zhang et al., 2004; Gan et al., 2007; ArRajehi et al., 2010), topography (z), and seismicity (M ≥ 6, colored by depth, z, from the catalog of Engdahl et al., 1998). B) Residual topography relative to a regional mean, estimated by correcting for isostatic adjustment using the CRUST2.0 model (Bassin et al., 2000) (see supplementary information). Fig. 2. Reference flow model for the Tethyan belt. The density structure used for the computations is based on the SMEAN composite tomography (Becker and Boschi, 2002). A) Plate boundaries are treated as weak narrow belts, as depicted by the log10 of surface viscosity relative to the reference of 1021 Pas. B) Large-scale comparison between predicted (white) and observed (NUVEL-1A model, blue DeMets et al., 1994) velocities, in a Eurasia fixed reference frame. C) Predicted surface velocities (white vectors), geodetic velocities (orange, averaged from Fig. 1A), and dynamic topography of the Tethyan plate system. NUVEL-1A plate boundaries are shown in green. D) Horizontal (white vectors) and vertical flow (background shading) field at 300 km and 700 km depth. E and F) Cross-sections from Afar to Iran (E), and from the Carlsberg ridge to Central Asia (F) (trace of cross section marked in C), with non-dimensional temperature in the background. Fig. 3. Additional flow models for the Tethyan belt. Surface velocities and dynamic topography, as in Fig. 2C. A) Flow generated by density anomalies only, holding the surface motions fixed to zero outside India and Arabia, for the reference model of Fig. 2; B) Flow generated by density anomalies only (i.e. surface held fixed outside the India and Arabia) based on MIT08 P wave tomography (Li et al., 2008a); C) Flow generated by density anomalies holding the Australian plate moving at NUVEL-1A model rates (DeMets et al., 1994); D) Flow generated by density anomalies holding the African plate moving at NUVEL-1A rates; E) Flow generated by density anomalies restricted within the Wadati Benioff zones; F) Flow generated by density anomalies from reference models, including the first 100 km, but with upper mantle positive density anomalies restricted to Wadati Benioff zones and slow tomographic velocities only between 100 and 660 km.

24

Fig. 4. Three-dimensional view of the conveyor belt beneath India obtained by visualizing a cross-section of the reference model shown in Fig. 2.

25

70˚ 3 cm/yr

60˚ 50˚ 40˚

z [km] 4 2 0 -2 -4

30˚ 20˚ 10˚

d [km] 0˚

400

-10˚

200 0

A

20˚

40˚

60˚

80˚

100˚

120˚

140˚

-20˚

60˚

40˚

δztwb [km] 2 1 0 -1 -2

20˚

0˚

B

20˚

40˚

60˚

80˚

- Figure 1 -

100˚

120˚

140˚

-20˚

70˚

log10(η) 2 1 0 -1 -2

60˚ 50˚

60˚

d D

40˚ 30˚ 20˚

b

10˚

40˚

0˚

A

-10˚ 20˚

30˚

40˚

50˚

60˚

70˚

80˚

-20˚ 90˚ 100˚ 110˚ 120˚ 130˚ 140˚ 150˚

20˚

60˚

a

30˚ 0˚

δz dyn [km] 0˚ 2 1 0 -1 -2

c

B

120˚

60˚

C

180˚

20˚

40˚

300 km 3 cm/yr

60˚

vr [cm/yr] 2

60˚

D

80˚

100˚

700 km 1.5 cm/yr

60˚

vr [cm/yr] 2 40˚

0 -2

40˚

0 -2

20˚

40˚

60˚

80˚

100˚

Afar plume

a

b

120˚

140˚

20˚

20˚

0˚

0˚

-20˚ 20˚

40˚

60˚

Tethyan subduction

80˚

100˚

120˚

140˚

-20˚

Carlsberg ridge upwelling

d

c E

140˚

120˚

F temperature

Tethyan subduction

- Figure 2 -

60˚

60˚

δzdyn [km] 2 1 0 -1 -2 20˚

δzdyn [km] 2 1 0 -1 -2

40˚

40˚

20˚

0˚

0˚

A 20˚

30˚

40˚

50˚

60˚

70˚

80˚

-20˚ 90˚ 100˚ 110˚ 120˚ 130˚ 140˚ 150˚

B

20˚

30˚

40˚

50˚

60˚

70˚

80˚

-20˚ 90˚ 100˚ 110˚ 120˚ 130˚ 140˚ 150˚

60˚ δzdyn [km] 2 1 0 -1 -2

60˚

δzdyn [km] 2 1 0 -1 -2 20˚ 40˚

40˚

20˚

0˚

C

20˚

40˚

60˚

80˚

100˚

120˚

140˚

0˚

-20˚

D 20˚

60˚

3 cm/yr δzdyn [km] 2 1 0 -1 -2

40˚

20˚

40˚

60˚

80˚

100˚

120˚

140˚

60˚

3 cm/yr δzdyn [km] 2 1 0 -1 -2

40˚

20˚

0˚

E

20˚

40˚

60˚

80˚

100˚

120˚

140˚

-20˚

Figure 3

-20˚

0˚

F

20˚

40˚

60˚

80˚

100˚

120˚

140˚

-20˚

Figure 4

Supplementary material for: A mantle conveyor belt beneath the Tethyan collisional belt

Thorsten W. Becker ∗ Department of Earth Sciences, University of Southern California, Los Angeles, CA

Claudio Faccenna Dipartimento Scienze Geologiche, Universit`a Roma TRE and IGAG, CNR Rome

1

1 Methods

2

The surface mechanical boundary conditions for most of our models are shear stress

3

free (free slip), such that plate velocities are driven self-consistently by density

4

contrasts within the mantle. The lithosphere in our models is defined rheologically

5

as a relatively stiff layer from the surface to 100 km depth that is “broken” by weak

6

zones of reduced viscosity, as in previous modeling (e.g. King and Hager, 1990;

7

King et al., 1992; Zhong et al., 2000; Yoshida et al., 2001; Faccenna and Becker,

8

2010). The location of the weak zones (Fig. 2A) was assigned for the major plate ∗ Address: Department of Earth Sciences, University of Southern California, MC 0740, 3651 Trousdale Pkwy, Los Angeles, CA 90089–0740, USA. Phone: ++1 (213) 740 8365, Fax: ++1 (213) 740 8801 Email address: [email protected] (Thorsten W. Becker).

Preprint submitted to Elsevier

26 April 2011

9

boundaries following the NUVEL-1 geometry, and, for the case of Asia, located

10

at the boundary of the actively deforming region. All velocities are plotted in an

11

Eurasia fixed reference frame (from a best-fit of the velocities within that plate,

12

which typically show large intraplate deformation), and the core-mantle boundary

13

is set to free-slip.

14

We have verified that individual models have converged to a stable velocity solu-

15

tion in the presence of lateral viscosity variations. All results shown here have a

16

mesh resolution of ∼17 km and ∼26 km in the horizontal and vertical direction,

17

respectively. Resolution tests indicate that model velocities are stable to within a

18

few percent under successive mesh refinement. The numerical resolution employed

19

is therefore sufficient such that the model uncertainties are mainly due to imperfect

20

knowledge of the input models, such as tomography velocity structure, and not due

21

to computational limitations.

22

As described in Faccenna and Becker (2010), dynamic topography from flow was

23

computed by dividing the radial normal stresses at the surface, as predicted by

24

the flow model, by the product of gravitational acceleration (10 m/s2 ) and density

25

contrast between the mantle (density ρm = 3350 kg/m3 ) and air over land, and the

26

contrast between mantle and seawater (density ρw = 1020 kg/m3 ) over water, a

27

standard approximation. Residual topography as inferred from crustal models was

28

estimated by correcting the observed topography for isostatic adjustment given the

29

density structure of CRUST2.0 (Bassin et al., 2000). All topography estimates are

30

shown relative to a regional mean after smoothing by convolution with a Gaussian

31

of 250 km diameter. 2

32

2 Tests with different tomography models

33

We test the effect of using different tomographic models (Fig. S1). Our reference

34

mantle density model is based on the global, composite S wave SMEAN model

35

(Becker and Boschi, 2002), (Fig. S1). Alternative tomography models explored in-

36

clude the upper mantle SV model LH08 (Lebedev and van der Hilst, 2008), merged

37

with SMEAN in the lower mantle (Figs. S1B and S3), and the higher resolution,

38

global P model MITP08 (Li et al., 2008) (Figs. S1C and S2).

39

The effect of a set of lateral variations in viscosity on the flow predictions is shown

40

in Fig. S4, where continental keels down to 250 km (from the 3SMAC mode Nataf

41

and Ricard, 1996) were assigned to be more viscous than the mantle by a factor

42

of 100, and temperature-dependent viscosity was applied following the simplified

43

rheological description

44

η(T ) = η0 exp (E(T0 − T ))

(1)

45

where η0 is the mean layer viscosity, T0 a reference temperature (the mean layer

46

temperature), T the non-dimensional temperature as depicted in the cross-sectional

47

figures and inferred from tomography, and E the Frank-Kaminetskii parameter

48

(here E = 10).

49

References

50

Bassin, C., Laske, G., Masters, G., 2000. The current limits of resolution for surface

51

52

wave tomography in North America (abstract). Eos Trans. AGU 81, F897. Becker, T. W., Boschi, L., 2002.

A comparison of tomographic and geody3

53

namic mantle models.

54

2001GC000168.

55

56

Geochem., Geophys., Geosys. 3, 1.

doi:10.1029/

Faccenna, C., Becker, T. W., 2010. Shaping mobile belts by small-scale convection. Nature 465, 602–605.

57

King, S. D., Gable, C. W., Weinstein, S. A., 1992. Models of convection-driven

58

tectonic plates: a comparison of methods and results. Geophys. J. Int. 109, 481–

59

487.

60

King, S. D., Hager, B. H., 1990. The relationship between plate velocity and trench

61

viscosity in Newtonian and power-law subduction calculations. Geophys. Res.

62

Lett. 17, 2409–2412.

63

Lebedev, S., van der Hilst, R. D., 2008. Global upper-mantle tomography with the

64

automated multimode inversion of surface and S-wave forms. Geophys. J. Int.

65

173, 505–518.

66

Li, C., van der Hilst, R. D., Engdahl, E. R., Burdick, S., 2008. A new global model

67

for P wave speed variations in Earth’s mantle. Geochem., Geophys., Geosys. 9,

68

Q05018. doi:10.1029/2007GC001806.

69

Nataf, H.-C., Ricard, Y., 1996. 3SMAC: An a priori tomographic model of the

70

upper mantle based on geophysical modeling. Phys. Earth Planet. Inter. 95, 101–

71

122.

72

Yoshida, M., Honda, S., Kido, M., Iwase, Y., 2001. Numerical simulation for the

73

prediction of the plate motions: effects of lateral viscosity variations in the litho-

74

sphere. Earth Planets Space 53, 709–721.

75

Zhong, S., Zuber, M. T., Moresi, L., Gurnis, M., 2000. Role of temperature-

76

dependent viscosity and surface plates in spherical shell models of mantle con-

77

vection. J. Geophys. Res. 105, 11,063–11,082.

4

Fig. S1. Comparison of velocity anomaly maps of tomographic models used. A) Reference, global composite S wave model, SMEAN (Becker and Boschi, 2002); B) Upper mantle SV model LH08 (Lebedev and van der Hilst, 2008), merged with SMEAN in the lower mantle; C) Global, high resolution P wave model MITP08 (Li et al., 2008), all shown at 150, 250, 550 and 950 km depths.

Fig. S2. Additional flow solution based on different tomography, the temperature structure is now inferred from the MITP08 tomography model (Li et al., 2008) (Fig. S1C). A) Predicted surface velocities and dynamic topography (boundary in green) of the Tethyan plate system (white vectors), geodetic velocities (orange), all shown in a best-fit, Eurasia fixed reference frame. B) Global predicted (blue arrows) and observed (white arrows, NUVEL1A) large-scale velocities; C) Horizontal (white vectors) and vertical flow (background shading) field at 300 km and 730 km depth. Cross-sections from Afar to Iran (E) and from the Carlsberg ridge to Central Asia (F) (trace of cross section as in A), with non-dimensional temperature in the background.

Fig. S3. Additional flow solution based on the SV wave tomography model LH08 (Lebedev and van der Hilst, 2008) (Fig. S1B). All plots shown depict the same quantities as in Fig. S2, see there for explanation.

Fig. S4. Additional flow solution for reference model (as in Fig. 2) but including stiff keel beneath cratons and lateral viscosity variations as inferred from temperature (Frank Kaminetskii approximation). All plots shown depict the same quantities as in Fig. S2, see there for explanation.

5

20˚

a) SMEAN

60˚

40˚

60˚

80˚

100˚

120˚

140˚

z = 100 km

20˚

60˚

40˚

40˚

20˚

20˚

0˚

0˚

40˚

60˚

80˚

100˚

δvS [%] -2 20˚

40˚

0

60˚

δvS [%]

2

80˚

-2

100˚

120˚

140˚

z = 500 km

20˚

60˚

40˚

40˚

20˚

20˚

0˚

0˚

-20˚

40˚

60˚

0 80˚

2 100˚

40˚

0

60˚

δvS [%]

2

80˚

-2

100˚

120˚

140˚

z = 100 km

20˚

60˚

40˚

40˚

20˚

20˚

0˚

0˚

-20˚

40˚

60˚

0

80˚

2

100˚

40˚

-4 60˚

-2

0 80˚

2

4 100˚

δvS [%]

6

-2 120˚

140˚

z = 500 km

20˚

60˚

40˚

40˚

20˚

20˚

0˚

0˚

-20˚

40˚

60˚

0 80˚

2 100˚

140˚

-20˚

20˚

40˚

60˚

0

80˚

δvS [%]

2

-2

100˚

120˚

140˚

z = 100 km

20˚

60˚

40˚

40˚

20˚

20˚

0˚

0˚

40˚

60˚

0

2

80˚

100˚

0

1

80˚

100˚

0

1

120˚

140˚

z = 250 km

-20˚

-20˚ δvP [%] -1 20˚

60˚

120˚

z = 800 km

δvS [%] -2

60˚

140˚

-20˚

20˚

c) MITP08

120˚

z = 250 km

δvS [%] -6

60˚

140˚

-20˚

20˚

60˚

120˚

z = 800 km

δvS [%] -2

b) LH08_SMEAN

140˚

-20˚

-20˚

60˚

120˚

z = 250 km

40˚

60˚

0

1

80˚

100˚

δvP [%] -1

120˚

z = 500 km

140˚

20˚

60˚

40˚

40˚

20˚

20˚

0˚

0˚

-20˚

40˚

60˚

120˚

z = 800 km

-20˚ δvP [%] -1

0

1

δvP [%] -1

140˚

A)

60˚

d D δz dyn [km] 2 1 0 -1 -2

B)

b

40˚

20˚

60˚

a

30˚

0˚

0˚

c 60˚

-20˚

120˚

180˚ 20˚

40˚

60˚

C)

D)

z = 300 km 60˚

3 cm/yr

100˚

120˚

140˚

z = 730 km 60˚

1.5 cm/yr vr [cm/yr] 2

vr [cm/yr] 2 40˚

0

40˚

0 -2

-2 20˚

20˚

0˚

0˚

-20˚

E) a

80˚

20˚

40˚

60˚

80˚

100˚

120˚

140˚

-20˚ 20˚

40˚

60˚

80˚

100˚

120˚

140˚

b F) c temperature

d

A)

60˚

δz dyn [km] 2 1 0 -1 -2

d D b

40˚

20˚

B)

60˚

a

30˚ 0˚

0˚

c

-30˚

C)

120˚

60˚

-20˚

180˚ 20˚

40˚

z = 300 km

60˚

60˚

3 cm/yr

100˚

120˚

140˚

z = 730 km

D)

vr [cm/yr] 2

60˚

1.5 cm/yr vr [cm/yr] 2

40˚

0 -2

E) a

80˚

40˚

0 -2

20˚

20˚

0˚

0˚

-20˚ 20˚

40˚

60˚

80˚

100˚

120˚

140˚

b F) c

temperature

-20˚ 20˚

40˚

60˚

80˚

100˚

120˚

140˚

d