Nucleation and evolution of ridge-ridge-ridge triple junctions: Thermomechanical model and geometrical theory

Tectonophysics xxx (xxxx) xxx–xxx

Contents lists available at ScienceDirect

Tectonophysics journal homepage: www.elsevier.com/locate/tecto

Nucleation and evolution of ridge-ridge-ridge triple junctions: Thermomechanical model and geometrical theory Taras Geryaa,⁎, Evgueni Burovb,1 a b

ETH-Zurich, Institute of Geophysics, Department of Earth Sciences, Zurich, Switzerland CNRS, UMR 7193, Institut des Sciences de la Terre Paris (iSTeP), F-75005 Paris, France

A R T I C L E I N F O

A B S T R A C T

Keywords: Triple junction Quadruple junction Rifting Oceanic spreading Mid-ocean ridges Plate tectonics

Ridge-ridge-ridge triple junctions are among the most remarkable features of global plate tectonics but their nucleation and evolution remains incompletely understood. Here, we employ 3D numerical models to study the processes of the nucleation and evolution of triple junctions induced by multi-directional lithospheric extension. The simulations show that two major classes of junctions develop: (i) transient quadruple and triple plate rifting junctions formed by the initial plate breakup that are gradually converted into (ii) stable triple oceanic spreading junctions formed by the accretion of new lithosphere. Quadruple junctions break into two diverging triple oceanic spreading junctions connected by a linear spreading center lengthening with time. This process gradually decreases the length of deforming boundaries between four diverging rigid plates and thus the integral mechanical resistance of these boundaries to the spreading. The geometry of triple oceanic spreading junctions varies from asymmetrical T-junctions to ideal 120° junctions. The structure of these junctions includes two crucial tectonic elements: oceanic ridges (spreading centers) and intra-plate ranges (triple junction traces) dividing crust accreted from different spreading centers. The geometrical steady state is achieved within several million years. We propose a new geometrical theory of a migrating steady state triple junction, which describes its structure as a function of the relative plate velocities in a moving triple junction reference frame. The relative plate velocities define the orientations and growth rates of respective intra-plate ranges. Orientations and lengthening rates of oceanic ridges are in turn given by the average relative velocity of two adjacent plates. The steady state triple junction migration velocity satisfies the condition of minimal cumulative ridge lengthening rate weighted by the square root of the spreading velocity. Our analytical theory agrees well with the results of the numerical simulations and natural data.

1. Introduction Triple junctions are among the most remarkable features of global plate tectonics (McKenzie and Morgan, 1969; Patriat and Courtillot, 1984; Kleinrock and Morgan, 1988), but their nucleation and evolution, as well as their connection to magmatism and underlying mantle structure, remain a topic of debate (e.g., Nakanishi et al., 1999; Ligi et al., 1999; Bird et al., 1999; Larson et al., 2002; Pockalny et al., 2002; Georgen and Lin, 2002; Kotelkin et al., 2004; Verzhbitskii et al., 2006; Georgen, 2008; Mitchell et al., 2011 and references therein). Divergent ridge-ridge-ridge (R-R-R) triple junctions are distinct because they are stable at all ridge configurations and spreading rates (McKenzie and Morgan, 1969). Their geometry can vary from symmetrical, almost ideal 120° junctions to very asymmetrical T-junctions, in which one (slower spreading) ridge joins another (faster spreading) ridge at nearly

⁎

1

right angles without crossing it. Triple junction geometry is defined by both rigid plates' motions and (potentially asymmetric) plate accretion at the surrounding mid-ocean ridges (McKenzie and Morgan, 1969; Patriat and Courtillot, 1984). It is commonly accepted (although not quantitatively tested) that the geometry and stability of R-R-R triple junctions may be related to the intuitive geometric considerations that 3-branch configurations should be more “stable” compared to > 3branch configurations (e.g. quadruple R-R-R-R junctions) under conditions of long-term multi-directional extension on a 3D Earth surface (McKenzie and Morgan, 1969; Hey and Milholland, 1979). Indeed, it has been suggested based on limited natural observations that some quadruple junctions can be formed by continental plate breakup and exist for an extended period of time (Hey and Milholland, 1979), yet a consistent mechanical explanation or experimental/numerical demonstration of this process has not been provided. Observations also suggest

Corresponding author. E-mail address: [email protected] (T. Gerya). Deceased 9 October 2015.

http://dx.doi.org/10.1016/j.tecto.2017.10.020 Received 15 February 2017; Received in revised form 25 September 2017; Accepted 19 October 2017 0040-1951/ © 2017 Elsevier B.V. All rights reserved.

Please cite this article as: Gerya, T., Tectonophysics (2017), http://dx.doi.org/10.1016/j.tecto.2017.10.020

Tectonophysics xxx (xxxx) xxx–xxx

T. Gerya, E. Burov

that plume impingement on a stationary continental lithosphere results in three rifts oriented at angles of about 120° to each other, perhaps because this configuration requires the least work (Burke and Dewey, 1973). However, this apparently contradicts the results of recent numerical experiments suggesting that mantle plume impingement on a non-pre-stressed lithosphere should produce axisymmetric domalshaped features with multiple radiating and circular fractures (e.g., Gerya, 2014; Burov and Gerya, 2014). Plume impingement on a unidirectionally pre-stressed lithosphere results in the development of linear rift systems (Burov and Gerya, 2014; Koptev et al., 2015, 2016). Therefore, the initiation of R-R-R triple or R-R-R-R quadruple junction systems (Burke and Dewey, 1973; Hey and Milholland, 1979) should rather be related to more complex (e.g. bi-directionally or multi-directionally) pre-stressed lithosphere. Reconstructing the evolution of R-R-R triple junctions on the basis of observations is a complex task (e.g., Apotria and Gray, 1985, 1988; Kleinrock and Morgan, 1988; Patriat and Parson, 1989; Mitchell and Parson, 1993; Dyment, 1993; Honsho et al., 1996; Ligi et al., 1997; Nakanishi et al., 1999; Ligi et al., 1999; Bird et al., 1999; Larson et al., 2002; Pockalny et al., 2002) since currently existing junctions represent a very fragmentary record and early evolutionary stages are poorly known. From this viewpoint, the Afar triple junction between the Red Sea rift, the Gulf of Aden rift and the Main Ethiopian rift, represents the only incipient R-R-R triple junction example in the globe. The associated complex motions of Arabian, Nubian, and Somalian plates (McKenzie et al., 1970; Le Pichon and Gaulier, 1988; Fournier et al., 2001; McClusky et al., 2003) imply a multi-directional extension of the lithosphere in this region, which may thus play a crucial role in the nucleation of the Afar triple junction. The limited availability of direct observations on the nucleation and evolution of triple junctions thus calls for a systematic 3D numerical modeling effort of associated lithospheric and mantle processes (e.g., Georgen and Lin, 2002; Georgen, 2008; Georgen and Sankar, 2010; Beutel et al., 2010; Mitchell et al., 2011; Brune and Autin, 2013; Brune et al., 2014; Gerya, 2014; Burov and Gerya, 2014; Koptev et al., 2015, 2016; Dordevic and Georgen, 2016). In particular, the spontaneous appearance of ridge-ridge-ridge and ridge-ridge-transform triple junctions during multi-directional extension processes has been recently documented in numerical experiments on plume-induced subduction initiation (Gerya et al., 2015). Taking into account the complex tectonomagmatic nature of rifting and spreading processes, a systematic highresolution 3D magmatic-thermomechanical approach (e.g., Gerya, 2013, 2014; Gerya et al., 2015) is needed to better understand the physics of triple junctions. Here, we numerically study how triple junctions form under bi-directional and multi-directional lithospheric extension. We use high-resolution 3D numerical magmatic-thermomechanical experiments that take into account a realistic thermorheological structure of the lithosphere and account for crustal and lithospheric accretion processes. Based on the results of our experiments we discriminate between (i) plate rifting junctions formed by the initial plate fragmentation and subsequently emerging (ii) oceanic spreading junctions controlled by the new oceanic crust and lithosphere accretion. We document major tectonic elements and the evolution of the junctions and develop a new geometrical theory of a migrating steady state triple oceanic spreading junction, which we then generalize for multibranch junctions.

a)

b)

Fig. 1. Two numerical model setups explored in this study: a) bi-directional extension of a single plate with a thermal perturbation and/or weak inclusion; b) bi-directional extension imposed on one side of pre-existing mid-ocean ridge. Top views (98 × 98 km) are shown for 98 × 98 × 50 km models. Arrows show velocity boundary conditions. Star with SWNE symbols show polar coordinates used for visualization of velocity direction.

implemented through a 5 km thick “sticky” water layer (e.g., Gerya, 2013; Crameri et al., 2012; Schmeling et al., 2008), allows for large strains and spontaneous crust growth by magmatic accretion. The employed numerical technique (Gerya, 2010a, 2010b, 2013) is based on a combination of a finite difference method applied on a uniformly spaced staggered grid with the marker-in-cell technique. 2.1. Model setup Two initial model setups were explored that correspond to two different geodynamic settings leading to the generation of oceanic R-RR triple junctions in nature. The first setup (Fig. 1a) corresponds to a bidirectional breakup of an oceanic plate with some initial perturbation (thermal anomaly and/or weak inclusion). This setup is (to some degree) inspired by previous hypotheses for the formation of continental and oceanic triple junctions due to an impingement of a plume on a (stressed) plate (e.g., Burke and Dewey, 1973; Kotelkin et al., 2004; Verzhbitskii et al., 2006), although the size of our regional numerical model is too small to investigate the complete process of plume-lithosphere interaction (Burov and Gerya, 2014; Koptev et al., 2015, 2016).

2. Numerical model This work employs high-resolution 3D Cartesian magmatic-thermomechanical numerical models of oceanic rifting and spreading based on I3ELVIS code (Gerya and Yuen, 2007; Gerya, 2013). Although a spherical geometry for a global plate mosaic is important for creating multi-directional extension, it is not crucial for our regional-scale models governed by prescribed plate kinematics. The Eulerian-Lagrangian visco-plastic model, which contains an internal free surface 2

Tectonophysics xxx (xxxx) xxx–xxx

T. Gerya, E. Burov

where D/Dt is the substantive time derivative, CP is isobaric heat capacity (J/K/kg), T is temperature (K), qi are heat flux components (W/ m2) and Hr, Hs, HL and Ha are respectively radioactive, shear, latent and adiabatic heating terms (W/m3).

