Use of a high-precision gravity survey to understand the formation of oceanic crust and the role of melt at the southern Red Sea rift in Afar, Ethiopia

Abstract The Red Sea arm of the triple junction in northeastern Ethiopia provides an opportunity to investigate rift-forming processes at divergent boundaries. In an attempt to study the subsurface, especially the distribution and role of melt in the rifting process, we carried out a high-precision gravity survey with a mean-square error of 0.011 mgal, assisted by differential global positioning system measurements. The profile is 162 km long and strikes ENE–WSW across the southern part of the Red Sea rift at a latitude of approximately 11.75° N. Modelling of the Bouguer anomaly, constrained by a priori information, showed detailed in-rift variations in the crustal structure and the distribution of melt beneath the rift axis. Our interpretation suggested that the process of continental break-up is governed by crustal stretching and rifting accompanied by the emplacement of melt into the lower crust above a lower density upper mantle. In addition, we interpreted the thickness of the crust beneath this part of the rift axis to be 25 km. The subsurface distribution of density beneath the profile shows that the south-central part of the Red Sea rift has modified thinned crust, intruded by high-density material, which resembles the crust formed during seafloor spreading.

Extension of the crust during continental break-up is accommodated by ductile stretching and faulting (McKenzie 1978) as well as by the intrusion of magma (e.g. Ebinger & Casey 2001). The intrusion of magma modifies the composition of the crust and is most commonly manifested globally at ancient rifted margins as high seismic velocities and high densities in the upper and lower crust, spatially coincident with crustal thinning across the continentocean transition (e.g. White et al. 2008;Ahmed et al. 2013). Despite the importance of the magma intrusion process during continental break-up, it remains unclear how molten rock is delivered into the lower crust and how it is transported through the crust to sometimes erupt at the surface. The southern Red Sea rift of Afar is the ideal study locale to address this problem as it sub-aerially exposes late-stage magmatic rifting caused by extension above anomalously warm (Ferguson et al. 2013) and seismically slow (Stork et al. 2013) mantle.
We obtained high-precision and high-resolution gravity measurements on a 162 km-long profile traversing the axis of the southern Red Sea rift in Afar. These data were used in a study of the crustal structure and to refine our understanding of the role of melt in the rifting process at divergent boundaries ( Fig. 1). To achieve this goal, the interpretation of the measured gravity data was constrained using information from boreholes, previous geophysical investigations and petrology. In addition to these data, results from spectral and Euler deconvolution methods were used as an input into a 2.5dimensional (2.5D) forward modelling process.

Tectonic setting
Afar marks a triple junction between the Nubian, Somalian and Arabian plates (Fig. 1), which are diverging as a result of extension in the Red Sea, Gulf of Aden and East African rifts (e.g. McKenzie et al. 1970;Mohr 1970;Tazieff et al. 1972;Ghebreab 1998;Beyene & Abdelsalam 2005). The excellent fit of the southern coast of Arabia into the Horn of Africa was one of the earliest case studies used to substantiate plate tectonic theory (e.g. McKenzie et al. 1970). Border faults on the SE and SW flanks of the Afar depression mark the abrupt transition from the rift valley floor to the 2-3 kmhigh Ethiopian and Southeastern plateaus. The conjugate rift flanks are located c. 350 km to the NE and define the southern tip of Arabia in Yemen (e.g. Bosworth et al. 2005). Geochronological constraints in Ethiopia suggest rifting began 29-31 Ma ago on the western Afar margin (e.g. Wolfenden et al. 2005;Ayalew et al. 2006).
Since the Oligocene, extension has migrated away from the border faults towards the current southern Red Sea rift axis in the Dabbahu-Manda Harraro rift and in the Tendaho Graben. The synrift volcanic geology of the western Afar margin shows that during the period 27 -25 Ma, extension and volcanism in the southern Red Sea rift migrated 40-50 km westwards of the main border fault. The morphology and spatial extent of lava flows suggest that intrusion and volcanism localized to narrow magmatic segments akin to those currently active in the Main Ethiopian Rift (Wolfenden et al. 2005;Ayalew et al. 2006). Over time, the locus of extension and magmatism has migrated eastwards in the Red Sea rift system (Wolfenden et al. 2005;Ayalew et al. 2006;Keir et al. 2011), leaving a suite of near-surface volcanic rocks and sedimentary basins that young towards the Quaternary-Recent axial magmatic segments, where faulting and magmatism are now localized (e.g. Barberi & Varet 1977;Manighetti et al. 2001).

