Afar triple junction triggered by plume-assisted bi-directional continental break-up

Divergent ridge-ridge-ridge (R-R-R) triple junctions are one of the most remarkable, yet largely enigmatic, features of plate tectonics. The juncture of the Arabian, Nubian, and Somalian plates is a type-example of the early development stage of a triple junction where three active rifts meet at a ‘triple point’ in Central Afar. This structure may result from the impingement of the Afar plume into a non-uniformly stressed continental lithosphere, but this process has never been reproduced by self-consistent plume-lithosphere interaction experiments. Here we use 3D thermo-mechanical numerical models to examine the initiation of plume-induced rift systems under variable far-field stress conditions. Whereas simple linear rift structures are preferred under uni-directional extension, we find that more complex patterns form in response to bi-directional extension, combining one or several R-R-R triple junctions. These triple junctions optimize the geometry of continental break-up by minimizing the amount of dissipative mechanical work required to accommodate multi-directional extension. Our models suggest that Afar-like triple junctions are an end-member mode of plume-induced bi-directional rifting that combines asymmetrical northward pull and symmetrical EW extension at similar rates.

slab-pull force that developed along the subduction boundary between the Eurasian and Afro-Arabian plates 18 . The Mediterranean and Biltis segments of that plate boundary progressively transitioned to continental collision between 30 and 20 Ma, with a slow-down of the Africa-Eurasia convergence 19 whereas active subduction continued along the remaining eastern segment of the Neo-Tethys 20 . The eastward closing of that ocean 21 led to a gradual westward decrease of slab pull forces from the Makran subduction to the collisional segment (Fig. 1b). The early collision stage coincides with the eruption of voluminous flood basalts linked to the emplacement of the Afar plume at ~30 Ma 22 . Thus, the Africa-Arabia separation was likely favoured not only by appropriate intraplate stress generated by lateral variation of slab-pull force along the Neo-Tethys subduction, but also by the thermal weakening effect of the plume impingement that facilitated strain localization inside the Afro-Arabian plate 18 .
At the same time, anomalous positive topography started to form in East Africa 23 , with the first evidence of doming in Ethiopia at ~30-40 Ma likely related to the Afar plume emplacement 22 . The current long-wavelength elevation of eastern and southern Africa, dynamically supported by the African Superplume, a broad low-velocity seismic anomaly imaged in the lower mantle 24,25 , is acquired at ~10 Ma 23 although the exact timing and dynamics of uplift in the Horn of Africa are still controversial [26][27][28] . The resulting lateral gradients of gravitational potential energy generate EW-directed extensional deviatoric stresses corresponding to forces on the same order as slab pull forces when integrated over the thickness of the lithosphere 29 . Given that the African continent is surrounded mostly by oceanic ridges, far-field extension can only be driven by a combination of lateral gradients of gravitational potential energy and viscous coupling with horizontal mantle flow at the base of the lithosphere [29][30][31] . Recent results of spherical shell finite element modelling of regional stress and strain field in Africa has shown that compressional stresses applied at the mid-ocean ridges surrounding Africa result in rift-perpendicular deviatoric extension in the East African rift valleys that, being combined with upwelling mantle plume(s), leads to localization of extensional deformation and explains the initial lithosphere break-up in the East Africa and the Afar 32 . Likewise, recent 3D thermo-mechanical deformation models show that pre-stressed lithosphere subjected to EW extension is a necessary initial condition for plume-induced localized rifts to develop in the central part of the East African Rift 8,9,33 .
We note that the detailed reproduction of the the timing of triple junction development is beyond the scope of this paper. We do not target to include such aspects as a delayed rifting in Northern Ethiopia with respect to that in the Red Sea and Gulf of Aden [34][35][36][37][38] or a westward propagation of the Gulf of Aden opening in the Miocene due to lateral propagation of the Sheba Ridge from the Indian Ocean into the African continent 15,39,40 . In contrast, our main objective is to quantify the general consequences of thermal and buoyancy-driven mechanical effects of the plume head in the context of laterally homogenous lithosphere subjected to bi-directional, asymmetric far-field stresses. To this purpose we first present the generic features of the experiments, followed by a comparison of our model inferences with observed present-day configuration of break-up zones in the Afar triple junction.