The Afar triple junction formed by multidirectional extension of continental lithosphere is the obvious modern prototype for this model setup. However, we initially define thin oceanic rather than thick continental lithosphere in our relatively small-scale highly simplified plate breakup models. The initial thermal structure is defined either by a half-space cooling age or by a linear geotherm (273 °C at the surface and 1315 °C at 30 km depth) in a 30 km thick oceanic lithosphere with 8 km thick mafic crust. Such relatively thin and young lithosphere is characteristic for oceanic triple junctions and allows the inclusion of a 5 km thick top water layer and 15 km thick hot asthenospheric mantle into our relatively shallow (50 km) high-resolution models. The asthenospheric mantle is adiabatic with a 0.5 °C/km vertical temperature gradient. The second setup (Fig. 1b) corresponds to bi-directional extension which is imposed on one of the lithospheric plates formed at a pre-existing linear oceanic ridge. This setup investigates the generation of asymmetrical oceanic triple junctions (or T-junctions), such as the Galapagos triple junction, the Cocos-Nazca Rift or the Rodrigues Triple Junction. Those are thought to form as the result of multi-directional stresses imposed on some pre-existing oceanic plates by changes in global plate motions. The thermal structure of this setup corresponds to unidirectional oceanic spreading with lithosphere cooling ages increasing outward from a pre-existing linear ridge. The 98 × 98 × 50 km Eulerian computational domain used in both types of models (Fig. 1) is resolved with a regular rectangular grid of 197 × 197 × 101 nodes and contains 34 million randomly distributed Lagrangian markers. Larger 202 × 202 × 50 km models with 405 × 405 × 101 nodes and 130 million markers have also been tested (Table 1), which demonstrated that the main conclusions of this paper based on smaller scale models are robust. The momentum, mass and heat conservation equations are solved with the thermomechanical code I3ELVIS (Gerya and Yuen, 2007) on a non-deforming Eulerian grid whereas the advection of transport properties including viscosity, plastic strain, temperature etc. is performed by advecting the Lagrangian markers. Extensional velocity boundary conditions are prescribed along the entire vertical boundaries of the model (Fig. 1). Compensating vertical influx velocities through the upper and lower boundaries are chosen to ensure conservation of volume of the model domain and a constant average 5 km thickness of the seawater layer (Gerya, 2013). Lagrangian markers leave the Eulerian model domain through the lateral boundaries. Material enters through the top and the bottom of the model as water and mantle markers, respectively. This Eulerian-Lagrangian numerical modeling scheme with open boundaries allows for an infinitely long plate separation with the use of a laterally limited Eulerian computational domain (e.g., Gerya, 2013).

2.3. Modeling of rifting and spreading processes We adopted the tectono-magmatic numerical model of rifting and spreading developed by Gerya (2013), which accounts for the following four key physical processes: (i) thermal accretion of the oceanic mantle lithosphere resulting in growth of the plate thickness, (ii) partial melting of the asthenospheric mantle, melt extraction and percolation toward the ridge resulting in crustal growth, (iii) magmatic accretion of the new crust beneath the ridge and (iv) hydrothermal circulation at the ridge axis, resulting in excess cooling of the crust. These physical processes are included in our numerical model in a simplified parameterized manner. Thermal accretion of the mantle lithosphere is modeled by solving the heat conduction equation combined with a temperature-dependent viscosity for the non-molten mantle (dry olivine flow law, Ranalli, 1995). This accounts for conductive cooling of asthenospheric mantle, which becomes rheologically strong and accretes spontaneously to the bottom of the lithosphere. Hydrothermal circulation at the ridge axis induces rapid cooling of the new crust (e.g., Theissen-Krah et al., 2011), which is parameterized with an enhanced thermal conductivity of the crust according to (Gregg et al., 2009)

k eff = k + k 0 (Nu − 1) × exp{A[2 − (T − 273)/(Tmax − 273) − y / ymax ]}, (4) where k is the thermal conductivity of dry rocks (W/m/K) [k = 0.73 + 1293 / (T + 77) for the mantle and k = 1.18 + 474 / (T + 77) for the crust (Clauser and Huenges, 1995)], k0 = 3 W/m/K is the reference thermal conductivity, Nu = 2, the Nusselt number for the hydrothermal circulation (Gregg et al., 2009; Gerya, 2013), A = 0.75 is a smoothing factor, Tmax is the cutoff maximum temperature (873 K), y is depth (km), and ymax = 6 km is the cutoff maximum depth of hydrothermal circulation. Also, in order to ensure an efficient heat transfer from the upper surface of the plate in the Eulerian model, the thermal conductivity of the water layer above this plate is taken to be a hundred times higher (200 W/m/K) than that of the dry mantle (1–4 W/m/K). Partial melting of the asthenospheric mantle, melt extraction and percolation toward the ridge. According to our model, crustal growth at the ridge is balanced by the melt production and extraction from the mantle. However, melt percolation toward the ridge (e.g., Katz, 2010) is not modeled directly and assumed to be nearly instantaneous (Connolly et al., 2009). The standard (i.e. without melt extraction) volumetric degree of mantle melting M0 changes with pressure and temperature according to the parameterized batch melting model of Katz et al. (2003). Lagrangian markers track the amount of melt extracted during the evolution of each experiment. The total amount of melt, M, for every marker takes into account the amount of previously extracted melt and is calculated as (Nikolaeva et al., 2008).

2.2. Governing equations We modeled 3D creeping flow using the extended Boussinesq approximation. The incompressible continuity equation takes the form

div (v ) = 0,

(1)

where v is the velocity vector (m/s). The Stokes equation of slow flow in Einstein notation takes the form

∂σij′ ∂x j

−

∂P = −ρgi , ∂x i

M = M0 − Σn Mext , (2)

where Σn Mext is the total melt fraction extracted during the previous n extraction episodes. The rock is considered non-molten (refractory) when the extracted melt fraction is larger than the standard one (i.e. when Σn Mext > M0). If M > 0 for a given marker, the melt fraction Mext = M is extracted and Σn Mext is updated. Similarly to previous rifting/spreading models (Gregg et al., 2009; Gerya, 2013), melt is moved upward and then along the mantle solidus surface topography until it reaches some local maximum in this topography where it can accumulate. In order to ensure conservation of melt volume and account for mantle compaction and subsidence in response to melt

where xi, xj are coordinates (m), σij′ are components of the deviatoric stress tensor (Pa), P is pressure (Pa), ρ is density (kg/m3) and gi are gravitational acceleration components (m/s2). The energy conservation equation is solved in a Lagrangian reference frame using temperature advection with markers (Gerya and Yuen, 2007)

ρCP

∂q DT = − i + Hr + Hs + HL + Ha, ∂x i Dt

(5)

(3) 3

1

1

trac

trad

4

98

98

202

1

1

2c

2c

1

tram Fig.5

tranb Fig.4d

trao

trap

traq

98

98

202

1

98

98

tral

1

1

Fig.4a

98

98

98

98

98

98

98

Size x=y (km)

trak

traj

1

Fig.3

trai

1

1

Fig.2

trah

trag

1

traf

Fig.6

1

trae

Fig.4c

Setup

Model figure

1.89/1.89/1.89/ 1.89

1.89/1.89/1.89/ 1.89

1.89/1.89/1.89/ 1.89

1.89/1.89/1.89/ 1.89

1.89/1.89/1.89/ 1.89

1.89/1.89/1.89/ 1.89

1.89/1.89/1.89/ 1.89

1.89/1.89/1.89/ 1.89

1.89/1.89/1.89/ 1.89 1.89/1.89/1.89/ 1.89

1.89/1.89/1.89/ 1.89

2.52/2.52/1.26/ 1.26

1.89/1.89/1.89/ 1.89

1.26/1.26/0.63/ 0.63 0.95/0.95/0.95/ 0.95

Outward velocity at W/E/S/N boundaries (cm/yr)

Table 1 Conditions and results of numerical experiments.

Cooling age 3 → 0. 1 → 3 Ma Cooling age 5 → 0. 1 → 5 Ma Cooling age 3 Ma

Cooling age 3 → 0.01 → 3 Ma

Cooling age 3 → 0.01 → 3 Ma

Cooling age 5 Ma

Cooling age 3 Ma

Cooling age 3 Ma

Cooling age 2 Ma Cooling age 3 Ma

Linear 30 km

Linear 30 km

Linear 30 km

Cooling age 2 Ma Cooling age 2 Ma

Plate geotherm

WI, two 20 × 20 km, at x = 30 km, y = 70 km and x = 70 km, y = 30 km, cooling age 3 → 1 Ma

2 TRJ connected by a NW-SE rift, NW-SE spreading initiation

TJ with an orthogonal EW rift, EW rift → EW ridge, TJ → TRJ 3 NW-SE rifts forming 3 TJs along the NS ridge

Central NS ridge, cooling age 0.5 → 0.1 → 0.5 Ma Central NS ridge, cooling age 0.5 → 0.1 → 0.5 Ma

Symmetrical QJ surrounded by ridges

Symmetrical QJ surrounded by ridges

Central QJ with NS and EW ridges

Central QJ with NS and EW ridges

Linear N-S ridge with orthogonal rifts forming TJ

Linear N-S ridge, QJ formed by the N-S ridge and W-E rift

LINEAR N-S rift with elliptic SC

Symmetrical QJ with circular SC surrounded by rifts

Linear NE-SW rift

Asymmetrical QJ surrounded by rifts

Linear N-S rift

Symmetrical QJ with circular SC surrounded by rifts

Linear N-S rift with elliptic SC Asymmetrical QJ surrounded by half-grabens

Initial rifting patterna (0–2 Myr)

WI, central 40 × 202 km, cooling age 5 → 0.1 Ma

WI, central 20 × 98 km, cooling age 3 → 0.1 Ma

WI, central 20 × 98 km, cooling age 3 → 2 Ma

WI, central 20 × 20 km, cooling age 3 → 2 Ma

WI

WI

centraL 10 × 10 km at 25–30 km depth, T = 1315 °C

Central 10 × 10 km at 25–30 km depth, T = 1315 °C

Central 10 × 10 km, cooling age 1 Ma

Central 10 × 10 km, cooling age 1 Ma

Perturbationa

Central NW-SE rift → EW ridge, deactivation of other rifts, central TJ → TRJ Two migrating TRJ separated by lengthening NW-SE ridge

E-ward migrating nearly ideal 120° TRJ

QJ breakup, two migrating TRJ separated by lengthening NW-SE ridge

QJ breakup, two migrating TRJ separated by lengthening NW-SE ridge

Two migrating TRJ separated by lengthening NNE-SSW ridge

Two migrating TRJ separated by lengthening NW-SE ridge

Single N-S ridge with orthogonal rifts forming TJ, TJ transform to TRJ

Two migrating TRJ separated by lengthening NNW-SSE ridge

Single NE-SW ridge

Two migrating TRJ separated by lengthening NE-SW ridge

Single N-S ridge with orthogonal rifts forming TJ

Two migrating TRJ separated by lengthening NE-SW ridge

Single N-S ridge with orthogonal rifts forming TJ Two migrating TRJ separated by lengthening NW-SE ridge

Young spreading patterna (2–4 Myr)

0.51/0.00/7.10

