Mantle Convection Patterns Reveal the Mechanism of the Red Sea Rifting

We present a new model of the stress state and present‐day tectonics of the Red Sea Rift (RSR) based on an instantaneous geodynamic mantle flow model. The initial density and viscosity variations in the mantle are derived from a joint inversion of gravity, residual topography, and tomography, which provides higher resolution than existing models. The calculated mantle flow shows clear distinctions along the rift axis. The tectonics of the southern part of the Red Sea is mainly controlled by the Afar plume and characterized by divergent mantle flow. The passive rifting along the central part of the RSR can be explained either by asthenospheric upwelling due to the Red Sea floor spreading or by the plume, rising from the transition zone and not directly related to the Afar plume. We also observed ridge‐axis‐aligned flow in the uppermost mantle in the northern part of the RSR.


Introduction
The Red Sea is a narrow basin extending from the junction of the Gulf of Suez and Gulf of Aqaba in the northwest and the Gulf of Aden and Afar in the southeast. This relatively modern rift system broke up the African and Arabian plates in the late Oligocene-early Miocene (e.g., Bosworth et al., 2005;Mohriak & Leroy, 2013). It has a high level of seismicity with earthquakes magnitudes lower than 6.5 according to recent catalogues of the Egyptian National Seismic Networks (Kaban, El Khrepy, & Al-Arifi, 2018b) ( Figure 1a). The earthquakes along the axial depression of the Red Sea are clustered on oblique transform faults (e.g., Bosworth et al., 2005;Ghebreab, 1998) and on the Sinai triple junction in the northwestern edge. The Red Sea Rift (RSR) is usually divided in three zones with different tectonic regimes: (I) active rifting at the southern segment, where the seafloor spreading is developing; (II) the central segment, where incipient rifting is accompanied by abundant volcanic activity; and (III) passive extension in the northern segment, where continental lithosphere is extending with no clear signs of oceanic crust. Here, there is only minor volcanic activity (Figure 1a) during the entire history of the RSR (e.g., Almalki et al., 2015;Ghebreab, 1998;Stern & Johnson, 2010). The Arabian plate is surrounded by the Dead Sea transform fault in the west, Bitlis-Zagros suture zone in the north, Makran subduction zone in the northeast, and the Owen fracture zone with the Gulf of Aden ridge in the east (Figure 1a). These boundaries are controlled by slab pull in the Eurasia-Arabia continental collision zone, and the Afar plume below the junction of the Ethiopian, Red Sea, and Gulf of Aden ridges (e.g., Bellahsen et al., 2003;Bosworth et al., 2005;Koptev et al., 2018). Previous seismic tomography models suggest that the Afar plume is likely a continuation of the African ©2020. The Authors. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

10.1029/2019TC005829
Key Points: • The Afar plume acts only in the southern part of the Red Sea Rift forming two flow branches in the northern and nearly eastern directions • The central part of the Red Sea Rift is a passive rift not directly related to the Afar plume • The modeled stress field shows an off-axis extension with a strike-slip component along all segments of the Red Sea Rift Supporting Information: • Supporting Information S1 • Data Set S1 • Data Set S2 • Figure S1 Correspondence to: superplume in the upper mantle (Ritsema et al., 1999;Simmons et al., 2007Simmons et al., , 2009Forte et al., 2010;Hansen & Nyblade, 2013), which is in accordance with geochemical studies (e.g., Hilton et al., 2011;Pik et al., 2006).
It has been demonstrated that the global-scale forces produced by mantle flow play a key role in the regional tectonics as well (e.g., Glišović & Forte, 2016;Moucha & Forte, 2011). In previous studies, numerical modeling suggested that the mantle flows caused by the Tethyan slab subduction below the Bitlis-Zagros suture zone (Becker & Faccenna, 2011;Faccenna et al., 2013) in the north, and upwelling of the African superplume in the south (Forte et al., 2010;Moucha & Forte, 2011), play a major role in the regional tectonics of the Arabian plate. These models incorporated continental-scale mantle dynamics and used density distributions derived from a low-resolution tomography model (e.g., SMEAN, Becker & Boschi, 2002). However, a much higher resolution density distribution is required to model the 100-km scale structural features below the RSR and adjacent regions (Kaban et al., 2015).
During the past decades significant amount of geological, geophysical, and geochemical data have been collected for the RSR area. However, there remains a debate concerning timing of rifting, evolution, and current tectonic settings of different parts of the RSR (e.g., Almalki et al., 2015). A geodynamical model for the whole region, which takes into account far-field stresses induced at the plate boundaries, traction forces at the base of the lithosphere caused by mantle convection, and incorporating various independent observables, could be an appropriate tool to explain the tectonic heterogeneities observed in the three segments of the RSR. Unfortunately, existing models do not provide the resolution sufficient for geodynamic reconstruction of the RSR segments.
To investigate the driving forces responsible for the along-strike partitioning of the RSR, we developed a global mantle convection model with a focus on the Arabian plate and surroundings. We use recent highresolution density and viscosity models of the upper mantle obtained for this region in previous studies ( Kaban et al., 2016;Kaban, Petrunin, et al., 2018a). These models were based on seismic tomography (Ritsema et al., 2011;Schaeffer & Lebedev, 2013), the resolution of which has been refined in a joint inversion with the residual mantle gravity anomalies and residual topography. This method provides a significantly higher resolution as compared to the direct conversion of seismic wave anomalies from global-scale models to density and viscosity (e.g., Kaban et al., 2015). The final model of mantle flow takes into consideration both the traction forces at the base of the lithosphere induced by global mantle convection and the regional forces caused by local heterogeneity of the lithosphere and upper mantle.