Modelling Approach and Results
The Afar triple junction forms in a context where a plume interacts with a continental lithosphere subjected to bi-directional far-field forcing that combines northward pull from the Neo-Tethys slab subduction and EW extension. We use this configuration to design a model geometry that we implement into high-resolution 3D thermo-mechanical numerical calculations using the viscous-plastic I3ELVIS code 41 , which is based on a combination of a finite difference method with a marker-in-cell technique (see Methods). We use a 1500 × 1500 × 635 km rectangular model domain with 297 × 297 × 133 nodes that offers a spatial resolution of ca. 5 × 5 × 5 km per grid cell. The laterally homogenous lithosphere consists of a bi-layer, 36-km-thick crust and 114-km-thick lithospheric mantle. We initiate a mantle plume by seeding a 200 km-radius hemispheric temperature anomaly at the base of the upper mantle, 300 K warmer than the surroundings (Supplementary Fig. 1). In contrast to previous studies that use an arbitrary pre-defined triple junction geometry 4,42 , we impose no pre-existing structuration within the crust or lithospheric mantle in order to investigate the spontaneous initiation of a divergent triple junction in response to the simultaneous action of the mantle plume and far-field extensional stresses. We simulate tectonic forcing by applying a constant divergent velocity normal to the "eastern" and "western" model boundaries combined with laterally varying pull along the northern side of the model ( Supplementary Fig. 1). We perform three groups of models with EW extension half-rates of 0, 3, and 6 mm/yr (Supplementary Table 1). These boundary velocities are derived from the Neogene kinematics of the Nubia-Somalia plate system 43 . Each group consists of six experiments that use the following kinematic boundary conditions along northern side of model domain: free-slip (i.e. absence of northward pull), constant northward pull of 4 mm/yr, and northward pull linearly increased from west to east from 0-6 mm/yr to 3-6 mm/yr, 6-12 mm/yr, and 12-18 mm/yr (Supplementary Table 1). The laterally-varying northward pull is meant to mimic the geodynamics of the northern convergent margin of the Afro-Arabian plate at the time of the Afar plume impingment 18 (Fig. 1b).
The results summarized in Fig. 2 show that far-field extensional stresses play a key role in the style of system development. A non-pre-stressed lithosphere expectedly results in axisymmetric surface deformation in response to radial lateral spreading of the plume head while radial crustal cracks are not distinguishable because of the limited resolution of the model (Fig. 2b). In contrast, the presence of any far-field stress triggers the development of non-axisymmetric features that are governed by the value and orientation of the boundary velocities applied. EW linear rifts (Fig. 2c) occur in experiments where northward pull dominates EW extension. On the contrary, faster EW extension leads to continental break-up localized along extension-perpendicular (i.e. NS oriented) structures (Fig. 2d). In certain cases, two en échelon overlapping spreading zones trending in a NS direction are connected by a >200 long orthogonal transform zone, forming a complex transform-ridge pattern (Fig. 2e). An Afar-like triple junction end-member (Fig. 2f) develops when EW extension (3 mm/yr) and northward pull (3-6 mm/yr) are close. Note that full EW extension of 6 mm/yr in this model is consistent with present-day geodetic estimates for the northernmost segment of the Main Ethiopian rift 44 . We observe that the synchronous increase in both EW extension (up to 6 mm/yr) and northward pull (up to 6-12 mm/yr) does not change significantly this pattern ( Supplementary Fig. 2e). A further augmentation of the northward pull (up to 6-18 mm/yr) results in a ring-like structure with four radial axes, which are combined in a series of four interconnected triple junctions (Fig. 2g).
Four main rifting modes can thus be distinguished: (i) axisymmetric deformation (Fig. 2b), (ii) linear rift(s) stretched in EW to WNW-ESE (Fig. 2c) or NS to NW-SE (Fig. 2d) directions, (iii) ridge-transform pattern (Fig. 2e) and (iv) single triple junction (Fig. 2f) or four interconnected triple junctions resulting in a ring-like structure (Fig. 2g). Despite the theoretical stability of R-R-R triple junctions 1 , their spontaneous initiation in numerical experiments appears to be controlled by far-field tectonic forcing. In most our experiments, extensional deformation is localized along one or two linear structures oriented perpendicularly to the dominant forces, whereas more complex patterns such as divergent triple junctions occur only within a limited range of extension/pull combinations ( Fig. 2a; Supplementary Fig. 2).
We extracted shear heat production (H s ) from each successive stage of model 10 ( Supplementary Fig. 2). H s is a measure of the dissipation of mechanical energy during irreversible non-elastic (e.g., viscous) deformation ( Supplementary Fig. 3). We observe two episodes of rapid decrease of mechanical work in the system, one at 0-10 Myr when upper crustal deformation dominates the initial stage of the model evolution, the other at 40-55 Myr as the model transitions from single break-up axis to the final triple junction configuration. We interpret this as showing that a triple junction provides the optimal geometrical structure for continental break-up under the conditions of model 10 -that most resembles the Afar triple junction -as it maximizes the rate of decrease of the dissipative mechanical work required for multi-directional extension. The stability of Afar-like triple junctions thus satisfies the principle of minimum mechanical dissipative work 5 .