− 0.25/0.25/ 7.27 0.25/− 0.25/ 7.27 0.25/0.25/7.27 −0.25/ − 0.25/7.27 0.25/0.25/7.27 −0.25/ − 0.25/7.27 0.25/0.25/7.27 −0.25/ − 0.25/7.27 − 0.25/0.25/ 7.27 0.25/− 0.25/ 7.27 − 0.25/0.25/ 7.27 0.25/− 0.25/ 7.27 0.51/0.00/7.10

− 0.12/0.12/ 2.59 0.12/− 0.12/ 2.59 0.25/0.25/7.27 −0.25/ − 0.25/7.27 0.16/0.27/6.77 −0.16/ − 0.27/6.77 0.25/0.25/7.27 0.00/0.00/4.35 −0.07/ − 0.81/4.22 No TJ/TRJ

No TJ/TRJ

vxm/vym/Rmf (cm/yr)/(cm/ yr)/(cm3/2/ yr3/2)

− 0.25/0.25/ 7.27 0.25/− 0.25/ 7.27 (continued on next page) Two migrating TRJ separated by lengthening NW-SE ridge

E-ward migrating TRJ

E-ward migrating nearly ideal 120° TRJ

Two migrating TRJ separated by lengthening NW-SE ridge

Two migrating TRJ separated by lengthening NW-SE ridge

Two migrating TRJ separated by lengthening curved NE-SW ridge

Two migrating TRJ separated by lengthening NE-SW ridge

Two migrating TRJ separated by lengthening NE-SW ridge

Two migrating TRJ separated by lengthening NW-SE ridge

Single NE-SW ridge

Two migrating TRJ separated by lengthening NE-SW ridge, after 10 Ma - single migrating TRJ

Single NNE-SSW ridge with two migrating and one stationary TJ

Single migrating TRJ

Single migrating TRJ with coarse topography patterns

Single NNE-SSW ridge

Mature spreading patterna (> 4 Myr)

T. Gerya, E. Burov

Tectonophysics xxx (xxxx) xxx–xxx

1

trar

98

98

98

98

98

2d,e

2d,e

2d

2d

2d

trawa

trawb Fig.10

trawc Fig.7

trawd Fig.13b

trawe

5

98

1.89/0.95/1.64/ 1.64

1.89/1.89/1.89/ 1.89

3.16/3.16/0.95/ 0.95

3.16/3.16/0.95/ 0.95

3.16/3.16/0.95/ 0.95

3.16/3.16/0.95/ 0.95

1.89/1.89/0.95/ 0.95

1.89/1.26–2.52/ 1.89/1.89

1.89/1.26–2.52/ 1.89/1.89

1.89/1.89/1.89/ 1.89

1.89/1.89/1.89/ 1.89

Outward velocity at W/E/S/N boundaries (cm/yr)

Cooling age 2 → 0. 1 → 4 Ma Cooling age 2 → 0. 1 → 4 Ma Cooling age 2 → 0. 1 → 4 Ma Cooling age 2 → 0. 1 → 4 Ma Cooling age 2 → 0. 1 → 4 Ma Cooling age 2 → 0. 1 → 4 Ma Cooling age 2 → 0. 1 → 4 Ma Cooling age 2 → 0. 1 → 4 Ma Cooling age 2 → 0. 1 → 4 Ma

Cooling age 3 Ma

Cooling age 3 Ma

Plate geotherm

NS ridge at x = 25 km, cooling age 0.5 → 0.1 → 0.5 Ma

NS ridge at x = 25 km, cooling age 0.5 → 0.1 → 0.5 Ma

TJ, NS ridge at x = 25 km, cooling age 0.1 Ma, EW ridge at x = 25–98 km, y = 49 km, cooling age 0.1 → 3 Ma TJ, NS ridge at x = 25 km, cooling age 0.1 Ma, EW ridge at x = 25–98 km, y = 49 km, cooling age 0.1 → 3 Ma NS ridge at x = 25 km, cooling age 0.5 → 0.1 → 0.5 Ma

NS ridge at x = 25 km, cooling age 0.5 → 0.1 → 0.5 Ma

TJ, NS ridge at x = 25 km, cooling age 0.1 Ma, EW ridge at x = 25–98 km, y = 49 km, cooling age 0.1 Ma NS ridge at x = 25 km, cooling age 0.5 → 0.1 → 0.5 Ma

NS ridge at x = 25 km, cooling age 0.5 → 0.1 → 0.5 Ma

WI, one 20 × 20 km at x = 70 km, y = 40 km, cooling age 3 → 1 Ma

WI, two 20 × 20 km, at x = 30 km, y = 60 km and x = 70 km, y = 40 km, cooling age 3 → 1 Ma

Perturbationa

TJ with an orthogonal EW rift, EW rift → EW ridge, TJ → TRJ TJ with an orthogonal EW rift, EW rift → EW ridge, TJ→ TRJ

TJ with an orthogonal EW ridge

TJ with an orthogonal EW rift, EW rift → EW ridge

TJ with an orthogonal EW rift, EW rift → 2 EW ridges with 2 TJs TJ with an orthogonal EW rift, EW rift → EW ridge

TJ with an orthogonal EW rift, EW rift → EW ridge

NW-SE and SW-NE rifts forming 2 TJs along the NS ridge TJ → TRJ, NE-ward migration of TRJ

Stationary nearly ideal 120° TRJ

Bending of EW ridge, slow E-ward migration of TJ, development of triangular slow spreading region Bending of NS ridge, slow E-ward migration of TJs, deactivation of one TJ Bending of NS ridge, slow E-ward migration of TJ, development of triangular slow spreading region Bending of NS ridge, slow E-ward migration of TJ, development of triangular slow spreading region Bending of NS ridge, slow E-ward migration of TJ, development of triangular slow spreading region E-ward migrating nearly ideal 120° TRJ

NE-ward migration of TRJ

Deactivation of SW-NE rift, NW-SE rift → EW ridge, TJ → TRJ

Two migrating TRJ separated by lengthening NW-SE ridge

c

b

Stationary nearly ideal 120° TRJ

Bending of EW ridge, nearly stationary TJ and triangular slow spreading region Bending of EW ridge, nearly stationary TJ and narrow triangular slow spreading region Bending of NS ridge, nearly stationary TJ and triangular slow spreading region Bending of NS ridge, nearly stationary TJ and triangular slow spreading region Bending of NS ridge, nearly stationary TJ and triangular slow spreading region E-ward migrating nearly ideal 120° TRJ

NE-ward migration of TRJ

No result (model stopped at 3.5 Ma)

No result (model stopped at 2.6 Ma)

Two migrating TRJ separated by lengthening NW-SE ridge

Two migrating TRJ separated by lengthening NW-SE ridge

2 TRJ connected by a NWWSEE rift, NW-SE spreading initiation Asymmetrical QJ surrounded by rifts, QJ breaks into 2 TRJ

Mature spreading patterna (> 4 Myr)

Young spreading patterna (2–4 Myr)

Initial rifting patterna (0–2 Myr)

Abbreviations: SC = spreading center, TRJ = triple junction, TJ = T-junction, QJ = quadruple junction, WI = central weak 2 × 2 × 2 km inclusion (γ =10) at 13 km depth. No strain weakening of the crust and mantle. Free slip condition at 50% of N,S boundaries (Fig. 1b). d Free slip condition at 25% of N,S boundaries (Fig. 1b). e Step-wise increase (doubling) of the outward velocity in the middle of the E-boundary (Fig. 1b). f Values computed with Eq. (24) for observed TJ and TRJ. g Nu = 3.

a

Fig.9

2d

2d

trax

tray

98

2d

trawgg

Fig.8

98

2d

trawfg

98

98

1

98

Size x=y (km)

tras

Fig.4b

Setup

Model figure

Table 1 (continued)

0.00/0.00/5.14

0.51/0.00/7.10

0.13/0.00/6.67

0.13/0.00/6.67

0.13/0.00/6.67

0.13/0.00/6.67

0.18/0.00/4.36

0.58/0.31/7.15

− 0.25/0.25/ 7.27 0.25/− 0.25/ 7.27 − 0.25/0.25/ 7.27 0.25/− 0.25/ 7.27 0.58/0.31/7.15

vxm/vym/Rmf (cm/yr)/(cm/ yr)/(cm3/2/ yr3/2)

T. Gerya, E. Burov

Tectonophysics xxx (xxxx) xxx–xxx

Tectonophysics xxx (xxxx) xxx–xxx

T. Gerya, E. Burov

extension),

extraction, melt addition to the bottom of the magma region is performed at every time step by converting the shallowest markers of hot, partially molten mantle into magma markers. The total volume of these magma markers matches the total volume of extracted melt computed for the time step. Magmatic accretion of the new crust is modeled by spontaneous cooling and crystallization of melts at the walls of the lower-crustal magma regions (e.g., Wanless and Shaw, 2012). This simple crust accretion algorithm does not account for volcanic and plutonic (e.g. dyking, Buck et al., 2005; Wanless and Shaw, 2012) processes above the magma regions, neither does it account for internal convection, melt segregation and crystal differentiation inside these regions (Wanless and Shaw, 2012). The effective density of the molten crust in the magma region is calculated as

ρ ρeff = ρsolid ⎛⎜1 − M + M 0molten ⎟⎞, ρ0solid ⎠ ⎝

ϕγ = 0.6(1 − γ/γ0) for γ ≤ γ0 and ϕγ = 0 for γ > γ0 when P ≥ Pf (confined fracture)

M = (T − Tsolidus)/(Tliquidus − Tsolidus) when Tsolidus < T < Tliquidus,

(7b)

M = 1 when T > Tliquidus,

(7c)

α eff

Q ∂M ⎞ =α+ρ L⎛ , T ⎝ ∂P ⎠T = const

∫

1 (ε̇ij (plastic ) )2 dt − 2

∫ ε̇healing dt,

1 (σij′ )2 , 2

(15)

(16)

We performed 26 numerical experiments (Table 1) by varying different model parameters for both model setups. The results of our numerical experiments systematically predict the development of quadruple, triple and T-junctions as a result of plate rifting and subsequent oceanic spreading. In the following, we will describe the numerical results and the influences of various model parameters on the formation and evolution of different types of junctions.

(9)

3.1. Bi-directional extension of a single plate 3.1.1. Reference model development Fig. 2 shows the development of the reference model (trag, Table 1) with a central weak inclusion prescribed as a 2 × 2 × 2 km highly strained region (γ = 10) in the lithospheric mantle at 13 km depth. Extension rates prescribed along SN and EW directions are identical for this model. The geodynamic evolution of the model can be subdivided into four stages with gradual transitions between them: (1) radial extension of the plate (0–1 Ma, Fig. 2a), (2) quadruple junction formation by plate rifting (1–1.5 Ma, Fig. 2b), (3) breakup (splitting) of the quadruple junction into two triple junctions by incipient oceanic spreading (1.5–2.5 Ma, Fig. 2c), (4) outward migration of two stable triple junctions (> 2.5 Ma, Fig. 2d). During the first stage, symmetrical radial extension of the plate dominates in accordance with the prescribed homogeneous velocity boundary conditions mimicking lithospheric deformation above a spreading plume (Burov and Gerya, 2014; Gerya, 2014). Tensile deformation is focused in the center of the model above the initial weak inclusion, which produces a ring-like depression (Fig. 2a). During the second stage, the symmetry of the radial deformation breaks by rapid plate rifting, which produces an irregular cross-like quadruple junction (Fig. 2b). The third stage is manifested by the onset of oceanic spreading, which rapidly destabilizes the initially formed quadruple junction by breaking it into two diverging triple junctions. A small T-junction forms within one of the moving plates by