Data and Methods
The method used in this study comprises two main stages. In the first stage, a high-resolution 3-D density model was developed for the study region from a joint inversion of gravity, topography, and seismic tomography data with subsequent integration into the global model (Kaban et al., 2016;Kaban, Petrunin, et al., 2018a). The second stage comprises a numerical calculation of instantaneous mantle flow velocities, stresses, and strain rates acting on the lithosphere with the use of the numerical code ProSpher 3D (Petrunin et al., 2013).
The 3-D density model is based on a joint inversion of three data sets (supporting information section S1). The first one is the residual mantle gravity field obtained after eliminating the effect of the crust from the observed one. In the same way, the residual topography, which represents the second data set, was estimated (Kaban et al., 2016). In the inversion, we have limited the model to the 700-km depth. Below, the method loses resolution. For the greater depth, we use the unmodified model based on seismic tomography. Therefore, the effect of the lower mantle (both static and dynamic) has also been eliminated from these fields based on global dynamic models (Kaban & Trubitsyn, 2012;Petrunin et al., 2013). As a result, the residual gravity and topography reflect the density heterogeneity of the upper mantle. It has been showed that using both fields in the inversion provides a possibility of reliably resolving 3-D density structure of the upper mantle (Kaban et al., 2015). The inversion is constrained by the initial density model based on the seismic tomography of Schaeffer and Lebedev (2013). For the uppermost mantle (<300 km), the initial density variations, which represent the third input parameter in the inversion, are estimated from seismic velocities using a mineral physics approach (Stixrude & Lithgow-Bertelloni, 2005;Tesauro et al., 2014). At depths greater than 300 km, a constant scaling factor is applied as done for the global model (e.g., Steinberger & Calderwood, 2006).
The final 3-D density model (supporting information Figure S1) has been obtained from the inversion with the aim of minimizing the least square differences with the initial fields (Kaban et al., 2015). The dynamic effects related to the density variations in the upper mantle have also been taken into consideration. This method has been extensively tested using numerical examples (Kaban et al., 2015, supporting information). It has been demonstrated that the input density structure of the upper mantle is successfully retrieved even when some anomaly is not present in the initial tomography model although the amplitudes of the obtained density anomalies might be reduced by damping in the last case. Thus, the final 3-D density model implies the effect of the crustal structure, topography, and mantle heterogeneity, which define the main driving forces as input for the geodynamic model. The details of the method used and discussion of the potential uncertainties can be found in the studies by Kaban et al. (2015Kaban et al. ( , 2016 and supporting information section S1. The detailed density model of the upper mantle obtained by the joint inversion of gravity, seismic, and topography data has been merged with the global density model of the lower mantle based on a direct conversion of S wave velocity perturbations according to the global seismic tomography model S40RTS (Ritsema et al., 2011) to the density variations by applying the scaling factor d (lnΔρ)/d (lnΔVs) = 0.28 (e.g., Karato, 1993;Mitrovica & Forte, 2004;Steinberger & Calderwood, 2006). Possible small-scale variations of the physical parameters below the transition zone, which are not reflected in the seismic model, should not affect principally the convection patterns in the upper mantle because of the strong viscosity contrast at the bottom of the transition zone. The viscosity variations were calculated using the homologous temperature approach (e.g., Yamazaki & Karato, 2001). Although this approach does not recognize specific types of viscosity such as diffusion, dislocation, or Peierls creep, it is sufficient for calculation of velocities and stress field without having to solve a nonlinear problem. The 3-D temperature variations were calculated from the density variations with the depth-dependent thermal-expansion coefficient (e.g., Paulson et al., 2005). This, in combination with the depth-dependent adiabatic temperature profile (Katsura et al., 2010) and radial distribution of the homologous temperature (Yamazaki & Karato, 2001), allows for a direct conversion of the density anomalies to the 3-D viscosity model (supporting information Figure S2). In addition to the viscosity model based on the seismic tomography and high-resolution density model in the study area, we implemented weak plate boundaries based on an integrated global strain rate map (Kreemer et al., 2003). At the plate boundaries, the viscosity was reduced by 4 orders of magnitude compared to the characteristic viscosity of the lithosphere. Such contrast suffices to provide efficient decoupling between tectonic plates . The global strain rate map model provides a global distribution of the strain rates derived from global positioning system measurements and geologic/geodetic observations. This gives an opportunity for introducing local viscosity changes in the plate boundaries depending on the provided strain rates. The typical width of the weak boundaries is one to two elements. Further details regarding this procedure are provided in Kaban et al. (2014) and in supporting information section S2.
To calculate instantaneous mantle flows and corresponding stresses, we used the numerical code ProSpher 3D Petrunin et al., 2013). This code allows for solving the Navier-Stokes equation coupled with the Poisson equation for the gravity field in a spectral domain and can handle high-lateralviscosity variations, compressibility, and self-gravitation at a high resolution. The resulting model includes both the upper and lower mantle and has a resolution of 120 harmonics laterally (approximately 160 km in the space domain) and 50 km in depth. Although the vertical resolution is not sufficient to accurately resolve vertical stress variations within the lithosphere, it is enough to show an integrated effect of the mantle convection on the stress state of the lithosphere. The model implies the free-slip condition for the upper and lower boundaries and suggests a solution in the no-net-rotation reference frame.

Results
The result of the calculations is an instantaneous whole mantle convection model that allows us to estimate the far-field stresses controlling the tectonic processes in the study area. In the current study, we focus on the area surrounding the Arabian plate (Figure 1a), for which the model is based on high-resolution data. The calculated plate velocities are in good agreement (both in amplitude and direction) with the observed velocities, which proves that the dynamic model is robust (Figure 1b). In addition, the vertical flow velocities in the uppermost mantle (uprising flow marked by the red and yellow regions in Figure 2a) are in good correspondence with the distribution of magmatic fields and volcanoes (black triangles, Figure 2a), which also supports the reliability of the model.

The Afar Triple Junction and Southern Segment of the RSR
Obviously, the Afar plume is the major factor that shapes the mantle flow pattern in the southern part of the RSR. The dynamic model shows that the plume head spreads from the Afar triple junction toward the southwest and is located along the Ethiopian Rift Valley (Figure 2). The resulting convection pattern of the Afar plume is typical for an active rift. The profiles across the center of the Afar Plume head (Figures 3c and 3d) reveal a complicated density structure of the upper mantle. It should be noted that the density anomalies are calculated for each depth relative to a 1-D reference model, which implies an increase in the density with the depth (Kaban et al., 2016;Kaban, Petrunin, et al., 2018a). The most significant negative density anomaly is located at shallow depths (<150 km) and probably indicates partial melting associated with the current volcanic activity (Figure 2a). The prominent negative density anomaly in the transition zone (between 410 and 660 km) might be a primary source for the ascending mantle flow owing to the alleged high water content within the layer (e.g., Katz et al., 2003). Before, the possible hydrous upwelling in the transition zone under the Afar hotspot and a possible mechanism of the water transport in the mantle were considered by Thompson et al. (2015) and Hermann and Mookherjee (2016).

10.1029/2019TC005829
The model also shows that the plume material is distributed along two main conduits in the uppermost mantle. In the western branch, the mantle material flows along the southern segment of the RSR, then northward below the Arabian platform , whereas in the main eastern branch, the mantle flow is directed to the Gulf of Aden and partially underlies the Yemen volcanic fields (Figure 2a). It is noteworthy that the path of the eastern branch does not coincide with the axis of the rift in the Gulf of Aden. Possibly, the formation of the Aden rift is controlled by large-scale tectonic processes, and the hot mantle material from the Afar plume is directed by local thinning of the lithosphere.  (Storchak et al., 2013). The upper 50-km layer was intentionally excluded from the figures since it is characterized by very strong density variations that need an extended color scale.

The Central Part of the RSR
The central segment of the RSR is usually considered as a transition zone between the southern part characterized by the ongoing active seafloor spreading and the northern segment with some signs of continental rifting (e.g., Lazar et al., 2012). These parts are separated by transverse accommodation zones allegedly inherited from preexisting weak zones (Bosworth, 1994;Jarrige et al., 1990;Makris & Henke, 1992). Remarkably, the central segment of the RSR represents the mantle flow pattern typical of active rifting since the mantle material rises in the upper mantle like a plume and spreads below the RSR segment. However, unlike the southern segment, it does not have a deep source of buoyant material, which forces the continental lithosphere to break apart (Figures 2, and 3b and 3d). The low-density zone at depths up to 150-200 km likely represents the asthenosphere material with partial decompression melting filling the gap between the moving apart plates. This supports the idea that the rifting in the central segment of the RSR is rather passive and is mainly forced by the active processes at the Afar triple junction along with far-field forces. This is also confirmed by the orientations of the maximal principal stresses, which are oblique to the ridge axis and coincide well with the transform faults direction within the central part of the RSR (Figure 4).
Recent studies have also supported the idea that the opening of the central segment of the RSR is associated with intrusion of the hot asthenospheric material into the lower part of the continental lithosphere owing to the extensional tectonic environment (Ligi et al., 2011(Ligi et al., , 2012Mooney et al., 2017). This process results in further formation of initial nuclei for oceanic crust accretion, where the ascending hot asthenospheric material forms a chain of along-axis basins with signs of the initial stages of oceanic spreading. Although it is reasonable to assume that the hot material from the Afar plume is conducted along the RSR due to much thinner lithosphere, there exists evidence for propagation of the plume material in the northern direction, oblique to the RSR axis. Recent tomography studies do not confirm negative seismic velocity anomaly below 100 km beneath the central segment of the RSR (e.g., Koulakov et al., 2016). Based on the tomographic studies, SKS splitting data, and volcanic rocks distribution, Chang et al. (2011) hypothesized the offset of the plume-related mantle flow presumably channelized by the prerifting structure (Ebinger & Sleep, 1998).

The Northern Segment of the RSR
The mantle flow in the northern part of the RSR differs significantly. This part is characterized by the nearly horizontal flow along the ridge axis at shallow depths of 50-150 km (Figures 2a and 2b, and 3d), which turns downward at greater depths in the northern edge of the RSR (Figures 2a and 2d). This change in the flow direction is likely caused by the positive density anomaly at the depths 200-400 km (Figures 3a and 3d). Two upward streams of the hot material in the upper mantle at the central RSR and Levant (Dead Sea Rift), which are associated with high volcanic activity (Figure 2a), might also contribute to the compensating descending flow below the northern segment of the RSR and Gulf of Aqaba. Unlike the southern and central segments of the RSR, the mantle flow pattern revealed by our model does not show any sign of ongoing rifting north of 24°N. However, the maximal stress orientations in the lithosphere show a normal-to-axis extension with the south-north oriented oblique component (Figure 4a) similar to the stress field in the rest of the Red Sea basin, which might be the effect of the interaction between far-field forces and the Afar plume impact.

The Stress State of the RSR and Adjacent Areas
Previous models revealed that the stress state of the Arabian plate is chiefly controlled on a regional scale by the mantle convection cell associated with subduction of the Arabian plate in the north and impingement of the Afar plume's hot material in the south (Faccenna et al., 2013;Forte et al., 2010;Moucha & Forte, 2011). Our alternative model based on the direct conversion of the tomography model to density also shows results very similar to the previous studies (supporting information Figure S3) that confirms existence of the convecting cell below the Arabian plate on a regional scale. However, this model does not recognize local mantle flow patterns in different segments of the RSR. Direct comparison of the high-resolution model and the alternative model, solely based on seismic tomography (Figures 2 and  S3), clearly demonstrates that identification of different tectonic regimes in the RSR segments requires a relatively high-resolution model.
With the use of the joint inversion technique, the resolution of the present model becomes sufficient to recover regional-to-local pattern of the mantle flow and stress fields. The color map illustrating tectonic regimes ( Figure 4a) shows a strong variation of the stress even within the RSR basin, in contrast to the low-resolution model. According to the new model, the Ethiopian rift, southern and partially central segments of the RSR are undergoing by maximum extension forces, whereas the northern part is under moderate or weak tensile stress. This also might be evidence for tectonic segmentation of the RSR.
The maximum principal stress orientations at a depth of 50 km (Figure 4a) show the normal-to-axis extension within all segments of the RSR with a remarkable strike-slip component that is likely reflected in the transform fault series in the central segment of the RSR. It is worth noting that the SKS splitting data showing the fast seismic velocity azimuth are in a good agreement with the calculated maximum principal stress directions within the Arabian plate and the surrounding areas ( Figure 4a). However, in the Ethiopian rift, the SKS splitting data are inconsistent with the model (Kendall et al., 2005(Kendall et al., , 2006. At the depth 150 km, the stress orientation becomes significantly different (Figure 4b). It can be explained by the fact that horizontal forces are generally maximal in the lithosphere, whereas in the asthenosphere, especially in the vicinity of plumes, the maximum principal stress directions might be arbitrary and mostly controlled by the asthenospheric mantle flow. Therefore, a relation of the stress orientations and the SKS splitting data is not so obvious at these depths.
This inconsistence in the narrow zone along the Ethiopian ridge might have two main reasons. First, the lithosphere is very thin in this area and the polarization of seismic waves depends on mantle flow direction, that is, on spatial variations of the strain rates in this region. The resolution of the model which is about 100 km laterally and 50 km by depth might be insufficient to reproduce this pattern. The second possible reason is that we use the maximum principal stress orientation as a proxy for the azimuth of seismic wave's polarization. For divergent boundaries, the maximum principal stress direction is to be near-perpendicular to the fault and does not reflect the hot mantle material propagation in a shallow depth just beneath the rift crest.
The flow along the rift axis is also confirmed by the model (Figure 2a). Within the RSR, at the depth 150 km, the stress orientation is changed to almost pure normal-to-axis extension and agrees well with the directions of the mantle flow at this depth ( Figure 2b).

Conclusions
In this study, we present a dynamic model of the upper mantle and use it to elucidate the origin of the tectonic processes in the Red Sea, Arabian plate, and the surrounding areas. This model is based on the highresolution 3-D density and viscosity distributions obtained in previous studies, which are embedded into a global convection model. In this way, it became possible to analyze the combined effect of the global convection and regional factors such as the Afar plume and Zagros subduction zone. As a result, the mantle flow velocities and maximum principal stress directions have been estimated for the entire area. The calculated surface velocities are in good agreement with the observed plate velocities, which supports the consistency of the model. In addition, the zones of ascending flow in the uppermost mantle correspond well with the distribution of recent volcanic activity. Furthermore, the maximum principal stress directions are in good agreement with the SKS splitting and World Stress Map project data. The main results are summarized as follows: 1. The calculated mantle flow velocities are very high in the Afar plume area, thus resulting in a strong dynamical support of the Ethiopian valley similar to the Reykjanes rift. The diverging flow originating from the Afar plume causes strong extension in a narrow area and forms two streams. The first stream flows obliquely to the Red Sea axis in the northward direction, while the second one forms the conduit below the Yemen volcanic fields and Gulf of Aden. 2. According to the obtained model, the Afar plume has no significant influence on the central segment of the Red Sea, which demonstrates the classical flow pattern of a young mid-ocean ridge with upwelling of the asthenospheric material below its axial part. 3. In contrast, the northern segment of the Red Sea Valley is characterized by an along-axis mantle flow. For this structure, we can assume an initial stage of passive continental rifting driven by the relative motion of the African and Arabian plates accommodated along the Dead Sea transform fault and Gulf of Aqaba. At the depths greater 200 km, the mantle flow descends at the northern edge of the Red Sea. 4. The maximum principal stress directions in the lithosphere point to the extension along the entire RSR.
They are generally characterized by a south-north direction, i.e., a direction oblique to the RSR axis, which implies a strike-slip component in the basin tectonics.