Geology of the Tendaho Graben
At the southernmost extent of the Red Sea rift, the Dabbahu-Manda Harraro rift and the Tendaho Graben ( Fig. 1) define a series of NNW-SSEstriking, c. 80 km-long fault-bound basins (e.g. Manighetti et al. 2001;Acocella et al. 2008). The graben, which initiated at about 1.8 Ma (Acton et al. 1991(Acton et al. , 2000, is 30 -60 km wide (Gianelli et al. 1998) and represents the largest graben in central Afar. The basins are bounded by flanks 200-300 m higher than the lowest elevation of the basin floor, which is 300-500 m above mean sea-level. A single clear border fault is generally lacking. Instead, the graben is bordered with numerous closely spaced, relatively small-offset faults distributed across the rift margin (e.g. Manighetti et al. 2001;Rowland et al. 2007;Acocella 2010). In general, normal faulting is the most dominant structural displacement mechanism, with some sinistral strike-slip faults (Abbate et al. 1995;Ayele et al. 2015).
The major lithological unit that forms the footwalls of the graben is the Stratoid Basalts, the age of which is progressively younger towards the middle of the graben. The basins are filled with 1-2 km of lacustrine and fluvial deposits and by Quaternary-Recent volcanic rocks (Acocella et al. 2008), the distribution and age of which constrain the evolution of the basin. Faulting and volcanism have progressively localized throughout the Quaternary to a c. 10 km-wide central axis where the youngest (0.2-0 Ma) fissure basalts crop out (Lahitte et al. 2003a, b;Acocella et al. 2008;Acocella 2010). On a regional scale, the c. NNW-trending Dabbahu-Manda Harraro rift and Tendaho Graben are 50-150 km from the c. 2 km elevation margin of western Afar and c. 200 km from the c. 500 m elevation Danakil microplate on the eastern margin of Afar. The Dabbahu-Manda Harraro rift and Tendaho Graben are therefore relatively narrow basins set within a 250-450 km-wide Afar depression and have experienced deformation during the last c. 30 Ma.

Crustal structure
The structure and composition of the crust beneath Afar have been the subject of much debate during the last c. 35 years (e.g. Berckhemer et al. 1975;Makris & Ginzburg 1987;Mohr 1970Mohr , 1989Hammond et al. 2011). Discussion has often focused on whether it is completely new oceanic crust, or whether it is stretched and thinned continental crust, heavily intruded by magma (transitional crust). Constraints on crustal structure primarily come from wide-angle seismic studies (Berckhemer et al. 1975;Makris & Ginzburg 1987;Prodehl et al. 1997;Maguire et al. 2006), teleseismic receiver function studies (Ayele et al. 2004;Dugda et al. 2005;Stuart et al. 2006;Hammond et al. 2011) and gravity analyses (Redfield et al. 2003;Tessema & Antoine Fig. 1. Location of the survey profile. The profile was 35.0, 15.5, 10.5 and 10.5 km north of the towns of Mile, Logia, Semera and Serdo, respectively. All the CORS stations were within 125 km of any of the measurement points and there was at least one station within 64 km radius of the furthest points on the profile and 10 km for most of the points in the profile located in the central part. DMHR, Dabbahu-Manda Harraro rift; GNSS, Global Navigation Satellite Systems Service; MG-n, controlled source seismic profile n from the work of Makris & Ginzburg (1987), where n stands for III, IV or V; TG, Tendaho Graben; TGD, Tendaho-Goba'ad discontinuity.
The crust in Ethiopia is thickest beneath the Ethiopian Plateau (40-45 km); it is c. 35 km thick beneath the Southeastern Plateau. The crust thins abruptly to 15-30 km into Afar (e.g. Maguire et al. 2006;Stuart et al. 2006). Throughout Ethiopia the crust appears to have a consistent, laterally continuous layering consisting of lower crust with a P-wave seismic velocity of V p ¼ 6.7-7.0 km s 21 , an upper crust with V p ¼ 6.0-6.3 km s 21 and cover rocks of lava flows and sediments with V p ¼ 2.2-4.5 km s 21 (Makris & Ginzburg 1987;Prodehl & Mechie 1991). The cover rocks are generally thickest where the crust is thinnest, suggesting a strong link between crustal thinning, rift valley subsidence and the resultant accumulation of basin infill. The seismic velocity structure of the Afar crust contrasts with that observed in the submarine oceanic ridges in the Red Sea and Gulf of Aden, where the crust is less than c. 10 km thick and displays a different internal seismic architecture (Prodehl & Mechie 1991). The crust beneath Afar is considered to be variably stretched and intruded continental crust and not yet fully oceanic in nature (Makris & Ginzburg 1987).
In central Afar, the crust is 20-25 km thick beneath the current locus of volcanism and strain in the Dabbahu-Manda Harraro rift and the Tendaho Graben (Makris & Ginzburg 1987;Hammond et al. 2011). The crust is of similar thickness in a c. 40 km-wide strip near the western Afar margin (Hammond et al. 2011), where extension is inferred to have localized in the Miocene during the development of the southern Red Sea rift (Dahla basalts) (Wolfenden et al. 2005). The regions of crustal thinning in the present and past loci of strain are separated by 30 km-thick crust, similar to the situation beneath the Danakil microplate (Makris & Ginzburg 1987;Hammond et al. 2011).
In addition to crustal thickness, receiver function studies constrain the bulk crustal V p /V s ratios. In general, the plateaus have V p /V s ratios of 1.7-1.9, whereas localized regions of thicker (c. 30 kmthick) crust within Afar have V p /V s ratios of 1.8-1.9 (Dugda et al. 2005;Hammond et al. 2011). V p /V s ≤ 1.85 can be explained compositionally, but observations of V p /V s . 1.85 strongly suggest the presence of partial melt. Near the active Quaternary-Recent magmatic segment and beneath the Miocene Dahla basalts, where the crustal thickness is about 30 km, V p /V s is very high at 1.9-2.1further strong evidence for a partial melt in the crust (Dugda et al. 2005;Hammond et al. 2011).
Information about the crustal structure in Afar also comes from magnetic and gravity studies, particularly from a 50 km-long profile across the Tendaho Graben (Bridges et al. 2012). Here, the magnetic data show near-symmetrical anomalies across the centre of the rift, with the scale, shape and magnitude of the anomalies similar to those observed in nearby seafloor spreading centres in the Gulf of Aden. These data, together with the Bouguer gravity field, are consistent with a symmetrical sheeted mafic dyke complex intruded into the upper c. 10 km of the crust and capped by a 2-3 km-thick basin-fill of fissure basalts and sediments (Bridges et al. 2012).
Our work aimed to acquire further information to constrain across-rift variations in crustal thickness as well as the intra-crustal distribution and geometry of melt and solidified magma intrusions. For this purpose, high-precision gravity data were collected along a profile across the Tendaho Graben, supported by differential global positioning system (GPS) measurements to attain the highest possible accuracy in height measurement. Modelling of the variations in crustal thickness and density structure was constrained using information from boreholes, natural and controlled source seismic studies, magnetotelluric (MT) data and petrology.

Data collection
The survey profile was laid out in a straight line from the foot of the western margin of the Red Sea rift in Afar, across the axis of its southern end to the northeastern part of the town of Serdo (Fig. 1). The total length of the profile was 161.68 km. The average station interval was 0.55 km near the interpreted rift axis of the Tendaho Graben and 2-5 km at the flanks of the profile, with a single point at a distance of 7 km at the western end of the profile.
Pre-measurement work included the selection of station intervals, the location of measurement points on a map, the transfer of their coordinates to a handheld GPS instrument, and prediction of the tidal variations at five-minute intervals for the planned field session so that on-the-spot correction could be made during the measurements. In addition, the functionality of the gravimeter was checked between two benchmarks in Digdiga and Semera (Fig. 1), the absolute values of which were known to an accuracy of +0.003 mgal (Le Moigne et al. 2010).
The instrument used for the gravimetric measurements was a Lacoste & Romberg G781 gravimeter. This instrument was fully serviced and upgraded with the ALIOD100 utility in July 2007. As part of the routine check-up procedure, both the spirit and electronic levels were calibrated using a standard technique before the field sessions. During the survey, a routine gravimeter handling and check-up procedure was followed to prevent large drifts and tares in the relative gravity values.
A Garmin eTrex hand-held GPS receiver was used to transfer the planned locations of measurement points to the ground. Dual-frequency Trimble NetRS and Trimble R7 receivers with Trimblemade Zephry geodetic antennae were used for the differential GPS survey. The roaming receiver was mounted on spike bipod mounts for all measurements. The differential GPS data were referenced to the continuously operating stations DASM, DAYR, DAFT, DA35 and DA45 (Fig. 1). All the reference stations were equipped with Trimble NetRS receivers with Zephyr geodetic antennae.
The gravity profile was partially measured from October to November 2009 and extended to the east and west from October to November 2012. To tie the measurements at different epochs, measurements were repeated on four points, two at each end of the October 2009 profile. The average gravity difference from measurements taken during the different epochs at these four points was +0.0056 mgal.
From the total of 81 points measured on the profile, 40 were revisited at least twice, with a total of 213 revisits. The primary base station used for the survey was the absolute gravity station at Semera (Le Moigne et al. 2010), which is about 300 m west of the DASM continuous GPS station ( Fig. 1). In general, secondary base stations were always established for the measurements each day and the absolute station was measured only two or three times per day. The secondary base stations were mostly revisited three to four times per day. The revisits to stations were carried out systematically to include measurements from different days.
During any visit to a station, at least two readings at two-minute intervals were taken to check the repeatability of the measurements. Moreover, if a station was reoccupied for a second time, on-thespot comparisons were made between the current and previous readings, taking into account tidal variations. The maximum tolerance for discrepancies was 0.010 mgal. To guarantee the accurate reduction of temporal variations, the time stamp was taken from the Garmin eTrex hand-held GPS instrument controlled by a synchronized wrist-watch.
Together with the gravity measurements, the air pressure was observed to remove its effect on the measured gravity value. The air temperature was also measured to monitor any unusual behaviour of the gravimeter due to external temperature variations. The height of the gravimeter above the measurement surface was always registered to allow the reduction of the data to the surface of the Earth. During the two field campaigns, the gravimeter behaved very well, with an average drift of 0.0035 mgal/day and no tare was observed.
The GPS measurements at the gravity stations were made in the differential mode for an average duration of 30 minutes using dualfrequency instruments with a sampling interval of 15 seconds. As the area has little vegetation and as a result of the systematic choice of measurement locations, there were always 9-11 satellites visible at each point.

Data reduction and processing
After the field campaign, the locations of the five continuously operating reference GPS stations ( Fig. 1) were determined using GAMIT/GLOBK software (Herring et al. 2010) relative to the International Terrestrial Reference Frame 2005 by taking into account the station velocities and jumps due to tectonic events. For this purpose, a total of 13 International Global Navigation Satellite Systems Service stations, including ADIS in Addis Ababa, Ethiopia, were used. The five CORS stations were later used as reference stations to determine the locations of the gravity measurement points in the differential mode, referenced to the World Geodetic System 1984 (WGS84). The best accuracy in location was dX ¼ 0.00382 m, dY ¼ 0.00413 m and dZ ¼ 0.00534 m, while the maximum error observed was dX ¼ 0.02741 m, dY ¼ 0.02992 m and dZ ¼ 0.03722 m.
The first step in the processing of the gravity data was the removal of temporal variations from the data. This was achieved by first removing the tidal and then later the atmospheric pressure variations. The former used the Cartwright and Eden model (Cartwright & Eden 1973) and the latter the factor suggested by Merrian (1992) and Boedecker & Richter (1984). The output was then reduced to the ground surface from the instrument measuring position using a gradient term of 0.3086 mgal m 21 . Later, repeated measurements at multiple locations were adjusted using a least-squares adjustment technique by removing the gravimeter drift (Becker 1984). The resulting mean-square error of the adjustment of the data from the two epochs was 0.0113 mgal. This puts the survey in the highprecision category (Becker 1984;Rymer 1989).
The ellipsoid used to determine the normal gravity value was the Geodetic Reference System 1980 (GRS80) (Moritz 1980). Apart from being the current best global geodetic reference system for estimating the normal gravity value, GRS80 is, for any practical work, the same as WGS84 (Li & Götze 2001), which is an ellipsoid widely used by the geodetic community. WGS84 was also used as the reference ellipsoid for the location of measurement points in the GPS survey. The Somigliana closed-form formula (Somigliana 1930) was used to compute the normal gravity values.
As discussed by Moritz (1980), the mass of the Earth's atmosphere is included with the mass of the solid Earth in determining the normal gravity on GRS80. However, the mass of the atmosphere above a gravity station does not affect the station's measured gravity (Hinze et al. 2005). Therefore a correction is made for this effect using the relation by Ecker & Mittermayer (1969), as discussed by Wenzel (1985) and used by Hinze et al. (2005).
As the survey was carried out along a relatively long profile and there was a maximum height difference of about 672 m, the second-order approximation formula (Heiskanen & Moritz 1967;Hinze et al. 2005) was used for the free-air reduction. This has led to a gravity difference of about 0.026 mgal along the entire profile compared with the implementation of the first-order (0.3086 mgal m 21 ) gradient term, a difference of about 2.36 times the measurement accuracy.
For the simple Bouguer reduction, which is also known as the intermediate layer reduction, we applied the Bullard B correction, which takes into account the effect of the curvature of the Earth using the closed-form formula for a spherical cap of radius 166.7 km, as given in LaFehr (1991). This was used instead of the traditional term, which assumes an infinite horizontal slab. The use of the average crustal density (2.67 g cm 23 ) for the topographic reduction re-established the dependence of the gravity value on the topography. The density that would have been used in the simple Bouguer reduction to de-correlate the output gravity data from the topography was about 0.71 g cm 23 . As this value is considerably less than the density of the different lithological units mapped in the area and the density used in this work, it clearly indicates that there is local isostatic compensation (Hildebrand et al. 1990).
The final step we took to remove the topographic effect was the implementation of terrain reduction. To implement this reduction, it was necessary to test which topographic models best fitted the purpose. Although the accuracies of the chosen terrain models were not ideal for our purpose, the Shuttle Radar Topography Mission (SRTM) data with a 3 arc second sampling rate (SRTM3/CGIAR-CSI) (Jarvis et al. 2008) and the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER), which is a product of METI and NASA, were compared with the height from differential GPS measurements along the profile. The comparison was made by interpolating pixels of the elevation models at the location of the gravity measurement points as determined by differential GPS measurements. The result of the comparison showed that the SRTM3 data fitted better, with a standard deviation of 3.2 m, than the ASTER data, which had a standard deviation of 7.5 m, in agreement with previous findings (e.g. Kääb 2005; Zhao et al. 2011;Rexer & Hirt 2014). Based on these results, the SRTM3/CGIAR-CSI data were used for terrain reduction. The resulting complete Bouguer anomaly is presented in Figure 2, together with the raw gravity data and the topography of the profile.

Analysis and interpretation
The complete Bouguer anomaly (Fig. 2a) was dominated by a long wavelength increase in gravity from an average value of 2100 mgal in the west to 243 mgal towards the eastern part of the profile. There were also noticeable medium wavelength negative anomalies of about 29 and 212 mgal around the axis of the Tendaho Graben and 45 km west of the rift axis, respectively, that were superimposed on this long wavelength field. The negative anomaly about the axis of the rift was, in turn, superimposed by a sharp and positive anomaly of 8 mgal.

Power spectrum analysis of the gravity data
The first step in the interpretation of the gravity data was to estimate the depth to the horizontal surfaces with major density variations that were responsible for the different anomalies. For this purpose, we used the spectral analysis technique first implemented by Spector & Grant (1970) and widely used thereafter (e.g. Bhattacharyya & Leu 1975;Studinger et al. 1997;Ates & Kearey 2000;Gómez-Ortiz et al. 2011). In this technique, the power spectrum of the gravity data is usually plotted against the angular frequency and the major density variation zones are interpreted from the slopes of the linear segments.
To compute the energy spectrum we used the fast-Fourier transform technique after multiplying the Bouguer anomaly by the Hann window function. The natural logarithm of the spectrum was then computed and smoothed using a weighted moving average. Linear curves were fitted to the different segments using a piece-wise least-squares linear curve-fitting approach. As the computation of the slope of the different segments was sensitive to the addition of one or two points from one segment to the other, tests were carried out by varying the number of points forming a particular line segment. The combination of points that gave the average depth estimate were then considered as belonging to that particular segment. The standard errors were computed from the resulting depths obtained by varying the number of points forming a particular line segment.
The results (Fig. 3) indicated that a major contribution to the gravity anomaly came from vertical density variations at a depth of 24.8 + 1.4 km.
The second segment, which is also in the long wavelength spectrum but with a gentler slope, was caused by a body located at a depth of about 10.5 + 0.17 km. As the source that generated this long wavelength feature lay at a shallower depth, it was interpreted to have a broader horizontal spatial extent. On the other hand, there was another source of anomaly, the third segment, which was responsible for the signal in the medium wavelength spectrum and lay at an approximate depth of 16.7 + 1.6 km. Typically, the slopes of successive segments decrease with increasing angular frequency. However, in our case, the slope of the third segment was steeper than the slope of the second segment. There have been a few examples where similar results have been found (e.g. Ates & Kearey 2000, fig. 9b). As the interpreted depth to the source was larger than that from the second segment and the frequency was in the medium wavelength spectrum, with a relatively smaller amplitude, we interpreted the anomaly to be caused by a body that has a small horizontal spatial extent relative to the horizontal extent of the anomaly itself, with a finite depth and small density contrast. The shallowest depth at which major density variations were interpreted was 1.57 + 0.46 km. The interpretation of information from this segment needs to be made carefully, especially if the data are erroneous, for it is within the spectrum of random noise. Many studies commonly avoid the interpretation of the spectrum in the short wavelength range (e.g. Gelisli & Maden 2006) or consider it as white noise (e.g. Gómez-Ortiz et al. 2011). In our case, the uncertainty was 29% of the computed depth and the information was therefore not reliable. However, we have presented this depth for comparison with the borehole data and it will not be used as a constraint in the final model.

Euler deconvolution method
The spectral analysis technique can be used to identify depths to major vertical density variation zones, but it does not constrain the location of these zones along the profile. To estimate the location of the major horizontal density variation zones, we used the Euler deconvolution method (Thompson 1982).
As a result of the availability of information about the shallow structures from deep boreholes, and to avoid the interpretation of anomalies from the high-frequency spectrum that often causes instability during the numerical differentiation required in Euler deconvolution, the very high-frequency effects were first filtered using a low-pass filter with a cut-off frequency of 0.3 km 21 . We then computed the vertical and horizontal derivatives numerically to apply the Euler deconvolution method. We used a data window of 31 points and first identified the plateaus (Silva & Barbosa 2003) to arrive at the best estimates of the horizontal source positions. After identifying three distinct source positions at 290, 237 and 22 to 12 km, we then used the approach described by Melo et al. (2013) to identify the best structural indices in the three regions. This was achieved by testing tentative structural indices for the three source bodies separately and computing the correlation coefficient between the complete Bouguer anomaly and the base level estimates for each of the structural indices. The structural indices that gave the minimum correlation were 0.72, 0.84 and 0.48 for the three horizontal source positions at 290, 237 and 22 to 12 km, respectively. We then rounded these structural indices to integer values to take care of the cautions of Reid et al. (2014). We assigned the nearest integer value 1 for the sources at 290 and 237 km. We tested both 0 and 1 for the source between 22 and 12 km. Previous studies have shown that there is a dyke-like structure (Bridges et al. 2012) at this horizontal position, from the surface of the Earth to about 8 km depth, which is underlain by a highly electrically conductive body from 15 to 28 km depth (Johnson 2012). Therefore it is logical that the method is not able to differentiate between the dyke and the melt, which can be approximated by a dyke and a horizontal cylinder that usually give structural indices of 0 and 1, respectively (Reid et al. 2014).
We then used the method suggested by FitzGerald et al. (2004) with a conservative reliability factor of 0.91 to discriminate between reliable and spurious solutions. The choice of this reliability factor rejected the source at 290 km with a structural index of 1 and that between 22 and 12 km with a structural index of 0, but passed the source bodies at 240 km and between 22 and 12 km, both with structural indices of 1.
The final result (Fig. 4) revealed a tight assembly of points 37 km west of the interpreted rift axis at an average depth of 27.5 km. This type of assemblage for a structural index of 1 usually arises from causative bodies approximating a line, cylinder or thin-bed fault (Stavrev 1997;Reid et al. 2014). The comparison of our results with previous work (Makris & Ginzburg 1987;Hammond et al. 2011) led us to the conclusion that the causative body was similar to a thin-bed fault representing a sudden change in mantle topography. Similarly, the assemblage of points under the rift axis between 22 and 12 km was interpreted to be caused by major horizontal density variation zones between 8 and 15 km depth. As the high-frequency effects were filtered out during processing for the Euler analysis, the effects resulting from sources shallower than 3.3 km were not expected to be resolved. This might also be the cause of the poor reliability value when a structural index of 0 was used for the anomaly at the rift axis. The value of 1 for the structural index that led to the assemblage of points at a depth of 8-15 km might signify the presence of a causative body resembling a horizontal cylinder (Cooper 2008;Reid et al. 2014) with its axis perpendicular to the direction of the profile. The selection of the cylindrical source was made after comparing our results with those of Johnson (2012).

Information from boreholes and seismic profiles
There are a number of shallow (maximum depth 600 m) and deep (maximum depth 2.2 km) boreholes in the Afar region that have been drilled for groundwater and geothermal exploitation, respectively. The current gravity profile traverses very near to three (about 7.0 km to the nearest borehole) deep, geothermal exploration boreholes (Fig. 1). The lithological logs (Aquater 1996;Gianelli et al. 1998;Battistelli et al. 2002) from these boreholes show that there are 1400, 1200 and 600 m-thick sedimentary layers, with intercalations of thin basaltic layers, in boreholes TD1, TD2 and TD3, respectively. In general, these boreholes show that the thickness of sediments decreases towards the north and increases towards the centre of the Tendaho Graben.
The seismic profiles III and IV of Makris & Ginzburg (1987) are near to our profile and their profile V crosses it (Fig. 1). Their interpretation has put the boundary between the upper and lower crust at an average depth of 10 km and the Moho at a depth of 26 km. This agrees well with the major horizontal density variation zones at depths of 10.5 and 24.8 km in the spectral analysis method. The information from profile V (Makris & Ginzburg 1987) suggests that the sediments are 2.5 km thick in the south around Assaita and thin to zero about 65 km north of our profile. Their profile crosses ours approximately 53 km north of the town of Assaita (Fig. 1), where we expect the depth to the base of the sedimentary layer to be a few hundred metres.
Based on the information from the lithological logs from nearby boreholes and the published density values for such deposits (e.g. Parasnis 1952;Telford et al. 1990), a density value of 2.27 g cm 23 can be considered as an average value for the lithological units dominated by the sedimentary rocks.
As there were no direct density measurements on rock samples available, the constraints on the densities of the other lithological units were mainly made using values derived from seismic velocities obtained by Makris & Ginzburg (1987). The average velocities of the different layers, taken from their various profiles that either cross our profile or run very near to it, are given in Table 1.
Although densities have been derived from P-wave velocities using a variety of numerical relations, we only present (Table 1) density values derived using those of Brocher (2005) and Nafe & Drake (1957). Nevertheless, the density values derived from the numerical relation of Nafe & Drake (1957) agree with published average density values of the different layers (Pembroke 1969) and these were used to constrain the density values in the modelling process.

Modelling and interpretation
The complete Bouguer anomaly is dominated by a major positive jump of 40 mgal at distances between 245 and 230 km (Fig. 2a), centred on the Tendaho -Goba'ad discontinuity (Fig. 1) at about 237 km. A gravity low centred at about 250 km has noticeably affected the smooth transition from the low Bouguer anomaly values in the west to the high anomaly values in the eastern part of the profile. On the other hand, a low amplitude broad negative anomaly asymmetrically centred at 5 km and another relatively higher frequency positive anomaly symmetrically centred at the interpreted rift axis are superimposed on the gravity high in the eastern part of the profile. There are also very highfrequency anomalies along the profile that have a low amplitude and are not noticeable unless the dominant low-frequency and high-amplitude effects are filtered out.
The long wavelength anomaly with a major jump between 245 and 230 km is interpreted to be caused by an abrupt change in the topography of the Moho (Fig. 5). This is especially well captured in the Euler deconvolution results and is in good Table 1. Densities of the different formations derived from P-wave velocities (Makris & Ginzburg 1987) using the numerical relation reported by Nafe & Drake (1957) and Brocher (2005) Formation P-wave velocity (m s 21 ) Density (g cm 23 ) Profiles used as sources of P-wave velocities Brocher (2005) Nafe & Drake ( The density of the upper crust was derived from the average velocity of what was taken as the volcanic and upper crust velocity in Makris & Ginzburg (1987). agreement with the recent receiver function study by Hammond et al. (2011) and the controlled source seismic experiment of Makris & Ginzburg (1987). In our study, the shallowest depth to the Moho from the surface of the Earth is about 25 km beneath the interpreted rift axis, unlike the 20 km depth inferred from the gravity data of Makris & Ginzburg (1987, fig. 11). However, our Moho depth estimate is in good agreement with that from their seismic profile at the place where it crosses our profile. However, our 3.2 g cm 23 density value resulting from the modelling process is lower than the upper mantle density, which is 3.25 -3.38 g cm 23 (Pembroke 1969). Our result maps the topography of the Moho with a significantly improved resolution compared with previous studies and demonstrates that the crust in this area is thinned and stretched. Although the location of the relatively broad, lower amplitude, negative anomaly around the central part of the profile coincides with an area of sedimentary cover, especially alluvial and lacustrine deposits (see geological map in Aquater 1996; Gianelli et al. 1998), it was not possible to ascribe the full magnitude of the anomaly solely to the sediments (Fig. 5). This would have been possible only if the sediments had a thickness of about 3.4 km or a density of 1.5 g cm 23 . Such large thicknesses are unlikely because the seismic work of Makris & Ginzburg (1987) and the lithological logs from boreholes TD1, TD2 and TD3 have shown that the thickness of the sediments is not more than 1.2 km. These studies also showed that the aggregate thickness of the sediments decreases northwards and should only be a few hundred metres at the location of our profile, before disappearing further to the north. Moreover, the average density of the sediments, as derived from the seismic profiles of Makris & Ginzburg (1987) and using the relation of Nafe & Drake (1957), is about 2.29 g cm 23 (Table 1), whereas that based on the lithological logs from the boreholes is about 2.27 g cm 23 . In addition, the width of the anomaly covers about 60 km and the surface exposure of the sediments is not more than 30 km. Sediments with thicknesses ,1.2 km and a width of 30 km cannot cause such an anomaly. To explain the relatively broad, lower amplitude, negative anomaly around the central part of the profile, we suggest the presence of a deeper source. Euler deconvolution showed an assemblage of points from a depth of 8-15 km and the spectral analysis method found an interface in the intermediate spectrum window. This requires the presence of a body with a lower density than the surrounding lower crustal country rock to partially model the negative anomaly. Given the geological context (Ferguson et al. 2013), the most viable interpretation of this low-density material is a basaltic melt in a magma reservoir.
To test this assumption, we considered the physical properties of a magma reservoir. We initially assumed 100% basaltic melt of the same composition as that supplying recent dyke episodes and axial fissure eruptions and calculated the density using the model of Ochs & Lange (1999). Following Ferguson et al. (2013), we adopted the composition of the axial basalt erupted in 2009 with a ferric to total iron ratio of 0.15, a pre-eruptive H 2 O content of 0.5 wt% and a pre-eruption storage temperature of 12008C. We calculated a density of 2.748 g cm 23 at a pressure of 0.5 GPa and an average depth of 15 km. A similar density (2.754 g cm 23 ) was calculated for the 2007 axial basalt lava under the same conditions. Changing the storage temperature by +508C changes the density by only 0.010 g cm 23 and changing the storage pressure by +0.1 GPa (3 km) changes the density by 0.016 g cm 23 . As greater storage depths would probably result in higher temperatures, the two effects tend to cancel each other out. The presence of crystals within the reservoir would increase the magma density. In terms of our modelling, there may be some trade-off between the crystal content of the magma reservoir and its dimensions. Optimization during modelling of the gravity data yielded an excellent fit with a density of 2.824 g cm 23 (Fig. 5). This is higher than the calculated density of pure melt, suggesting a zone of partial melt rather than 100% melt. The partial melt may be distributed through a higher density crystalline framework in the form of a partially molten mush, or concentrated into vertical (or horizontal) ribbons separated by solid rock in the form of a dyke (or sill) complex. We cannot distinguish these two possibilities from gravity data alone because they both have the same average density. However, if we take the melt density to be 2.748 g cm 23 and that of the solid lower crust to be 2.959 g cm 23 , the proportion of melt required to give a bulk density of 2.824 g cm 23 is 65%, which is more consistent with a plexus of closely spaced dykes or sills than a crystal-rich magmatic mush.
Our estimated melt fraction is larger than that calculated by Johnson (2012) based on the MT profile across the Tendaho Graben, which lies close to our gravity profile. The highest electrical conductivity material occupies a region approximately 13 km wide, at depths of 15-28 km beneath the rift axis. This body is very similar in size and depth to the body (2.824 g cm 23 ) used in our gravity modelling (13 km width, depth range 8.5-25 km). On the basis of the high conductivity of the body (0.11-0.30 S m 21 ), Johnson (2012) concluded that it represented a region of partial melt. Using similar melt compositions, pressures and temperatures to those given here, and a variety of mixing rules for mixtures of partial melt and solid matrix, Johnson (2012) estimated melt fractions in the range 5-8%, appreciably lower than the values we used to model the gravity profile. There are several possible explanations for this discrepancy. First, the calculated physical properties (conductivity and density) are subject to uncertainty. The melt conductivity (Pommier & Le-Trong 2011) is currently the least well-known of the calculated parameters with a relative uncertainty of 0.4 log-units (factor of 2.5); the melt density (Ochs & Lange 1999) is known to within 0.03 g cm 23 , which is 1% of the estimated density. Consequently, the two melt fraction estimates may converge if due attention is paid to uncertainties in the model parameters. For example, the observed conductivity of the partially molten region lies between 0.11 and 0.30 S m 21 , while the range of calculated melt conductivities, accounting for the uncertainty in melt chemistry and the temperature, is 1.7-3.5 S m 21 , which expands to 0.7-14 S m 21 once the uncertainty in the SIGMELTS (Pommier & Le-Trong 2011) algorithm is considered. For the highest conductivity (0.3 S m 21 ) of the body and the lowest conductivity of the constituent partial melt (0.7 S m 21 ), the estimated melt fraction increases to 64% at the upper Hashin-Shtrikman bound (Johnson 2012). This melt fraction corresponds to a bulk density of 2.824 g cm 23 , which is the same as that obtained from the gravity modelling. On the other hand, using higher densities in the gravity modelling to arrive at a lower melt fraction requires a significant increase in body size, which is inconsistent with the MT result and other constraints used to model the gravity data. Thus, a value of 65% melt represents a satisfactory result given the assumptions and uncertainties associated with modelling data from the gravity and MT profiles.
The modelled magma reservoir has a geometry that is narrow at a depth of about 8.5 km and broadens at depths between 15 and 25 km (Fig. 5). The modelling of its shape between 8.5 and 15 km depth is guided by the Euler deconvolution result. The broader interface between the magma reservoir and the surrounding rocks at 15.3 km is in good agreement with the geodetic estimates of the source depth of dykes intruded during the current magmatic episode in Afar (Hamling 2010;Keir et al. 2013). The location of the magma reservoir is also consistent with the MT result (Johnson 2012). The results of other MT surveys in the area (Desissa et al. 2013;Didana et al. 2014) also showed the presence of a highly electrically conductive body that is best explained by the presence of partial melt in the area. The maximum thickness of the sediments underneath this negative anomaly is about 0.75 km. From this result we can conclude that the process of crustal stretching and thinning under the central part of the Tendaho Graben is accompanied by magma intrusion.
The peculiar short wavelength positive anomaly centred at the rift axis is interpreted as caused by a dyke complex extending from the underlying magma reservoir at 8.5 km depth to the base of the sedimentary units at 0.75 km (Fig. 5). Its horizontal location agrees with the work of Bridges et al. (2012). Their modelling, based mainly on a magnetic survey, was carried out 20 km south of our profile and was only interpreted to a depth of 8 km. For this body, we initially considered a bulk density equivalent to the density of the lower crust, which was later modified to a density of 3.0 g cm 23 to fit the gravity data using this structure. Note that the density needs to be significantly higher than that of basaltic melt. Thus, we interpret this body to be a solidified dyke (or dykes) supplied by the magma reservoir at depth, but which has crystallized as a result of more rapid cooling at shallower depths. The surface manifestation of the Ayrobera thermal springs (Aquater 1996; Gianelli et al. 1998), just above this 1.45 km-wide highdensity material, might signify the presence of a zone of high heat flow underneath. As a result of its high density, which is well above the average density of upper crustal material, and the possibility that it is within a higher temperature zone (Aquater 1996), we believe that this unit has the bulk composition of basaltic oceanic crust.
Another significant negative anomaly was found 250 km from the interpreted rift axis. We interpreted this gravity low using the sediments mapped at the surface. The composition of the sediments at this location is similar to that within the Tendaho Graben. The thickness of this unit is 0.95 km. However, because there is no independent information from boreholes or any other geophysical technique, further investigation is needed to ascertain the current interpretation. The anomaly lies west of the Tendaho -Goba'ad discontinuity, which marks the active and ancient boundary between the East African and the Red Sea rifts (Ebinger et al. 2008), and the possible existence of deeper low-density material cannot be ruled out as an alternative model -thus our interpretation in this particular area should be considered cautiously. We recommend further investigation to resolve this ambiguity.
The general positive linear trend from the west towards the east is partially modelled by the rise in the Moho towards the east (from a depth of 33 km in the west to 25 km about the interpreted rift axis) and partially by the linear and smooth depth variation between the lower and upper crust (Fig. 5). In general, the data derived from the 2.5D model agrees very well with the measured data, except at the ends of the profile. The misfit at the ends of the profile also results from missing structures beyond the ends of the profile that influence the gravity data near the end-points.

Conclusion
Modelling of new high-precision and high-resolution gravity data measured on a 162 km-long profile across the southern end of the Red Sea rift in Afar has revealed detailed in-rift variations in the crustal structure, the configuration of the magma plumbing system, and the distribution of ductile crustal stretching and thinning beneath the axis of the Tendaho Graben. Along the profile, the Moho gently rises from a depth of 33 km at the western margin of the Red Sea rift to a depth of 25 km underneath the axis of the Tendaho Graben, with a sharp gradient about 40 km west of the rift axis. Our model suggests that the process of continental break-up is dominated by magma intrusion occurring where crustal thinning is at a maximum. Magma intrusion is characterized by the emplacement of basaltic melt in the lower crust, either as dykes or sills, over an upper mantle that has a lower than average density. The magma reservoir in the lower crust is connected to the surface by dykes that supply the volcanism present in the axial region. The dykes solidified quickly as a result of heat loss close to the surface, whereas the deeper, lower crustal magma reservoir appears to have remained partially molten for longer periods.
In general, the results of the analysis of the gravity data showed that the upper crust beneath the axis of the Tendaho Graben is intruded by a high-density material similar to oceanic crust and is fed from below by a magma plumbing system. This implies that the crust in Afar is stretched and thinned continental crust, interspersed by high-density, solidified basalt, which might indicate that the Afar region is at an early stage of the style of extension that characterizes modern ocean basins.
The fieldwork for this research was funded by a Consortium Grant (NE/E005284/1) from the UK Natural Environment Research Council (NERC). We acknowledge the multi-fold support from Addis Ababa University, the Ethiopian Federal Government and the Afar Regional State Government. JB thanks the European Research Council (Advanced Grant CRITMAG) and the Royal Society Wolfson Research Merit Award for additional support. DK is supported by NERC grant NE/L013932/1. We are also grateful to the University of Edinburgh and the British Geological Survey for making available the modelling software used in this research. The support of Ethioder Tour and Travel Company is also appreciated. We also acknowledge the use of borehole information from the Geological Survey of Ethiopia. The discussion with Catherine Annen was informative. Esubalew Adem is thanked for assistance during the second field campaign.