(10)

where Cp = 1000 J/kg is the heat capacity of the solid crust and QL = 380 kJ/kg is the latent heat of crystallization of the crust (Turcotte and Schubert, 2002). 2.4. Rheological model The rheological model assumes a constant low viscosity (1018 Pa s) for the partially molten crustal and mantle rocks, and a visco-plastic rheology for the lithospheric plates, with a realistic temperature and strain rate dependent viscosity computed according to experimentally determined flow laws (Ranalli, 1995). We use the dry olivine flow law for the mantle and the plagioclase (An75) flow law for the crust. Cutoff limits of 1018 Pa s and 1024 Pa s were applied to limit minimal and maximal rock viscosity, respectively. Standard brittle/plastic rheology of solid rocks assumes fracture-related strain weakening (e.g., Lavier et al., 2000; Huismans and Beaumont, 2002; Choi et al., 2008; Allken et al., 2011, 2012; Gerya, 2010b, 2013), and is implemented by using a plate strength limitation of the form (Gerya, 2013)

σII ≤ Cγ + ϕγ (P − Pf )

(14)

3. Results of numerical experiments

(8)

where α = 3 × 10− 5 1/K and β = 10− 5 1/MPa are the thermal expansion and compressibility of the crust. The effect of latent heating due to equilibrium crystallization of the crust from the magma regions is included implicitly by increasing the effective heat capacity (Cpeff) and the thermal expansion (αeff) of the partially crystallized/molten rocks (0 < M < 1), calculated as (Gerya, 2010a):

∂M ⎞ Cpeff = Cp + QL ⎛ , ⎝ ∂T ⎠P = const

Pf = yρf gy ,

where σII is the square root of the second stress invariant (Pa), Pf is fluid pressure (Pa), y is the vertical coordinate (m), gy = 9.81 m/s2 is the vertical gravitational acceleration component, ρf = 1000 kg/m3 is water density, γ ≥ 0 is integrated plastic strain (γ0 = 1 is the upper strain limit for the fracture-related weakening), t is time (s), εij̇ (plastic ) is ̇ the plastic strain rate tensor (1/s), εhealing is the fracture healing rate (1/ s), ϕγ is the internal friction coefficient which depends on the plastic strain γ, Cγ is the rock strength at P − Pf = 0 (for both confined and tensile fracture), which also depends on the plastic strain γ, and C0 = 10 MPa and C1 = 3 MPa are the initial and final strength values for the fracture-related weakening, respectively. Partial healing of deactivated, slowly creeping fractures was applied in the experiments by ̇ using an imposed constant healing rate εhealing = 10− 13 1/s (Gerya, 2013) which reduces the accumulated plastic strain γ with time and allows for the reproduction of changes of mid-ocean ridge morphology with changing spreading rate (Püthe and Gerya, 2013).

where Tsolidus = 1327 + 0.091P and Tliquidus = 1423 + 0.105P are, respectively, solidus and liquidus temperature (K) of the crust (Hess, 1989) at a given pressure P (MPa). ρ0solid = 3000 kg/m3 and ρ0molten = 2800 kg/m3 are the standard densities of solid and molten crust, respectively and ρ0solid is the density of solid crust at given P and T computed from

ρsolid = ρ0solid × [1 − α(T − 298)] × [1 + β(P − 0.1)],

(13)

σII =

where local melt fraction M is computed from the simplified linear batch melting model (Gerya, 2010a) (7a)

Cγ = C0 + (C1 − C0 )γ/γ0 for γ ≤ γ0 and Cγ = C1 for γ > γ0

γ=

(6)

M = 0 when T < Tsolidus,

(12)

(11)

ϕγ = 1 when P < Pf (tensile-like fracture oriented perpendicular to 6

Tectonophysics xxx (xxxx) xxx–xxx

T. Gerya, E. Burov

Fig. 2. Evolution of the numerical model with a central 2 × 2 × 2 km weak inclusion (model trag in Table 1). Top views (98 × 98 km) are shown for 98 × 98 × 50 km models. Left column shows topography evolution (colors correspond to elevation relative to the initial crustal surface located at 5 km level). Right column shows evolution of the velocity field (arrows, colors correspond to velocity direction in polar coordinates, Fig. 1) at the level of 2.75 km below the initial crustal surface.

a)

b)

c)

d)

7

Tectonophysics xxx (xxxx) xxx–xxx

T. Gerya, E. Burov

A strain weakening rheology strongly influences faulting and topography patterns during both rifting and spreading (e.g., Gerya, 2010b, 2013) but does not seem to have a very strong effect on the stability and geometry of triple junctions. Models without strain weakening show smooth topographies without off-ridge normal faulting, but the overall evolution and geometry of diverging triple junctions is very similar to the model with strain weakening (cf. Figs. 5 and 4d). The symmetry of bi-directional extension velocities strongly affects all stages of the model development (cf. Figs. 6, 2 and 4c). Models with non-uniform bi-directional extension develop a single rift, which is oriented perpendicular to the direction of maximal extension velocities (Fig. 6a). At a later stage, a single gently curved ridge develops from this rift, which becomes truncated by several orthogonal secondary rifts into a series of spreading sections separated by T-junctions (Fig. 6b). The orthogonal rifts evolve into narrow triangular slow spreading centers surrounded by intra-plate ranges. Both migrating and stationary T-junctions related to these centers develop. The initial ridge becomes curved and its different sections spread at different rates and therefore have different topographic and faulting patterns (Fig. 6c,d). This result again supports our choice of the second model setup (Fig. 1b) simulating T-junction formation at a pre-existing straight ridge. It should also be noted that intra-plate ranges are less pronounced for cases with non-uniform bi-directional spreading compared to symmetrical extension models (cf. Figs. 3d and 6d). Decreased magnitude of the symmetrical bi-directional extension velocity produces coarser and more asymmetrical rifting and spreading patterns, which agrees well with results of previous uni-directional spreading modeling (Püthe and Gerya, 2013). Initially, pronounced half grabens form around an irregular quadruple junction (cf. Fig. 4c, 1.30 Ma and Fig. 3b). At a later stage, coarse spreading patterns develop, which include a (somewhat irregular) triple junction surrounded by poorly pronounced spreading centers and intra-plate ranges (cf. Fig. 4c, 5.76 Ma and Fig. 3d).

intra-plate rifting oriented nearly orthogonal to the incipient oceanic ridge bounding this plate (Fig. 2c). This result motivated our choice of the second model setup aimed at investigating how T-junctions form by plate motions rearrangement. During the subsequent oceanic spreading stage, continued plate accretion causes outward migration of two stable triple junctions. The migration is likely caused by the non-zero average velocity of the three plates surrounding each of these junctions (Fig. 2d). Evolving oceanic spreading centers (ridges) clearly differ in their spreading rates, which produces distinct faulting patterns on faster (smoother pattern) and slower (rougher pattern) accreting parts of the plates (e.g., Püthe and Gerya, 2013). Pronounced elongated intra-plate ranges (called triple junction traces, e.g., Georgen and Lin, 2002), growing from triple junctions separate crust accreted from different spreading centers. The orientation of these intra-plate ranges may deviate slightly from the velocity vectors of the respective plates (Fig. 2d). The orientation of mid-ocean ridges forming between different pairs of plates seems to be closely defined by the bisectors of angles formed by respective intra-plate ranges. The differences in the plate accretion velocities around triple junctions define their asymmetric geometry, which deviates significantly (Fig. 2d) from ideal 120° junctions (e.g., McKenzie and Morgan, 1969). 3.1.2. Influence of model parameters We tested the influences of several model parameters on the nucleation and evolution of triple junctions, namely (Figs. 3–6): (i) the geometry and origin (thermal/rheological) of the initial perturbation, (ii) strain weakening rheology, (iii) the symmetry and (iv) magnitude of bi-directional extensional velocities. The geometry and origin of the initial perturbation mainly affects the initial rifting and breakup pattern of the plate (i.e., stages 1–3) but has little influence on the long-term character of the mature spreading pattern (i.e. for stage 4). Prescribing a larger 20 × 20 × 5 km thermal anomaly (Fig. 3) below the small 2 × 2 × 2 km weak inclusion (Fig. 2) produces a more regular cross-like geometry of the quadruple junction with an active central upwelling driven by temperature- and meltinginduced mantle buoyancy (Fig. 3b). Oceanic spreading initiates earlier by combined passive and active extension, which also results in an earlier breakup of the quadruple junction (cf. Figs. 2b,c and 3b,c). The subsequent spreading evolution is quite similar to the reference model except that the direction of the faster spreading ridge in the center of the model changes by 90° (cf. Figs. 2d and 3d). This direction is chosen randomly and varies non-systematically between different models (Figs. 2d, 3d, 4a,b, 5d). Prescribing of an elongated 20 × 98 × 5 km thermal perturbation results in a single rift formation with an elliptic mantle upwelling (Fig. 4a, 0.53 Ma). However, at a later stage, a standard pattern of two outward migrating triple junctions forms as the result of plate motion rearrangement during oceanic spreading (Fig. 4a, 5.76 Ma). Prescribing two 20 × 20 × 5 km offset thermal perturbations results in the rapid development of two linked irregular triple junctions formed by the plate rifting (Fig. 4b, 0.42 Ma). Subsequently, these rifting junctions are transformed by the oceanic plates' accretion into more regularly shaped outwards moving oceanic spreading junctions (Fig. 4b, 2.60 Ma). We additionally tested the stability of quadruple junctions by performing a simulation with two pre-existing oceanic ridges crossing at right angles (Fig. 5a). This “ideal” initial configuration remains relatively stable over ca. 2 Ma (Fig. 5b), producing four intra-plate ranges from the quadruple oceanic spreading junction. Eventually the quadruple junction becomes unstable and breaks into two diverging triple junctions separated by a lengthening faster spreading ridge (Fig. 5c,d). The breakup is associated with a fork-like splitting of two opposite intra-plate ranges (Fig. 5c). It is clearly visible from this model that the breakup of the quadruple junction and the subsequent divergence of two triple junctions gradually decreases the length of deforming boundaries between four diverging rigid oceanic plates (cf. Fig. 5b and d) and thus the integral mechanical resistance of these boundaries to the spreading.