Implications for the Afar Triple Junction
The configuration of break-up zones in the triple junction experiment shown on Fig. 2f mimics the geometry of the Afar triple junction quite well, as shown on Fig. 3, despite the relative simplicity of the initial model setup: a NW-SE-trending northern rift branch (similar to the Red Sea rift) joins two other rift branches, one trending Thus, this model reproduces the first-order structures of the Afro-Arabian rift system from the Red Sea to the Kenyan Rift in the simple context of an initially homogenous lithosphere due to the thermal impact of the active mantle plume and its buoyancy-driven flow combined with bi-directional, asymmetric far-field stresses. The presence of pre-existing linear zones of weakness, proposed to explain the formation of the Red Sea 50,51 and Gulf of Aden 52 , appears not to be mandatory for deformation to localize and, ultimately, lead to the present configuration of the Afar triple junction.

Methods
We produced the simulations presented in this contribution using numerical code I3ELVIS that is based on a combination of a finite difference method with a marker-in-cell technique 41,53 . This parallel 3D code allows for non-diffusive numerical simulation of multiphase flow in a rectangular fully staggered Eulerian grid 54 . The momentum equations are solved in the form of Stokes flow approximation. The components of the deviatoric stress tensor are calculated using the viscous constitutive relationship between stress and strain rate for a compressible fluid. The mechanical equations are fully thermodynamically coupled with heat conservation equations accounting for mineralogical phase changes, adiabatic, radiogenic and frictional internal heating sources.
Non-newtonian viscous-plastic rheologies are implemented via evaluation of the effective viscosity of the material. The material deforms according to Newtonian diffusion creep and power law dislocation creep. In addition, Peierl's plasticity applies when stress reaches a specific limit ("Peierl's stress") that characterizes transition to plastic failure 55,56 . The contributions from dislocation, diffusion and Peierls creep are taken into account via computation of inverse average viscosity. The ductile rheology is combined with a brittle/plastic rheology to yield an effective visco-plastic rheology with a Drucker-Prager yield criterion. The visco-plastic rheology is assigned to the model by means of a "Christmas tree"-like criterion, where the rheological behaviour depends on the minimum viscosity attained between the ductile and brittle/plastic fields [57][58][59] .
The free surface topography is reproduced using the 'sticky air' technique enhanced by the introduction of a high-density marker distribution in the near-surface zone.
The implicit parallel numerical scheme uses an indirect multigrid solver 53 , that permits to accelerate 3D calculations.
In-depth description of the computer code I3ELVIS is provided in the book by T. Gerya 53 .