3.2. Extension of one of the plates formed at a pre-existing ridge 3.2.1. Reference model development Fig. 7 shows the development of the reference model (trawc, Table 1) with an initial single ridge separating two plates. The right plate is subjected to bi-directional spreading with NS directed full spreading velocity of 1.9 cm/year, which is twice smaller than the WE spreading velocity (3.8 cm/year) imposed across the pre-existing ridge. This model is inspired by the results from the non-uniform bi-directional extension model, which developed a series of T-junctions at a single NS ridge (Fig. 6). In order to localize deformation in the middle of the right plate, cooling ages along the ridge decrease linearly from 0.5 Ma at the boundaries to 0.1 Ma in the center. Cooling ages at the left and right boundaries are 2 Ma and 4 Ma respectively. The evolution of this model can be subdivided into four stages: (1) rifting of the right plate and formation of an incipient T-junction (0–1.5 Ma, Fig. 7a), (2) onset of the slow spreading ridge in the rifting zone (1.5–2.5 Ma, Fig. 7b), (3) rightward migration of the T-junction associated with bending of the pre-existing faster spreading ridge (2.5–5.0 Ma, Fig. 7c), (4) steady state spreading with a nearly-stationary T-junction geometry and position (> 5.0 Ma, Fig. 7d). During the first stage, a ridge-orthogonal deep narrow rift forms within the right plate, which gives rise to an incipient T-junction at the crossing between the rift and the ridge (Fig. 7a). During the second stage, slow oceanic spreading starts inside the rift in the proximity of the faster spreading ridge. The ridge starts to bend gradually in response to a slow rightward migration of the T-junction (Fig. 7b). The third stage is manifested by continued ridge bending and T-junction migration (Fig. 7c). Slow oceanic spreading produces a triangular region with low elevation and rough topography characteristic for (ultra)slow spreading velocities (e.g. Püthe and Gerya, 2013). Intra-plate ranges separate this region from smoother crustal areas accreting at the faster spreading 8

Tectonophysics xxx (xxxx) xxx–xxx

T. Gerya, E. Burov

Fig. 3. Evolution of the numerical model with a central rectangular 20 × 20 km thermal anomaly (model trai in Table 1). Colors and indications are as in Fig. 2.

a)

b)

c)

d)

9

Tectonophysics xxx (xxxx) xxx–xxx

T. Gerya, E. Burov

a)

b)

c)

d)

Fig. 4. Influence of model parameters for the evolution of numerical experiments. (a) elongated thermal perturbation (model traj in Table 1). (b) two 20 × 20 km thermal perturbations (model trar in Table 1). (c) slower bi-directional extension rates for the case of central 10 × 10 km thermal perturbation (model trad in Table 1). (d) no strain weakening for the case of two pre-imposed crossing oceanic ridges (model tran in Table 1). Colors and indications are as in Fig. 2.

3.2.2. Influence of symmetry and magnitude of extension velocities We investigated the influence of extension velocities on the evolution of T-junctions (Figs. 8–10). Imposing uniform bi-directional extension on the second plate produces a gradual transformation of the Tjunction formed by the plate rifting (Fig. 8a) into a much more symmetrical triple junction migrating rightward (Fig. 8b–d). The initial straight NS ridge becomes strongly curved and a symmetrical (almost 120°) pattern of moving oceanic ridges and mid-plate ranges forms (Fig. 8d). This process again gradually decreases the length of deforming boundaries between three diverging rigid oceanic plates (cf. Fig. 8b and d) and thus the integral mechanical resistance of these

ridge. A relatively poorly pronounced intra-plate range starts to accrete from the T-junction within the left plate (Fig. 7c). The orientation of the ranges deviates slightly from the respective plate velocity vectors, whereas the orientation of mid-ocean ridges is nearly-parallel to the bisectors of angles between respective intra-plate ranges. During the last stage of steady state oceanic spreading, the position and geometry of the T-junction stabilizes and intra-plate ridges become prominent within all three plates. Similarly to the first reference model (Fig. 2), the differences in the plate accretion velocities around the T-junction define its asymmetric steady state geometry, which deviates significantly from classical 120° angles (Fig. 7d).

10

Tectonophysics xxx (xxxx) xxx–xxx

T. Gerya, E. Burov

Fig. 5. Evolution of the numerical model with two preimposed crossing oceanic ridges (model tram in Table 1). Colors and indications are as in Fig. 2.

a)

b)

c)

d)

11

Tectonophysics xxx (xxxx) xxx–xxx

T. Gerya, E. Burov

Fig. 6. Evolution of the numerical model with non-uniform bi-directional extension (model traf in Table 1). Colors and indications are as in Fig. 2.

a)

b)

c)

d)

12

Tectonophysics xxx (xxxx) xxx–xxx

T. Gerya, E. Burov

Fig. 7. Development of a migrating T-junction by non-uniform bi-directional extension of one of the plates formed at a pre-existing ridge (model trawc in Table 1). Colors and indications are as in Fig. 2.

a)

b)

c)

d)

13

Tectonophysics xxx (xxxx) xxx–xxx

T. Gerya, E. Burov

Fig. 8. Development of a migrating triple junction by uniform bi-directional extension of one of the plates formed at a pre-existing ridge (model trax in Table 1). Colors and indications are as in Fig. 2.

a)

b)

c)

d)

14

Tectonophysics xxx (xxxx) xxx–xxx

T. Gerya, E. Burov

Fig. 9. Development of a stationary 120° triple junction by balanced 120° bi-directional extension of one of the plates formed at a pre-existing ridge (model tray in Table 1). Colors and indications are as in Fig. 2.

a)

b)

c)

d)

15

Tectonophysics xxx (xxxx) xxx–xxx

T. Gerya, E. Burov

Fig. 10. Development of a migrating triple junction from pre-existing T-junction of two ridges (model trawb in Table 1). Colors and indications are as in Fig. 2.

a)

b)

c)

d)

16

Tectonophysics xxx (xxxx) xxx–xxx

T. Gerya, E. Burov

boundaries to the spreading. By comparing the final stationary shapes of the junctions in different numerical experiments (cf. Figs. 7d and 8d) we can also find a physical explanation for the stability of asymmetrical T-junctions compared to more symmetrical triple junctions in the case of strongly heterogeneous multi-directional spreading velocities (e.g., Fig. 7). The asymmetrical geometry of a T-junction predominantly increases the length of the slowest (and thus less dissipative) spreading segment and decreases the lengths of the two other faster (and thus more dissipative) segments. A stationary position and nearly ideal 120° geometry of the triple junction can be achieved by choosing balanced extension velocities, implying that the three plates diverge at 120° at the same rate (Fig. 9). Less pronounced mid-plate ridges and a coarser uniform topography pattern are produced due to the overall reduction in the spreading rates caused by reduced prescribed plate velocities (cf. Figs. 8 and 9, Table 1). Finally, we tested multi-directional plate migration models with a step-wise increase (doubling) of velocity in the middle of the right boundary and an initial T-junction formed by two pre-existing ridges (Fig. 10). In this model the T-junction is transformed into a migrating triple junction surrounded by three spreading ridges – all having dissimilar spreading rates and respectively faulting and topography patterns. The patterns change in the clockwise direction from very smooth to intermediate and then to coarse (Fig. 10c).

vectors and lengthening rates of oceanic ridges forming between the three pairs of plates v ′12 = {v x ′12 v y′12}, v ′13 = {v x ′13 v y′13} , v ′23 = {v x ′23 v y′23} are defined by the average relative velocity of respective plate pairs (Fig. 11b).

v ′12 = (v ′1 + v ′2 )/2, v x ′12 = (v x ′1 + v x ′2)/2, v y′12 = (v y′1 + v y′2)/2,

(18a)

v ′13 = (v ′1 + v ′3 )/2, v x ′13 = (v x ′1 + v x ′3)/2, v y′13 = (v y′1 + v y′3)/2,

(18b)

v ′23 = (v ′2 + v ′3 )/2, v x ′23 = (v x ′2 + v x ′3)/2, v y′23 = (v y′2 + v y′3)/2.

(18c)

This orientation ensures symmetrical accretion of material to both plates, which is, however, non-orthogonal to ridges and given by the following accretion vectors (half spreading rates):

vi : ij ≠ i = (vi − vj )/2, vxi : ij = (vxi − vxj )/2, vyi : ij = (vyi − vyj )/2,

(19)

where vi : ij ≠ i = {vxi : ij vyi : ij} is the accretion vector (half spreading rate) of the plate i from the ridge ij. The postulated ridge orientation is based on the results of our numerical experiments and (in the general case) introduces some shear component along the ridge. In nature, this component could be accommodated by ridge-transform patterns spontaneously forming at ridges that are oblique to the extension direction (e.g., Gerya, 2010b). The postulated ridge orientation must not be necessarily perpendicular to the spreading direction and thus differs from the one suggested by McKenzie and Morgan (1969). As one can see in Fig. 11b, the angles between three oceanic ridges crossing at a triple junction (as well as between intra-plate ranges) depend on the relative plate migration velocities v ′1 , v ′2 , v ′3 , and ideal 120° angles between all of them are only realized in the peculiar case when the relative velocities of all three migrating plates are equal in magnitude and oriented at 120° relative to each other. Therefore, geometrical differences between triple junctions and asymmetrical Tjunctions are mainly related to differences in their relative plate migration velocities such that any transitional geometries between the two classes of junctions are stable (McKenzie and Morgan, 1969) and can be potentially realized in nature. The pattern of ridges is always more symmetrical than that of intra-plate ranges (Fig. 11b) since orientations of the former are defined by averages between orientations of the later (Eqs. (18a), (18b), (18c)). The key question for our (and for any other) geometrical theory of a triple junction is how to express its migration velocity from the velocities of surrounding plates (Fig. 11a). In some numerical models (Fig. 2c,d, 7c,d) the triple junction migration direction clearly correlates with the average velocity of the three surrounding plates, however in other models (Fig. 6c,d, 10c,d) this relationship is less obvious. It should be noted that oceanic spreading triple junctions migrate in response to both the plates' motion and accretion. Therefore, the influence of both processes needs to be parameterized and the exact mathematical form of the function vt = f (v1 , v2 , v3) is not obvious. The theory should, however, account for two modeling-based observations: (i) the superior stability of triple junctions is related to the reduction of the cumulative length of ridges (e.g., Fig. 5) and (ii) the asymmetrical geometry of a T-junction predominantly increases the length of the slowest (and thus less dissipative) spreading segment and decreases the lengths of two other faster (and thus more dissipative) segments. (cf. Figs. 7d and 8d). Here, we propose that the triple junction migration velocity vt should satisfy the condition of minimal cumulative ridge lengthening rate (Rt) weighted by the square root of the full spreading rate

4. Discussion 4.1. Geometrical theory of a migrating steady-state triple junction Our numerical experiments suggest that the geometry of triple junctions varies widely from irregular T-junctions to ideal 120° junctions as a function of time and plate velocities. During the model development this geometry reaches a steady state even for migrating triple junctions (Figs. 3–5,7–10). The steady state geometry presumably minimizes the lengths of surrounding ridges (and thus dissipative work) and can vary from an ideal 120° shape to very asymmetric T-junctions depending on the symmetry of the applied plate velocities (cf. Figs. 9d, 7d). The structure of the numerically produced triple junctions includes two crucial tectonic elements; namely oceanic ridges (spreading centers) and intra-plate ranges (triple junction traces), which can be more or less pronounced depending on model parameters. Oceanic ridges are oriented nearly parallel to the bisectors of angles formed by intra-plate ranges. On the other hand, the orientation of intra-plate ranges can slightly but systematically deviate from the respective plate velocities. This systematic deviation is more pronounced for rapidly migrating triple junctions (Figs. 2d, 3d, 5d, 8d, 10d). Based on these numerical results, we propose a new geometrical theory of migrating steady state triple junctions (Fig. 11). In an arbitrary (e.g., immobile) reference frame, the steady state of any triple junction (as well as of any T-junction) can be fully described by the migration velocity vector of the junction vt = {vxt , vyt } and three velocity vectors of interacting moving plates v1 = {vx1 , vy1} , v2 = {vx 2 , vy2} , v3 = {vx 3 , vy3} (Fig. 11a). It is then convenient to consider the moving triple junction reference frame x’, y’ (Fig. 11b). The internal geometrical steady state of the junction in this frame is independent on its position and migration velocity and is fully described by the relative velocity vectors of three interacting plates v ′1 = {vx ′ 1 , vy ′ 1} , v ′2 = {vx ′ 2 , vy ′ 2} , v ′3 = {vx ′ 3 , vy ′ 3} which are defined as (Fig. 11b).

v ′1 = v1 − vt , v x ′1 = vx1 − vxt , v y′1 = vy1 − vyt ,

(17a)

v ′2 = v2 − vt , v x ′2 = vx 2 − vxt , v y′2 = vy2 − vyt ,

(17b)

v x ′13 + v y′13 ×

4

4vx1:132 + 4vy1:132 +

v ′3 = v3 − vt , v x ′3 = vx 3 − vxt , v y′3 = vy3 − vyt .

(17c)

v x ′232 + v y′232 ×

4

4vx 2:232 + 4vy2:232 = min.

v x ′122 + v y′122 ×

Rt = 2

Intra-plate ridges (i.e., triple junction traces) accrete directly from the triple junction and therefore their orientation and accretion rate vectors match the respective relative plate velocity vectors. Orientation

2

4

4vx1:122 + 4vy1:122 +

(20)

Weighting of the ridge lengthening rate by the square root of the spreading rate has been found empirically by comparing several formalisms to both numerical results and (rather limited) natural data. 17

Tectonophysics xxx (xxxx) xxx–xxx

T. Gerya, E. Burov

a)

Fig. 11. (a), (b) geometrical theory of a migrating steady state triple junction. (c) - comparison of numerically estimated and calculated (Eq. (20)) vxt (orange squares) and vyt (green triangles) triple junction migration velocity components. (d), (e), (f) and (g) – comparison of computed velocity vectors (Eqs. (17a), (17b), (17c), (18a), (18b), (18c), (20)) with orientations of intra-plate ranges and oceanic ridges for models trai, trawc, traf and trawb (Table 1, Figs. 3d, 7d, 6d, 10d), respectively. Triple junction migration vectors, relative plate velocity vectors, ridge lengthening vectors, and accretion vectors are shown by red, blue, black and green arrows, respectively. See text for discussion. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

b)

c)

d)

e)

f)

g)

of respective intra-plate ranges and oceanic ridges (Fig. 11d–g). The surprisingly good agreement of Eq. (20) with numerical experiments exploring non-linear multi-physics spreading models implies that the proposed simple formalism captures the dominant mechanisms governing plate migration and accretion process at oceanic spreading triple junctions. Our geometrical theory differs significantly from the theory of McKenzie and Morgan (1969) who assumed that ridges spread symmetrically at right angles to their strike. Their orientations then correspond to perpendicular bisectors of the sides of spreading velocity triangle (cf. green dashed triangle in Fig. 11a). The perpendicular bisectors meet at

This comparison showed that Eq. (20) allows for notably better agreement with both sources of data than a non-weighted model. The proposed geometrical theory compares well with the results of numerical experiments for advanced stages of triple junctions migration (Fig. 11c–g). Using Eq. (20), we computed (Table 1) long-term migration velocity vectors and weighted minimal ridge lengthening rates for all documented triple junctions and T-junctions as a function of the surrounding plate velocities. Computed triple junction migration velocities agree well with their values evaluated from the numerical experiments (Fig. 11c). Computed relative plate velocities and ridge lengthening vectors are in a very good agreement with the orientation 18

Tectonophysics xxx (xxxx) xxx–xxx

T. Gerya, E. Burov

a point called the centroid, and this point in velocity space gives the triple junction migration velocity (McKenzie and Morgan, 1969). It should be noted that the geometrical theory of McKenzie and Morgan (1969) fails to correctly predict triple junction geometries and migration velocities obtained in our numerical experiments. This theory in particular erroneously suggests that no migration of triple junctions should occur in our reference model shown in Fig. 3. It is possible to demonstrate that the proposed geometrical model is physically meaningful. Weighting of the ridge lengthening rate by the square root of the spreading rate is equivalent to weighting by energy dissipation (Hs) per unit ridge length under assumption of constant brittle strength (σy) of deforming ridge crust.

Hs = 2σy vs

1 L

∫0

L

z b (l) dl,

a)

b)

(21)

where zb(l) is thickness of the deforming brittle crust at the ridge, vs is half spreading rate (Eq. (19)), L is the half-width of the deformation zone and the l coordinate is taken along one of the spreading directions at the ridge. According to the half-space cooling model (e.g. Turcotte and Schubert, 2002).

c)

z (l) ⎞ ⎛ z b (l) vs ⎞, Tb = Ta erf ⎜⎛ b ⎟ = Ta erf 2 κt ( l ) ⎝ 2 κl ⎠ ⎝ ⎠ ⎜

⎟

(22)

where x is horizontal distance from the ridge, t(l) = l/vs is cooling age, Tb is the temperature of the brittle/ductile transition, Ta is asthenospheric temperature and κ is thermal diffusivity. Combining Eqs. (21) and (22) gives.

Hs = 2σy vs

1 L

∫0

2erf −1 (Tb/ Ta ) κl 8 T dl = σyerf −1 ⎛ b ⎞ κLvs , vs 3 ⎝ Ta ⎠

L

⎜

⎟

(23)

which thus implies that the energy dissipation per unit ridge length is indeed proportional to the square root of the spreading rate vs (cf. Eq. (20)). The physical meaning of the weighting of ridge lengthening rate by Hs (respectively vs ) is that less dissipative slower ridges should grow proportionally faster. In order to test these theoretical considerations we compared (Fig. 12) the evolution of total energy dissipation with changes in the length of ridges for the model shown in Fig. 4d. This comparison shows that both energy dissipation and the cumulative length of ridges gradually decrease after splitting of a quadruple junction into two outward moving triple junctions (Figs. 4d, 12a,b). Moreover, energy dissipation during this period of model evolution (2–12 Ma) linearly correlates with the sum of lengths of ridges multiplied to square root of their spreading rates (Fig. 12c), which strongly supports our theory. In accordance with our theory, the number of intra-plate ranges (i.e., triple junction traces) in a single plate can be more than one and is equal to the number of triple junctions formed at that plate (Fig. 11d,f). Importantly, the orientation of different ranges in the same plate can also differ in accordance with differences in migration velocities of respective triple junctions (Fig. 11d,f). It should be stressed, however, that our geometrical theory applies to triple junctions in steady state, whereas the geometry of immature junctions changes with time and could deviate from steady state conditions (cf. Fig. 7b and d). The theory can be generalized for a multi-branch steady state oceanic spreading junction of an arbitrary dimension n > 1 (2 = ridge, 3 = triple junction, 4 = quadruple junction etc.) n−1

Rm =

n

∑ ∑ i=1 j=i+1

×

4

⎛ ⎝

vxi + vxj 2

2

− vxm ⎞ + ⎛ ⎠ ⎝

vyi + vyj

(vxi − vxj )2 + (vyi − vyj )2 = min,

2

Fig. 12. (a) Evolution of model topography after splitting of a quadruple junction into two outward moving triple junctions for the numerical experiment without strain weakening (2–12 Ma, model tran in Table 1, Fig.4d). (b) evolution of total length of ridges (solid squares) and energy dissipation (open triangles) with time. (c) linear correlation between the total weighted length of ridges [sum of lengths of ridges (km) multiplied to square root of their full spreading rates (km/Ma)] with energy dissipation.

also be used to analyze the relative stability of different configurations of ridges and junctions. For example, the solution of Eq. (24) for the quadruple junction shown in Fig. 5b (vxm = 0 cm/yr, vym = 0 cm/yr, Rm = 14.70 cm3/2/yr3/2) gives a larger weighted integral ridge lengthening rate than the sum of rates computed for two diverging triple junctions shown in Fig. 5c (vxm1 = − 0.246 cm/yr, vym1 = 0.246 cm/yr, Rm1 = 7.27 cm3/2/yr3/2 for junction 1 and vxm2 = 0.246 cm/yr, vym2 = −0.246 cm/yr, Rm2 = 7.27 cm3/2/yr3/2 for junction 2, Rm1 + Rm2 = = 14.54 cm3/2/yr3/2). This result explains the stationary position of the quadruple junction (Fig. 5a,b) and its breakup (splitting) into two migrating triple junctions (Fig. 5c,d) thus confirming the intrinsic instability (metastability) of quadruple junctions observed in the numerical experiments (Figs. 2, 3, 4c,d, 5). Estimates of the cumulative ridge lengthening rate Rm can also be used to evaluate the characteristic timescale (tm) of junction geometry reconfiguration toward steady state.

2

tm = La / Rm 2/3,

− vym ⎞ ⎠

(25)

where La is a characteristic ridge lengthening distance required for geometry reconfiguration. Based on the results of the numerical experiments (e.g., Fig. 8a–c) this length is on the order of 20–50 km, which implies a reconfiguration timescale of 1.5–8 Ma computed from estimated Rm values (2.6–7.3 cm3/2/yr3/2, Table 1).

(24)

where Rm is the cumulative weighted ridge lengthening rate for the junction, vm = {vxm , vym} is the junction migration velocity vector, vi = {vxi , vyi} , vj = {vxj vyj} are the velocity vectors of n interacting plates and summation is done over all potential plate boundaries. Eq. (24) can 19

Tectonophysics xxx (xxxx) xxx–xxx

T. Gerya, E. Burov

Fig. 13. Comparison of (a) the Rodriguez Triple Junction (Nakamura et al., 2007) with (b) numerical modeling results (model trawd in Table 1, 5.35 Ma). Numerical model is rotated by 180° for convenience. Comparison of schematic large-scale geometry (Georgen and Lin, 2002) of Rodriguez (c), Azores (d) and Galapagos (e) triple junctions with computed results from the geometrical model (Eq. (20)) based on relative plate velocities (black arrows, Georgen and Lin, 2002). Triple junction migration vectors, relative plate velocity vectors and ridge lengthening vectors are shown by red, blue and orange arrows, respectively. Abbreviations: RTJ = Rodriguez Triple Junction, CIR = Central Indian Ridge. SWIR – Southwestern Indian Ridge, SEIR – Southeastern Indian Ridge, N. MAR - Mid-Atlantic Ridge to the north of the Azores Triple Junction, S. MAR - MidAtlantic Ridge to the south of the Azores Triple Junction, Ter. R. - Terceira Rift, N. EPR - East Pacific Rise to the north of the Galapagos Triple Junction, S. EPR - East Pacific Rise to the south of the Galapagos Triple Junction, GSC = Galapagos Spreading Center, and AAR = American Antarctic Ridge. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

c)

a)

d)

b) e)

lithosphere and mantle plumes, a more specialized numerical modeling effort is needed to understand the initiation and dynamics of continental breakup patterns (e.g., Brune and Autin, 2013; Brune et al., 2014; Burov and Gerya, 2014; Koptev et al., 2015, 2016) including quadruple and triple rifting/breakup junctions. We will, therefore, focus our comparison with natural data on oceanic spreading junctions. The detailed comparison between the numerical results and natural data is limited by the availability of high-resolution bathymetry and magnetic data from the regions of modern triple junctions (e.g., Honsho et al., 1996; Ligi et al., 1999; Mitchell et al., 2011; Smith et al., 2011). Fig. 13a,b presents a first-order comparison for the Rodriguez Triple Junction in the Indian Ocean (e.g. Tapscott et al., 1980; Mitchell, 1991; Mitchell and Parson, 1993; Dyment, 1993; Honsho et al., 1996; Nakamura et al., 2007). As can be seen from this figure, models and data compare well in terms of the major structural elements. As in the numerical model, the (ultra)slow Southwestern Indian Ridge domain has a narrow triangular shape, coarser and deeper topography and is surrounded by pronounced intra-plate ridges (triple junction traces) separating coarse (ultra)slow spreading pattern of the Southwestern Indian Ridge from much smoother intermediate spreading patterns of the Southeastern Indian Ridge and Central Indian Ridge, which are oriented at a gentle angle to each other. Similarly to the model, no clear

4.2. Comparison to natural data Numerical experiments predict that (oceanic) plate rifting patterns under bi-directional extension can be represented by linear rifts (unbalanced extension and/or linear perturbations), quadruple junctions (balanced extension with point-like perturbations) and triple junctions (balanced extension with multiple perturbations). All these rifting patterns are transient and are strongly reworked and/or fully destabilized on the timescale of a few million years by subsequent oceanic spreading and crust accretion, which transforms them into either migrating or stationary triple junctions. Our experiments thus seem to support previous suggestions that R-R-R-R quadruple junctions can be formed by continental breakup (not modeled here) and exist for some time during oceanic spreading (Hey and Milholland, 1979). On the other hand, the results also suggest that the governing physics of plate rifting/breakup junctions (fragmentation structures) and oceanic spreading junctions (accretion structures) are distinctly different and they should be treated as separate classes of structures (e.g., Gerya, 2013). Our geometrical theory is developed for the stationary oceanic triple junctions and does not apply to transient plate breakup junctions (such as the Afar triple junction). Given the relatively small size and simplified geometry of our models, which do not include continental 20

Tectonophysics xxx (xxxx) xxx–xxx

T. Gerya, E. Burov

Fig. 14. Conceptual scheme of bi-directional lithospheric extension showing transition from continental breakup and formation of quadruple plate rifting junction to mature oceanic spreading with diverging triple junctions.

typically have a pronounced elevated topography (Bird et al., 1999; Larson et al., 2002; Pockalny et al., 2002; Georgen and Lin, 2002). A comparison of our geometrical theory with the large-scale geometry of triple junctions (Georgen and Lin, 2002) shows a relatively good agreement (Fig. 13c–e). Larger discrepancies occur for the Azores, but the steady state of this junction is not guaranteed because of the

elongated intra-plate ridge is visible to the east of the triple junction where only some broad region of positive topography is observed. It should be noted, however, that both natural evidence and physical controls for intra-plate ridges should be investigated more thoroughly. In contrast to the models, these tectonic elements (called triple junction traces) could be more difficult to recognize in nature and they do not

21

Tectonophysics xxx (xxxx) xxx–xxx

T. Gerya, E. Burov

spreading junctions of an arbitrary dimension n > 1 (2 = ridges, 3 = triple junctions, 4 = quadruple junctions etc.). The resulting equation can also be used to analyze the relative stability of different configurations of ridges and junctions. It shows, in particular, that a quadruple junction has a larger weighted integral ridge lengthening rate than the sum of rates for two diverging triple junctions. This explains the instability and breakup (metastability) of the quadruple junctions observed in numerical experiments. 7. Predictions from numerical models find some good analogies in nature but more systematic cross-disciplinary work is needed in the future, which should combine new high-resolution bathymetry and magnetic data in the regions of modern triple junctions and systematic numerical and analogue modeling.

very low spreading rates in the Terceira Rift (Fig. 13d). It should also be noted that triple junctions should only influence orientations of ridge portions that lengthen from them with time. The Mid-Atlantic Ridge to the north of the Azores should have lengthened very slowly (if at all) from the triple junction (Fig. 13d) and is thus likely in a transient state. It should also be noted that the steady state of other modern triple junctions is not guaranteed, although characteristic timescales of reaching the steady state should be relatively short (on the order of few Ma, Eq. (25)). 5. Conclusions Here, we investigate the nucleation and evolution of oceanic triple junctions by using high-resolution 3D numerical magmatic-thermomechanical experiments that take into account a realistic thermorheological structure of the oceanic lithosphere as well as crustal and lithospheric accretion processes. The following conclusions can be made on the basis of our study:

Acknowledgements This paper is based on our collaborative work with my colleague and friend Evgueni Burov who passed away unexpectedly 9 October 2015. We performed numerical experiments in 2014 and started writing the paper in 2015. The paper has been finalized recently by T.G. for this special Tectonophysics volume dedicated to Evgueni Burov who was particularly enthusiastic about the numerical modeling and understanding of triple junctions. This work was supported by an ERC ITNgrants 604713 ZIP, 642029 CREEP and 674899 SUBITOP (T.G.), by a SNF-project CRSII2_154434 Swiss-AlpArray (T.G.), an Advanced ERC grant 290864 RHEOLITH (E.B.), by INSU-CNRS, by a UPMC Invited Professor grant (T.G.) and by an ETH Invited Professor grant (E.B.). Numerical simulations were performed on the ETH Brutus and Euler clusters. Comments and suggestions from Boris Kaus, an anonymous reviewer, Ryan Stoner, Philippe Agard and Fernando Ornelas Marques, are greatly appreciated.

1. Two major classes of junctions form under a bi-directional or multidirectional extensional velocity field (Fig. 14): (i) transient quadruple and triple plate rifting junctions form from local perturbations by initial plate fragmentation and are subsequently converted into (ii) stable triple oceanic spreading junctions that are controlled by the addition of new oceanic crust and the accretion of mantle lithosphere. 2. Quadruple R-R-R-R plate rifting junctions are transient features that rapidly (within few Ma) evolve toward two diverging triple oceanic spreading junctions connected by a linear spreading center lengthening with time. This process gradually decreases the cumulative length of deforming boundaries between four diverging rigid oceanic plates and thus the integral mechanical resistance of these boundaries to spreading. The superior stability of triple oceanic spreading junctions thus satisfies the principle of minimum mechanical dissipative work. 3. The geometry of triple oceanic spreading junctions varies as a function of time and plate velocities. A geometrical steady state can be achieved on a timescale of several million years for both stationary and migrating triple junctions. The steady state geometry can vary from an ideal 120° shape to very asymmetric T-junctions depending on the symmetry of applied plate velocities. In the case of strongly heterogeneous multi-directional spreading velocities, the asymmetrical geometry of a T-junction predominantly increases the length of the slowest (and thus less dissipative) spreading segment and decreases the lengths of the two other faster (and thus more dissipative) segments. 4. The structure of the numerically produced oceanic triple junctions includes two crucial tectonic elements, which are oceanic ridges (spreading centers) and intra-plate ranges (triple junction traces), which can be more or less pronounced depending on the spreading parameters. Oceanic ridges in numerical experiments are oriented nearly-parallel to bisectors of angles formed by intra-plate ranges. 5. The new geometrical theory of a migrating steady state triple junction is proposed, which suggests that the steady state of the junction can be fully described by the migration velocity of the junction and velocity vectors of three interacting plates. The geometry of the junction is defined by the relative (i.e. in the moving frame of the migrating triple junctions) plate velocities, which define the orientations and growth rates of respective intra-plate ranges. The orientations and lengthening rates of oceanic ridges are in turn given by the average relative velocity of two accreting plates. The steady state triple junction migration velocity satisfies the condition of minimal cumulative ridge lengthening rate weighted by the square root of the spreading velocity. The analytical theory is in good agreement with the numerical experiments and natural data. 6. The theory is generalized to multi-branch steady state oceanic

References Allken, V., Huismans, R.S., Thieulot, C., 2011. Three-dimensional numerical modeling of upper crustal extensional systems. J. Geophys. Res. 116, B10409. Allken, V., Huismans, R.S., Thieulot, C., 2012. Factors controlling the mode of rift interaction in brittle-ductile coupled systems: a 3D numerical study. Geochem. Geophys. Geosyst. 13, Q05010. Apotria, T.G., Gray, N.H., 1985. Absolute motion and evolution of the Bouvet triple junction. Nature 316, 623–625. Apotria, T.G., Gray, N.H., 1988. The evolution of the Bouvet triple junction: implications of its absolute motion. Tectonophysics 148, 177–193. Beutel, E., van Wijk, J., Ebinger, C., Keir, D., Agostini, A., 2010. Formation and stability of magmatic segments in the Main Ethiopian and Afar rifts. Earth Planet. Sci. Lett. 293, 225–235. Bird, R.T., Tebbens, S.F., Kleinrock, M.C., Naar, D.F., 1999. Episodic triple-junction migration by rift propagation and microplates. Geology 27, 911–914. Brune, S., Autin, J., 2013. The rift to break-up evolution of the Gulf of Aden: insights from 3D numerical lithosphere-scale modelling. Tectonophysics 60, 65–79. Brune, S., Heine, C., Pérez-Gussinye, M., Sobolev, S., 2014. Rift migration explains continental margin asymmetry and crustal hyper-extension. Nat. Commun. 5, 4014. Buck, W.R., Lavier, L.L., Poliakov, A.N.B., 2005. Modes of faulting at mid-ocean ridges. Nature 434, 719–723. Burke, K., Dewey, J.F., 1973. Plume-generated triple junctions: key indicators in applying plate tectonics to old rocks. J. Geol. 81, 406–433. Burov, E., Gerya, T., 2014. Asymmetric three-dimensional topography over mantle plumes. Nature 513, 85–89. Choi, E., Lavier, L., Gurnis, M., 2008. Thermomechanics of mid-ocean ridge segmentation. Phys. Earth Planet. Inter. 171, 374–386. Clauser, C., Huenges, E., 1995. Thermal conductivity of rocks and minerals. In: Ahrens, T.J. (Ed.), Rock Physics and Phase Relations. AGU Reference Shelf. 3. American Geophysical Union, Washington D.C, pp. 105–126. Connolly, J.A.D., Schmidt, M.W., Solferino, G., Bagdassarov, N., 2009. Permeability of asthenospheric mantle and melt extraction rates at mid-ocean ridges. Nature 462, 209–212. Crameri, F., Schmeling, H., Golabek, G.J., Duretz, T., Orendt, R., Buiter, S.J.H., May, D.A., Kaus, B.J.P., Gerya, T.V., Tackley, P.J., 2012. A comparison of numerical surface topography calculations in geodynamic modelling: an evaluation of the ‘sticky air’ method. Geophys. J. Int. 189, 38–54. Dordevic, M., Georgen, J., 2016. Dynamics of plume-triple junction interaction: results from a series of three-dimensional numerical models and implications for the formation of oceanic plateaus. J. Geophys. Res. 121, 1316–1342. Dyment, J., 1993. Evolution of the Indian Ocean Triple Junction between 65 and 49 Ma (anomalies 28 to 21). J. Geophys. Res. 98, 13863–13877.

22

Tectonophysics xxx (xxxx) xxx–xxx

T. Gerya, E. Burov

Ligi, M., Bonatti, E., Bortoluzzi, G., Carrara, G., Fabretti, P., Penitenti, D., Gilod, D., Peyve, A.A., Skolotnev, S., Turko, N., 1997. Death and transfiguration of a triple junction in the South Atlantic. Science 276, 243–245. Ligi, M., Bonatti, E., Bortoluzzi, G., Carrara, G., Fabretti, P., Gilod, D., Peyve, A.A., Skolotnev, S., Turko, N., 1999. Bouvet Triple Junction in the South Atlantic: geology and evolution. J. Geophys. Res. 104, 29365–29385. McClusky, S., Reilinger, R., Mahmoud, S., Sari, D.B., Tealeb, A., 2003. GPS constraints on Africa (Nubia) and Arabia plate motions. Geophys. J. Int. 155, 126–138. McKenzie, D.P., Morgan, W., 1969. Evolution of triple junctions. Nature 224, 125–133. McKenzie, D.P., Davies, D., Molnar, P., 1970. Plate tectonics of the Red Sea and East Africa. Nature 226, 243–248. Mitchell, N.C., 1991. Distributed extension at the Indian Ocean triple junction. J. Geophys. Res. 96, 8019–8043. Mitchell, N.C., Parson, L.M., 1993. The tectonic evolution of the Indian Ocean Triple Junction, anomaly 6 to present. J. Geophys. Res. 98, 1793–1812. Mitchell, G.A., Montésia, L.G.J., Zhu, W., Smith, D.K., Schouten, H., 2011. Transient rifting north of the Galápagos Triple Junction. Earth Planet. Sci. Lett. 307, 461–469. Nakamura, K., Kato, Y., Tamaki, K., Ishii, T., 2007. Geochemistry of hydrothermally altered basaltic rocks from the Southwest Indian Ridge near the Rodriguez Triple Junction. Mar. Geol. 239, 125–141. Nakanishi, M., Sager, W.W., Klaus, A., 1999. Magnetic lineations within Shatsky Rise, northwest Pacific Ocean: implications for hot spot-triple junction interaction and oceanic plateau formation. J. Geophys. Res. 104, 7539–7556. Nikolaeva, K., Gerya, T.V., Connolly, J.A.D., 2008. Numerical modelling of crustal growth in intraoceanic volcanic arcs. Phys. Earth Planet. Inter. 171, 336–356. Patriat, P., Courtillot, V., 1984. On the stability of triple junctions and its relation to episodicity in spreading. Tectonics 3, 317–332. Patriat, P., Parson, L., 1989. A survey of the Indian ocean triple junction trace within the Antarctic plate. Implications for the junction evolution since 15 Ma. Mar. Geophys. Res. 11, 89–100. Pockalny, R.A., Larson, R.L., Viso, R.F., Abrams, L.J., 2002. Bathymetry and gravity data across a mid-Cretaceous triple junction trace in the southwest Pacific basin. Geophys. Res. Lett. 29. http://dx.doi.org/10.1029/2001GL013517. Püthe, C., Gerya, T.V., 2013. Dependence of mid-ocean ridge morphology on spreading rate in numerical 3-D models. Gondwana Res. 25, 270–283. Ranalli, G., 1995. Rheology of the Earth. Chapman & Hall, London, UK. Schmeling, H., Babeyko, A.Y., Enns, A., Faccenna, C., Funiciello, F., Gerya, T., Golabek, G.J., Grigull, S., Kaus, B.J.P., Morra, G., Schmalholz, S.M., van Hunen, J., 2008. Benchmark comparison of spontaneous subduction models - towards a free surface. Phys. Earth Planet. Inter. 171, 198–223. Smith, D.K., Schouten, H., Zhu, W.L., Montesi, L.G.J., Cann, J.R., 2011. Distributed deformation ahead of the Cocos-Nazca Rift at the Galapagos triple junction. Geochem. Geophys. Geosyst. 12, Q11003. Tapscott, C., Patriat, P., Fisher, R.L., Sclater, J.G., Hoskins, H., Parsons, B., 1980. The Indian Ocean triple junction. J. Geophys. Res. 85, 4723–4739. Theissen-Krah, S., Iyer, K., Rüpke, L.H., Phipps Morgan, J., 2011. Coupled mechanical and hydrothermal modeling of crustal accretion at intermediate to fast spreading ridges. Earth Planet. Sci. Lett. 311, 275–286. Turcotte, D.L., Schubert, G., 2002. Geodynamics. Cambridge University Press, Cambridge. Verzhbitskii, E.V., Lobkovskii, I., Kononov, V., Kotelkin, D., 2006. Genesis of Shatsky and Hess oceanic rises in the Pacific Ocean as deduced from geologic-geophysical data and numerical modeling. Geotectonics 40, 236–245. Wanless, V.D., Shaw, A.M., 2012. Lower crustal crystallization and melt evolution at midocean ridges. Nat. Geosci. 5, 651–655.

Fournier, M., Patriat, P., Leroy, S., 2001. Reappraisal of the Arabia-India-Somalia triple junction kinematics. Earth Planet. Sci. Lett. 189, 103–114. Georgen, J.E., 2008. Mantle flow and melting beneath oceanic ridge–ridge–ridge triple junctions. Earth Planet. Sci. Lett. 270, 231–240. Georgen, J.E., Lin, J., 2002. Three-dimensional passive flow and temperature structure beneath oceanic ridge–ridge–ridge triple junctions. Earth Planet. Sci. Lett. 204, 115–132. Georgen, J.E., Sankar, R.D., 2010. Effects of ridge geometry on mantle dynamics in an oceanic triple junction region: implications for the Azores Plateau. Earth Planet. Sci. Lett. 298, 23–34. Gerya, T.V., 2010a. Introduction to Numerical Geodynamic Modelling. Cambridge University Press, Cambridge, UK. Gerya, T., 2010b. Dynamical instability produces transform faults at mid-ocean ridges. Science 329, 1047–1050. Gerya, T.V., 2013. Three-dimensional thermomechanical modeling of oceanic spreading initiation and evolution. Phys. Earth Planet. Inter. 214, 35–52. Gerya, T.V., 2014. Plume-induced crustal convection: 3D thermomechanical model and implications for the origin of novae and coronae on Venus. Earth Planet. Sci. Lett. 391, 183–192. Gerya, T.V., Yuen, D.A., 2007. Robust characteristics method for modelling multiphase visco-elasto-plastic thermo-mechanical problems. Phys. Earth Planet. Inter. 163, 83–105. Gerya, T.V., Stern, R.J., Baes, M., Sobolev, S., Whattam, S.A., 2015. Plate tectonics on the Earth triggered by plume-induced subduction initiation. Nature 527, 221–225. Gregg, P.M., Behn, M.D., Lin, J., Grove, T.L., 2009. Melt generation, crystallization, and extraction beneath segmented oceanic transform faults. J. Geophys. Res. 114, B11102. Hess, P.C., 1989. Origin of Igneous Rocks. Harward University Press, London, UK. Hey, R., Milholland, P., 1979. Stability of quadruple junctions. Nature 277, 201–202. Honsho, C., Tamaki, K., Fujimoto, H., 1996. Three-dimensional magnetic and gravity studies of the Rodriguez Triple Junction in the Indian Ocean. J. Geophys. Res. 10, 15837–15848. Huismans, R.S., Beaumont, C., 2002. Asymmetric lithospheric extension: the role of frictional plastic strain softening inferred from numerical experiments. Geology 30, 211–214. Katz, R.F., 2010. Porosity-driven convection and asymmetry beneath mid-ocean ridges. Geochem. Geophys. Geosyst. 11, Q0AC07. Katz, R.F., Spiegelman, M., Langmuir, C.H., 2003. A new parameterization of hydrous mantle melting. Geochem. Geophys. Geosyst. 4, 1073. Kleinrock, M.C., Morgan, J.P., 1988. Triple junction reorganization. J. Geophys. Res. 93, 2981–2996. Koptev, A., Calais, E., Burov, E., Leroy, S., Gerya, T., 2015. Dual continental rift systems generated by plume-lithosphere interaction. Nat. Geosci. 8, 388–392. Koptev, A., Burov, E., Calais, E., Leroy, S., Gerya, T., Guillou-Frottier, L., Cloetingh, S., 2016. Contrasted continental rifting via plume-craton interaction: applications to Central East African Rift. Geosci. Front. 7, 221–236. Kotelkin, D., Lobkovskii, I., Verzhbitskii, E.V., Kononov, V., 2004. A geodynamical model for the formation of the Shatsky Rise (Pacific Ocean). Oceanology 44, 257–260. Larson, R.L., Pockalny, R.A., Viso, R.F., Erba, E., Abrams, L.J., Luyendyk, B.P., Stock, J.M., Clayton, R.W., 2002. Mid-Cretaceous tectonic evolution of the Tongareva triple junction in the southwestern Pacific Basin. Geology 30, 67–70. Lavier, L.L., Buck, W.R., Poliakov, A.N.B., 2000. Factors controlling normal fault offset in an ideal brittle layer. J. Geophys. Res. 105, 23431–23442. Le Pichon, X.T., Gaulier, J.M., 1988. The rotation of Arabia and the Levant fault system. Tectonophysics 153, 271–294.

23