In this study we image crustal structure beneath a magmatic continental rift to understand the interplay between crustal stretching and magmatism during the late stages of continental rifting: the Main Ethiopian Rift (MER). The northern sector of this region marks the transition from continental rifting in the East African Rift to incipient seafloor spreading in the southern Red Sea and western Gulf of Aden. Our local tomographic inversion exploits 172 broad-band instruments covering an area of 250 × 350 km of the rift and adjacent plateaux. The instruments recorded a total of 2139 local earthquakes over a 16-month period. Several synthetic tests show that resolution is good between 12 and 25 km depth (below sea level), but some horizontal velocity smearing is evident along the axis of the Main Ethiopian Rift below 16 km. We present a 3-D P-wave velocity model of the mid-crust and present the first 3-D Vp/Vs model of the region. Our models show high P-wave velocities (6.5 km s−1) beneath the axis of the rift at a depth of 12–25 km. The presence of high Vp/Vs ratios (1.81–1.84) at the same depth range suggest that they are cooled mafic intrusions. The high Vp/Vs values, along with other geophysical evidence, suggest that dyking is pervasive beneath the axis of the rift from the mid-crustal depths to the surface and that some portion of partial melt may exist at lower crustal depths. Although the crustal stretching factor across the Main Ethiopian Rift is ∼1.7, our results indicate that magma intrusion in narrow zones accommodates a large proportion of extensional strain, with similarities to slow-spreading mid-ocean ridge processes.
Current models of continental breakup differ in the mechanism of strain localization in time and space, but they often ignore the effects of magma intrusion (e.g. Ruppel et al. 1995; Huismans & Beaumont 2003). Numerical models of rifting that include magmatism (e.g. Buck 2004, 2006) show that dyking into the lithosphere reduces the yield stress by an order of magnitude and thus significantly aids continental breakup. Indeed, a large proportion of passive margins show evidence of magmatism: seaward-dipping reflector sequences, indicative of syn-rift extrusive lavas, and high-velocity lower crust interpreted as igneous underplating. However, the role of magmatism in the breakup process remains poorly understood (e.g. Holbrook & Kelemen 1993; Menzies et al. 2002).
The northern sector of the Main Ethiopian Rift (MER) marks the transition from continental extension in the East African Rift to incipient seafloor spreading in the southern Red Sea and western Gulf of Aden (Fig. 1). The Ethiopian Afar Geoscientific Lithospheric Experiment (EAGLE) project imaged crust and upper mantle structure beneath this area using controlled source and passive seismic, gravity and magnetotelluric experiments (Maguire et al. 2003). This work provides a strong framework in which to study crustal modification by stretching and magma intrusion in 3-D (e.g. Keranen et al. 2004; Bastow et al. 2005; Cornwell et al. 2006; Maguire et al. 2006; Keir et al. 2006a; Whaler & Hautot 2006). EAGLE built on the previous 2-D seismic refraction studies of Berckhemer et al. (1975) re-interpreted by Mechie & Prodehl (1998). Wide-angle seismic profiles across and along the MER (Mackenzie et al. 2005 and Maguire et al. 2006) constrain the crustal and uppermost mantle velocity structure. Receiver function studies map variations in crustal thickness between the rift and adjacent plateaus (Dugda et al. 2005 and Stuart et al. 2006). Keranen et al. (2004) used controlled source tomography to image crustal velocity variations. Their study based only on P-wave data, had good resolution in the uppermost 8–15 km in the central rift zone centred on the Boset magmatic segment (Fig. 1).
Here we develop a 3-D image of crustal structure in the transitional MER using local earthquakes recorded during the EAGLE passive experiment. We present a 3-D P-wave velocity model and the first 3-D Vp/Vs model for the mid-crust (12–25 km) beneath the MER. We compare our velocity variations and the new Vp/Vs patterns with those of wide-angle and teleseismic receiver functions studies over the same region. The aim of this study was to examine the crustal structure beneath the MER and the results have important implications for the role of magma intrusion during the transition from continental to oceanic rifting.
The NE–SW trending MER lies on the Ethiopian plateau that is thought to have developed above a mantle plume (e.g. Schilling 1973; Ebinger & Sleep 1998). Within the Afar depression, the MER intersects the Gulf of Aden and Red Sea Rifts to form a triple junction between the Nubian, Arabian and Somalian plates (Fig. 1). The plateau is covered by 2-km-thick flood basalts. They are thickest along the Red Sea margins in Ethiopia and Yemen, where they have been dated at 31–29 Ma (e.g. Hofmann et al. 1997; Baker et al. 1996). They were erupted roughly coeval with the onset of rifting in the southern Red Sea by 28 Ma (Wolfenden et al. 2005). Rifting within the southern and central MER initiated between 18 and 15 Ma (WoldeGabriel et al. 1991; Ebinger & Casey 2001) and at ∼11 Ma in the northern MER (Wolfenden et al. 2004).
The broad uplifted plateau regions of NE Africa and SW Arabia are underlain by anomalously hot asthenosphere, but debate continues on the morphology of the low velocity zones (e.g. Benoit et al. 2003, 2006; Montelli et al. 2004). Bastow et al. (2005) imaged a narrow (∼75 km wide) low velocity zone (δVp = 1.5 per cent, δVs = −4 per cent) that extends to <75 km beneath the MER; the relative magnitude of the delays indicate the presence of partial melt. Kendall et al. (2005) note the correlation of the shear wave splitting parameters with zones of Quaternary magmatism, as evidence for melt-filled cracks through the lithosphere. Similar patterns in crustal shear wave splitting suggest that anisotropy in the uppermost crust is controlled by aligned dykes and partial melt (Keir et al. 2005). Lithospheric xenoliths hosted in Quaternary rift lavas support the widespread dyking interpretation (Rooney et al. 2005).
Structural and geodetic data indicate that extensional strain in the MER has migrated from the Miocene border faults to ∼60 km long, <20 km wide en-echelon ‘magmatic segments’ located within the central rift valley (Boccaletti et al. 1998; Bilham et al. 1999; Ebinger & Casey 2001) (Figs 1 and 2a). The magmatic segments are the locus of magmatism and active faulting; they are marked by N–S striking fissures, aligned eruptive centres, and small offset faults (e.g. Wolfenden et al. 2004; Casey et al. 2006). Controlled source seismic studies show that the magmatic segments are underlain by high velocity (Vp > 6.5 km s−1) elongate bodies that are interpreted as cooled mafic intrusions (Keranen et al. 2004; Mackenzie et al. 2005). These magmatic segments are also characterized by relative positive Bouguer anomalies modelled as mafic intrusions (Mahatsente et al. 1999; Tiberi et al. 2005; Cornwell et al. 2006).
GPS data indicate that ∼80 per cent of extensional strain across the Nubian–Somalian plate boundary is accommodated beneath the magmatic segments (Bilham et al. 1999; Bendick et al. 2006). Keranen et al. (2004) propose that, below a depth of ∼7 km, extension is controlled by magmatic intrusion in a ductile middle to lower crust. Seismicity is concentrated near the top or above these dyke injection zones, indicating that stresses induced by the magma intrusion control faulting in the brittle crust above the dykes (Keir et al. 2006a).
Seismic and gravity data indicate a northward increase in extension and magmatic modification of the crust beneath the MER (Tiberi et al. 2005; Maguire et al. 2006; Stuart et al. 2006). Crustal thickness beneath the northern MER decreases from 38 km in the south beneath the caldera lakes (7–8.5° N) to 24 km beneath Fentale Volcano in the southern Afar depression (Dugda et al. 2005; Maguire et al. 2006). The crust beyond the rift zone is 35 km thick and 45–50 km on the southeastern and northwestern sides of the rift, respectively (Dugda et al. 2005; Mackenzie et al. 2005; Stuart et al. 2006). Vp/Vs ratios, estimated from receiver function analyses, increase both from the rift flanks to the Quaternary magmatic segments, and northward along the length of the MER (Dugda et al. 2005; Stuart et al. 2006). Thus, seismic and gravity data suggest both crustal underplate and mafic intrusion zones, but minor amounts of crustal stretching south of Fentale Volcano (Dugda et al. 2005; Mackenzie et al. 2005; Tiberi et al. 2005; Stuart et al. 2006) (Fig. 2a).
Local seismicity data from the EAGLE project
The EAGLE network was the largest seismic experiment yet deployed on the African continent (Maguire et al. 2003). Local earthquakes were recorded on 29 broad-band seismic stations that were operational between late 2001 October and 2003 January. The number of seismic stations was increased during the final 4 months of the experiment with the deployment of a further 50 broad-band instruments. These stations were deployed at ∼15 km spacing mainly in the rift valley, and were operational between 2002 October and 2003 February. Data from a further 93 instruments at 5 km spacing along the cross-rift controlled source profile, that recorded for the period between 2002 November and 2003 January, were also used in order to supplement data from the high density rift valley experiment, but note their deployment pattern was not designed for high resolution tomography. Data from the 17 controlled sources (shots) were used to assess hypocentral location accuracy.
A total data set of 13 388 P wave and 12 725 S-wave arrivals were picked from 2139 local earthquakes recorded at a minimum of 4 of the 172 stations (Fig. 2b). The rms of the data set was 0.6 s after the initial location estimation with EAGLE rift model (Fig. 5). The P-wave arrivals were picked on the vertical component and S-wave arrivals picked on the horizontal components. The quality of the P- and S-wave arrivals was assigned during initial data processing (Keir et al. 2006a). Arrival quality factors range from 0 to 4, where 0 is the highest quality pick and 4 is the lowest. 78 per cent of arrivals had a quality factor of 0 or 1. Local magnitudes where calculated with most events having a magnitude of less than ML ≈ 3 (Fig. 2b) (Keir et al. 2006b). The estimated average picking uncertainty for our P-wave arrivals was 0.125 and 0.165 s for the S-wave arrivals. Note that all depths mentioned and plotted in this paper are referenced to below sea level.
1-D Crustal Model for the Mer
Seismic tomography is a non-linear inverse problem and the results depend on the initial reference model. A common way to obtain a suitable initial model is to first establish an a priori 1-D model that accounts for the known geology of a region. For this purpose traveltime data are jointly inverted with hypocentre locations and station delays using Velest (Kissling 1995). The resultant ‘minimum 1-D model’ (Kissling 1988; Kissling et al. 1994) forms the starting point of the subsequent 3-D inversion. The minimum 1-D model is designed to locate the events with the smallest possible uniform location error and represents the velocity model that most closely reflects a priori velocity information. This model also has the advantage of being useful in future standard earthquake locations, which are still based on 1-D models.
Only good quality events may be used for the calculation of the minimum 1-D model and 3-D inversion due to the intrinsic coupling between hypocentre locations and seismic velocities. Therefore, only events with at least 8 P arrivals and recorded within an epicentral distance of twice the focal depth along with the greatest angle without observation as seen from the epicentre (referred to as GAP) of < 180°, were selected from the 2139 events mentioned above. This resulted in a subset of 283 events for the inversion of a minimum 1-D P-wave velocity model.
Following the recommendations of Kissling et al. (1994) we start with a large number of thin layers (2 km thick), and we combined the layers for which velocities converged to similar values during the inversion process. The initial velocities for each layer were chosen by examining previously published work (Mechie & Prodehl 1998), and new velocity information provided by the EAGLE cross-rift and along-axis wide-angle profiles (Mackenzie et al. 2005; Maguire et al. 2006). The inversion was stopped when the earthquake locations, station delays and velocity values did not vary significantly in subsequent inversions. The final minimum 1-D P-wave velocity model (Fig. 3b and Table 1), which includes station delays, produced a final rms error of 0.128 s, which is close to the a priori picking uncertainty of the data.
The accuracy of this minimum 1-D P-wave velocity model was assessed using several tests. First, the inversion was started using velocities in each layer significantly higher and lower than the minimum 1-D model (Fig. 3a). It can be seen that above 8 km and below 25 km there are no earthquakes, resulting in poor convergence of the initial models. Therefore, our models will only accurately image mid-crustal depths. The next test used to assess the suitability of the minimum 1-D P-wave velocity model was to randomly adjust hypocentre locations by 8–10 km in latitude, longitude and depth prior to input into the inversion. During the inversion the velocities are fixed. The aim of this test is to see if the minimum 1-D model is able to relocate the hypocentres in their original positions (Fig. 4). It is seen that most events are relocated close to their original positions and the largest errors occur in depth. Since 17 explosive shots were recorded by the network, they were located using the minimum 1-D model as a final test of the model. The errors on the relocated shots averaged 0.51 km and 0.70 km in longitude and latitude, with an average depth error of 6 km below sea level, which equates to ∼7 km subsurface. The large error in depth is not surprising since our model has little control on velocities above 8 km where the wide-angle results show considerable lateral variations (Maguire et al. 2006). The rms error of the observed versus predicted arrival times for the 283 events relocated with the min 1-D P-wave model was 0.20 s.
After the 1-D inversion of the P-wave velocities, S-wave arrival times were included to jointly invert for a P- and S-wave 1-D model. Therefore, all events were relocated with the minimum 1-D model and a Vp/Vs ratio of 1.75 estimated from a Wadati diagram of first arrival data from events selected for the inversion. Similar data selection criteria as the P-wave inversion were used but a minimum of 6 S arrivals and 8 P arrivals (215 events). In general focal depth errors obtained with a fixed Vp/Vs ratio are nearly twice as large as errors using P phases only (Maurer & Krandolfer 1996). Therefore, a minimum S-wave model is determined using an additional series of inversions. Four Vp/Vs ratios (1.60, 1.70, 1.80 and 1.90) were used to construct four initial S-wave velocity models (Fig. 3b). In each case the previous minimum 1-D P-wave velocity model incorporating station delays was taken as the initial P model. The poorly resolved upper 8 km of the P-wave model was held fixed in the joint P- and S-wave inversion. Therefore, for depths shallower than 8 km the Vp/Vs ratio is close to the initial value. Likewise the Vp/Vs ratios diverge at depths greater than 30 km as there are no recorded earthquakes at this depth. Fig. 3(b) shows the final P, S and Vp/Vs models for the MER. The accuracy of the resulting P- and S-wave model was tested using the same three steps as the P-wave velocity model discussed above (e.g. Table 2).
The minimum 1-D P model for the MER (Figs 3b and 5) shows a strong increase in velocity below 8 km depth. The wide-angle reflection/refraction profiles (Mackenzie et al. 2005; Maguire et al. 2006) show this increase as a sharp step at a much shallower depth of 0 km (sea level) beneath the southeastern plateau, 2 km beneath the northeastern plateau and 4 km beneath the rift (Fig. 5). Mackenzie et al. (2005) interpreted this step as the boundary between lower velocity (<5.0 km s−1) sedimentary and volcanic rocks and the 6.1–6.2 km s−1 cyrstalline crust. This sharp contrast is approximated by a gradual velocity gradient in the minimum 1-D model because there are few earthquakes at these depths. Upper crustal velocities increase from 6.09 km s−1 at 8 km to 6.35 km s−1 at 25 km depth, and the lower crust then has a velocity ∼6.70 km s−1 (Figs 3b and 5). As with the P-wave velocity model, the S-wave model shows a velocity increase at 8 km and then velocities increase uniformly with depth. The Vp/Vs model has an average value of 1.75, shows a negative gradient at ∼8 km and a gradual increase between 10 and 25 km.
All 2139 events relocated with the min 1-D P- and S-wave model are shown in Fig. 2(b). Note the clustering of events beneath the Fentale Dofan magmatic segment, along the Ankober border fault, and south of the Aluto Gedamsa magmatic segment (Keir et al. 2006a). The depth range of the resultant seismicity is approximately 8–14 km beneath the magmatic segments with some isolated deeper (∼20 km) events beneath the northwestern plateau (Keir et al. 2006a).
3-D Inversion of Local Earthquake Data
For the 3-D inversion all events were relocated with the minimum 1-D P- and S-wave models. Only events with at least six P-wave observations and six S-wave observations along with a GAP less than 180° and within an epicentral distance of twice the focal depth were selected from the 2139 events. The random shift test was used to check for stable earthquake locations (Table 2). The final data set for the 3-D inversion consists of 346 earthquakes with 4970 P-wave arrivals and 3832 S-wave arrivals (Fig. 6).
The inversion code simulps14 (Thurber 1983; Eberhart-Phillips 1990), extended by Haslinger & Kissling (2001) was used for full 3-D ray tracing to calculate traveltime residuals and simultaneous inversion for hypocentre locations and 3-D velocity structure (Vp and Vp/Vs). The method inverts S–P arrival times for Vp/Vs by projecting S–P arrival-time residuals into Vp/Vs variations (Thurber & Atre 1993).
The code solves the non-linear, coupled hypocentre-velocity problem by a linearized, iterative, damped least-squares method (Thurber 1983; Eberhart-Phillips 1990). Each iteration consists of an inversion for P-wave velocity and Vp/Vs variations and for hypocentre locations. New ray paths and new traveltime residuals are computed using the updated velocity models to start the next iteration. The complexity and resolution of the resultant velocity model will depend on the earthquake station spacing, model parametrization and parameter damping (Kissling et al. 2001).
With the uneven distribution of stations and distinct clustering of earthquakes (Fig. 6), we will have some artificial velocity smearing in our final velocity model. We therefore, analyse model parametrizations and resolution. A suitable model parametrization accounts for a priori knowledge of Earth structure and the resolution capability of the available data, with trade-offs between parametrization and resolution (e.g. Kissling et al. 2001). In this study, a graded inversion scheme is used to produce smoother velocity variations and fewer artefacts (e.g. Kissling et al. 1994; Eberhart-Philips & Reyners 1997; DeShon et al. 2006). A coarse grid (60-km-nodal spacing) inversion for P-wave velocity is performed, and the resultant velocity model is input into a finer inversion scheme with a 30 km nodal spacing. Extensive testing of alternate inversion pathways, that is, single step inversions with numerous nodal spacings and different coarse and fine grid nodal spacings indicated that the chosen pathway yielded the smoothest velocity model with the fewest artefacts. Grid nodes must be sampled by at least 10 rays to be included in the inversion and inversion, and velocities at grid nodes outside the station coverage were fixed.
Vertical grid spacing was taken from the minimum 1-D model with nodes at −4, 0, 8, 12, 16, 20 and 25 km below sea level (Table 1). This coarse vertical nodal spacing is a function of the number of earthquakes used and the station spacing coupled with the need to maintain an overdetermined inversion problem. The Velest minimum 1-D P-wave velocity model with its constant velocities in each layer was then converted (Fig. 3b) into a 3-D velocity model with linear interpolation between grid nodes needed for the 3-D inversion code. Our model consisted of 672 unknown P-wave parameters, 650 Vp/Vs parameters and 1384 hypocentre parameters. With 4970 P-wave observations and 3832 S-wave observations we get a ratio between known and unknown of around ∼2.4 for the P-wave inversion, and ∼3.3 for the Vp/Vs inversion, that is, the problem is over constrained.
Damping is introduced to handle the underdetermined nature of the inverse problem. Appropriate damping parameters are determined by examining the trade off curves between data variance and model variance (Eberhart-Phillips 1990). From a series of one iteration inversion tests, a damping of 30 is used for P-wave velocity and 50 for Vp/Vs. Station delays are also damped to allow data residuals to be incorporated into the velocity model, rather than the correction terms, during the first inversion step. Tests for appropriate station damping yielded a preferred value of 8. Note that care must be taken with the interpretation of absolute velocity anomalies when damping is used.
Once the final P-wave velocity model was obtained the S waves were included in the inversion for Vp/Vs. An initial Vp/Vs ratio of 1.75 was used as it represented an average of the minimum 1-D Vp/Vs velocity model. During the inversion the P-wave velocity model was heavily damped to ensure it remained effectively fixed (e.g. Husen et al. 2000).
A plot of ray coverage (Fig. 6) shows a high number of crossing rays for P and S pairs in the rift centred on the Fentale Dofen magmatic segment. The ray coverage deteriorates south of the EAGLE cross-rift wide-angle profile, where seismicity rates are much lower. Throughout the study region there are relatively few sources outside of the rift zone, and thus the resultant velocity will have limited sampling of the ‘background’ crust outside of the rift zone compared with the crust within the rift zone. Also note that the ray paths are strongly anisotropic, with many more ray paths along the rift than across it. This will cause smearing along the rift.
Resolution may be assessed using the hit count, the resolution diagonal element (RDE), the derivative weight sum (DWS), and spread function (e.g. Haslinger et al. 1999; Reyners et al. 1999; Husen et al. 2000). The hit count is the number of rays hitting each node. The RDE shows the amount of independence in the solution of one model parameter. The DWS is also used, as it quantifies the relative ray density in the volume of influence of a node, weighting the importance of each ray segment by its distance from a model node. It is more reliable than a hit count alone. Finally, the spread function, which is based on the relative values of the diagonal versus rows of the full resolution matrix, is used to constrain velocity smearing at each node (Toomey & Foulger 1989; Michelini & McEvilly 1991). We consider regions where the spread function is <2, the DWS is >50 and the RDE is >0.2 are the minimum requirements for classification of areas of moderate to good resolution.
Resolution tests using synthetic data are also useful in defining areas of good and poor resolution, along with testing of the amplitude and geometry recovery of the data set. In this study, chequerboard tests (e.g. Zelt 1998) were used to assess the overall resolution of the data. Computation of theoretical traveltimes through the checkerboard model was performed using the 3-D ray shooting method of Simulp14 with the same source–receiver distribution and inversion parameters as our experiment. Fig. 7 shows the test result with ±8 per cent variations in P-wave velocity and Vp/Vs ratio extending over 2 × 2 grid nodes in latitude and longitude and within each vertical layer of the velocity model. Lateral resolution is good in the centre of the rift and on the northwestern plateau between 12 and 25 km depth. Within the 25 km deep layer, resolution is good but NE–SW velocity smearing becomes apparent in both P-wave velocity and Vp/Vs models. The separation between right-stepping en echelon magmatic segments (∼10 km) is less than our velocity node spacing. Therefore, we will not be able to differentiate between a continuous or segmented high velocity zone beneath the en-echelon magmatic segments along the rift axis (e.g. Figs 1 and 2). Recovered amplitudes greater than 80 per cent are common in layers 12–20 km, but in the 25 km layer amplitudes of P-wave velocities anomalies and Vp/Vs ratio are ∼50–75 per cent. In the 12 km layer Vp/Vs ratios on the northwestern plateau are >100 per cent at a few nodes. This may be the result of vertical smearing from below with most rays travelling upward from the rift (Fantale Dofen magmatic segment) to the plateau. Consequently, care must be taken in interpreting Vp/Vs ratios on the plateau at this depth and above.
Leveque et al. (1993) showed that the ability of the checkerboard to resolve fine-scale structure does not necessarily imply good resolution of larger scale anomalies but it is useful in assessing the amount of image blurring in the data. Following Husen et al. (2000) and DeShon et al. (2006), we constructed a characteristic model that reflects the broad scale geometry of our 3-D solution. We put a low P-wave velocity and Vp/Vs anomaly with a magnitude of −8 per cent beneath the axis of the rift, and a positive anomaly of +8 per cent beneath the northwestern plateau (Fig. 8). We assigned opposite signs to the expected anomalies to mathematically decouple the synthetic model from the real data (Haslinger et al. 1999). We compute theoretical traveltimes using the 3-D ray shooting method of Simulp14 and the same source–receiver distribution and inversion parameters as our data.
Recovery of the input structure is good for the P-wave velocity model but with some horizontal smearing to the east above 9° N, particularly in the 12 km layer (Fig. 8). The sharp boundary between the northwestern plateau and the rift is well resolved at all depths. The majority of earthquakes occur between 12 and 20 km depth where the recovered P-wave velocity anomaly is 80–90 per cent of the prescribed anomaly. Below 20 km depth the recovery decreases to 50–70 per cent, but the general pattern of the prescribed anomalies remain. This implies that the magnitude of the anomalies are underestimates of the true velocity structure and this must be taken into account in any interpretation. Recovery of the Vp/Vs input structure is also good, but the output model contains positive anomalies along the eastern boundary of the rift and southeastern plateau where we have little or no ray coverage. Therefore, any such anomalies in our inversion results are likely to be artefacts.
3-D P Velocity and Vp/Vs Model of the Mer
The final P-wave tomographic model of the northern MER converged after three iterations with a data variance reduction of 58 per cent and a final rms of 0.092 s, which is on the order of the a priori picking uncertainty of 0.125 s. The P-wave velocity inversion was then used as input for the Vp/Vs inversion, for which the solution converged after five iterations with a 66 per cent reduction in data variance and a final rms of 0.159 s, of the order of the a priori P-wave picking uncertainty of 0.165 s. The mean absolute value changes in event locations (relative to the 1-D model locations), following the 3-D inversion of the combined data set of P and S–P traveltimes are, 0.23, 0.64 and 1.0 km, in latitude, longitude and depth, respectively. The algorithm provides formal 1σ hypocentre uncertainties, but these may significantly underestimate errors in real data sets with non-uniform ray distribution (e.g. Husen 1999). The formal 2σ uncertainties, 0.2 s for origin time and 0.25, 0.22 and 0.83 km for latitude, longitude and depth, respectively match the location shift in moving from the 1-D to 3-D model. The accuracy of the final hypocentre locations was tested by locating the 17 EAGLE explosive shots using the 3-D velocity model and randomly adjusting the hypocentres between 8 and 10 km in latitude, longitude and depth prior to locating them with the 3-D model. From these tests, the accuracy of the 346 earthquake locations is ∼0.6 km in latitude and longitude and ∼1.5–2.0 km in depth.
The final P-wave velocity and Vp/Vs models are shown in Figs 9–11. Areas where nodes are hit by less than 50 rays are masked in grey. The DWS of each layer is shown beside the P-wave velocity and Vp/Vs models, with final earthquake locations superimposed (Fig. 9). A vertical cross-section along the EAGLE cross-rift, wide-angle reflection/refraction profile shows absolute P-wave velocities and Vp/Vs ratios (Fig. 10). Fig. 11 shows an enlarged version of the 16 and 20 km depth slices through the 3-D model with the geological features of Fig. 2(a) superimposed and known artefacts outlined with ellipses.
At depths of 12–25 km, high P-wave velocity (6.3–6.5 km s−1) and high Vp/Vs ratio (1.80–1.84) anomalies are observed beneath the axis of the rift. Note the increase in Vp/Vs ratios from 1.72 to 1.76 beneath the northwestern plateau to 1.79–1.84 within the rift. In addition, a high P wave anomaly (∼6.5 km s−1) and high Vp/Vs ratio (∼1.83) is observed beneath the Debre Zeit volcanic lineament along the western rift margin at 8.6°N. This is particularly evident in the 16 km depth slice (Fig. 9), although the resolution is inadequate to determine its continuity with the anomalies beneath the rift axis. There is also a low P-wave velocity anomaly (12–20 km) at 39.8E, close to the Ankober border fault (Figs 9 and 11) and a sharp P-wave velocity and Vp/Vs contrast at depths of 16–25 km trends N–S from the southern tip of the Ankober border fault. Based on our synthetic tests (Figs 7 and 8) these anomalies are well resolved, whereas anomalies beneath the southeastern plateau and southeastern plateau margin are probably artefacts.
The average P-wave velocity of 6–6.4 km s−1 (from the min 1-D P-wave velocity model) in a depth range of 10–20 km is consistent with the global average values of Christensen & Mooney (1995) and Christensen (1996). Our results are also consistent with values derived for the Pan-African (Late Proterozoic) crust beneath the East African Rift in Kenya (e.g. Maguire et al. 1994), and derived from the independent EAGLE controlled source studies (Mackenzie et al. 2005 and Maguire et al. 2006). To assess bulk crustal composition, we compare our 1-D velocity models (which are average values over the entire study area) to predicted rock velocities from reference compositions (e.g. Christensen & Mooney 1995, Brocher 2005a; Figs 12b and c). The average P-wave velocity of ∼4.8 km s−1 at depths above 8 km are consistent with the results from the refraction studies (Mackenzie et al. 2005) and are interpreted as a mixture of Miocene–Recent sediments and volcanics. Brocher (2005b) shows that P-wave velocities of Cenozoic sedimentary rocks and volcanic rocks in Southern California have P-wave velocities less than ∼5 km s−1 at a depth range of 0–6 km (Fig. 12). At mid-crustal depth (where we have good model resolution) there is no obvious indication from the 1-D models that the mid-crust beneath the MER and adjacent plateaus is more mafic than global norms (Fig. 12a).
The Vp/Vs ratios at a depth shallower than 8 km (Table 1) would appear to be too high for the corresponding P-wave velocity, and therefore, S-wave velocities would appear to be slightly low (by ∼0.3 km s−1Fig. 12a), if derived using the empirical relationship between P-wave velocity and S-wave velocity of Brocher (2005a), but the trend is similar and the discrepancy is on the order of the errors (∼±0.25 km s−1) in our velocity measurements in this depth range (Table 1). At depths greater that 8 km the S-wave velocities are similar.
Along axis crustal structure
In our 3-D model we observe high P-wave velocities and Vp/Vs ratios at mid-crustal depths beneath the axis of the MER as delineated by the Aluto Gedamsa, Boset Kone and Fentale Dofan and Angelele magmatic segments (Figs 9 and 11). In a controlled source tomography experiment, Keranen et al. (2004) image an along-axis segmentation in the upper crust with high velocity (6.5 km s−1) zones beneath Boset Kone and the southern end of the Fentale Dofan magmatic segments. Mackenzie et al. (2005) also image a high P-wave velocity zone beneath the central Boset Kone magmatic segment.
The coincidence of high velocity zones in the upper and mid-crust, positive gravity anomalies and the aligned Quaternary eruptive centres within the magmatic segments supports an interpretation of the high velocity, high density bodies as frozen magma intrusions. Forward modelling of gravity profile coincident with the wide-angle refraction profile and comparisons with the results of Christensen & Mooney (1995) led Cornwell et al. (2006) to suggest that these intrusions contain at least 40 per cent gabbro.
Vp/Vs ratios provide additional constraints on rock composition. Below the melting point of many rocks, mineralogy is the important influence on Poisson's ratio and Vp/Vs (e.g. Christensen 1996). The abundance of quartz and plagioclase feldspar is the dominant effect in igneous rocks. For example, felsic quartz-rich rocks such as granite have a Poisson's ratio of 0.24 (Vp/Vs = 1.71); intermediate rocks have a Poisson's ratio of 0.27 (Vp/Vs = 1.78); and mafic plagioclase-rich rocks, such as gabbro, have a Poisson's ratio of 0.30 (Vp/Vs = 1.87). Brocher (2005b) also shows that melt is not needed to explain Vp/Vs ratios between 1.81 and 1.87; these values are consistent with solid mafic rocks. The observed increase in Vp/Vs ratios (as a result of the 3-D inversion) from beneath the margin of the northwestern plateau into the rift implies the mid-crust is more mafic beneath the rift. However, with our poorer resolution away from the rift axis, we cannot confidently interpret structure within the interior of the northwestern plateau.
In addition to the observed change in composition from a more felsic plateau to a more mafic crust beneath the rift axis, high Vp/Vs ratios and high P-wave velocity zones along the rift axis could indicate the presence of a molten fraction in fractures within solidified mafic intrusions. Independent observations from the MER suggest the presence of partial melt within the lower crust beneath the rift. Whaler & Hautot (2006) use magnetotelluric methods to image a low resistivity (0.8 Ωm) zone at ∼1 km depth beneath Boset volcano, located at the southern end of the Boset Kone magmatic segment, (Figs 1 and 2). Magnetotelluric results also show a zone of conductive material with a resistivity of ∼0.5 Ωm at depths of 15–20 km beneath Boset Volcano (Fig. 11). Whaler & Hautot (2006) interpret this conductive body as a cooled mafic intrusion that contains a small melt fraction at depth. Small amounts of interconnected partial melt have been shown to have a large effect on electrical resistivity with only a few percent partial melt needed to decrease resistivity significantly (Roberts & Tyburczy 1999). In addition, the marked increase in crustal seismic anisotropy at depth of 0–20 km beneath the rift axis compared to adjacent plateau and parallelism between anisotropy and strike of fissures and aligned eruptive centres, has been used as evidence that dyke intrusions in magmatic segments penetrate to near surface levels (Keir et al. 2005). Therefore, all seismic and MT results suggest that dyking is pervasive beneath the axis of the rift from mid-crustal depths to the surface and that some portion of partial melt may exist at lower crustal depths. Resolution is too poor to evaluate shallow magma reservoirs beneath individual volcanoes shown on Fig. 2(a).
The pattern of increased Vp/Vs ratios from the plateau into the rift is also observed in receiver function studies, but our values are lower than those of Stuart et al. (2006) and Dugda et al. (2005). Stuart et al. (2006) observe Vp/Vs of 1.85 beneath the northwestern plateau and values increasing to >2 beneath the Debre Zeit volcanic chain and the Fantale Dofen magmatic segment (Fig. 2a). The most likely reason for the difference in Vp/Vs ratios from the receiver function studies and the local tomography results presented here is that the receiver function studies (e.g. Dugda et al. 2005; Stuart et al. 2006) sample the entire crust; we, in contrast sample only the mid-crust. Lower crustal reflectors identified by wide-angle studies are interpreted as the top and bottom of a 10–15 km thick zone of high velocity (7.38 km s−1) material interpreted as mafic underplate beneath the northwestern plateau and parts of the rift (Mackenzie et al. 2005). High Vp/Vs ratios continue into the upper mantle; Bastow et al. (2005) used teleseismic traveltime tomography to image low velocity zones (δVp = −1.5 per cent, δVs = −4 per cent) in the upper mantle beneath the MER. Thus, our results are consistent with the presence of lower crustal underplating and some partial melt in the lower crust beneath the rift.
The most likely reason for the difference in Vp/Vs ratios from the receiver function studies and the local tomography results presented here is that receiver function analysis provides a depth averaged Vp/Vs for the entire crust, whereas we image the mid-crust. The results of Stuart et al. (2006) and Dugda et al. (2005), therefore, sample a greater depth range than this tomography study. This could indicate the presence of more mafic material and/or melt in the lower crust and/or upper mantle. Underplate has been identified beneath the northwestern plateau and parts of the rift from cross-rift controlled source studies that image 10–15-km-thick high velocity layer at the crust–mantle boundary (Mackenzie et al. 2005; Maguire et al. 2006). Thus, our results indirectly support lower crustal underplating and some partial melt in the lower crust beneath the axis of the rift.
Off axis crustal structure
High P-wave velocities and Vp/Vs ratios are imaged beneath the Debre Zeit chain of Quaternary eruptive centres. This region marks the intersection of the MER border faults and the reactivated Precambrian Ambo fault to the west of the main rift axis. Xenoliths derived from crustal depths exhibit evidence of melt beneath the Debre Zeit volcanic chain (Rooney et al. 2005). Bastow et al. (2005) used teleseismic traveltime tomography to image low velocity zones in the upper mantle (<100 km depth) which are laterally offset to the side of the rift with the highest flank topography. Kendall et al. (2005) and Bastow et al. (2005) suggest that a steep lithospheric–asthenospheric boundary beneath the rift margins enhances melt extraction, but the melt ponds beneath the rift axis where lithostatic stress is lowest. However the spatial correlation of melt in the upper mantle, crustal underplating beneath the cross-rift wide-angle refraction line (Mackenzie et al. 2005), and high P-wave velocities and Vp/Vs ratios beneath Debre Zeit suggest that a small portion of the melt intrudes the crust directly above the low velocity mantle anomaly. Geochemical data also suggests that the source of the melt creating our high velocity and Vp/Vs anomaly is at a greater depth than that of the magmatic segments along the rift axis. However the Debre Zeit region lacks the closely spaced, morphologically young fault scarps, and active seismicity of the rift axis (Wolfenden et al. 2005; Casey et al. 2006; Keir et al. 2006a). It also lacks the large relative positive Bouguer anomalies observed over the magmatic segments (e.g. Mahatsente et al. 1999; Mickus et al. 2004). Therefore, the region may represent either a failed or underdeveloped (i.e. incipient) zone of strain.
Sharp velocity gradients are visible at mid-crustal depths at the western boundary of the rift (Figs 9 and 11). The sharp boundary of this anomaly correlates with a steep gravity gradient and rapid decrease in crustal thickness into the rift (Tiberi et al. 2005; Stuart et al. 2006). The contrast, therefore, delineates the western boundary of the rift.
Comparison with non-volcanic rifts and subaerial oceanic rifts
Crustal structure within the northern MER contrasts with crustal structure in the Woodlark Rift in Papua-New Guinea, a highly evolved backarc rift without extrusive volcanism prior to breakup (e.g. Goodliffe et al. 1997; Taylor et al. 1999). P-wave velocities from a tomographic inversion of microearthquakes in the western Woodlark basin (Ferris et al. 2006) reveal a sharp lateral seismic velocity contrast in the crust between a continental region ahead of the ocean ridge spreading tip (P-wave velocities of 5.6–6.9 km s−1 at a depth range of 10–20 km) and an oceanic region (P-wave velocities of 6.5–7.2 km s−1 between 8 and 18 km depth), close to the modern zone of seafloor spreading. Ferris et al. (2006) believe that this boundary represents a first order change in crustal composition, from a felsic/intermediate crust in the region of active continental rifting, to a more mafic crust formed by seafloor spreading. Further east and closer to the active spreading tip velocities increase from 6.5 to 7.2 km s−1 between 8 and 18 km depth. These faster crustal velocities indicate a mafic crust similar to oceanic crust. Their results also indicate that crustal thinning has localized to a narrow zone marked by high velocities in the lower crust indicative of intrusive mafic material. The crustal stretching is much more pronounced in the Woodlark basin (crustal thickness ∼16–18 km) than in the MER (crustal thickness ∼25–40), but similar patterns are observed: extensive mafic addition to the crust within a narrow zone of localized strain. The MER magmatic segments, therefore, suggest that transitional crust can be produced at relatively small stretching factors where basaltic melt is available.
The crustal structure beneath the MER can also be compared to seafloor spreading zones in Iceland. The Northern Neovolcanic Zone consists of seven rift segments arranged en echelon along the Mid Atlantic ridge plate boundary. Each volcanic system consists of a central volcano, where productivity is highest, and a transecting fissure swarm for example Krafla, Askja and Bardarbunga (e.g. Brandsdottir et al. 1997). Crustal thickness is ∼19 km under the northern Neovolcanic Zone and 35 km outside the neovolcanic zone (Staples et al. 1997). Brandsdottir et al. (1997) image a 40-km-wide high velocity zone (∼6.7 km s−1) extending from the lower crust to a depth of ∼5 km beneath the Krafla rift segment and a width of ∼20 km. High Vp/Vs ratios (1.75–1.79) beneath the active rift zone indicate a more intermediate to mafic crustal composition (e.g. Staples et al. 1997; Menke et al. 1998; Allen et al. 2002). The along-axis segmentation, crustal velocity and Vp/Vs ratio patterns in the MER are comparable to patterns in the Iceland Neovolcanic zone. Thus, magma intrusion into a zone of strain localization in the MER may mark the transition from continental rifting to seafloor spreading without the reduction in crustal thickness seen in the Woodlark rift.
Our results, interpreted in light of independent observations and numerical models, indicate that a substantial portion of the mid-crust beneath the axis of the rift, the current loci of strain, is cooled mafic material. Our observations support the magma-assisted rifting model of Buck (2004, 2006) in which breakup is facilitated by rift axial injection of magma into the lithosphere and accompanying strength reduction via localized heating. Thus progressive rise of the mafic material during successive rifting episodes produces transitional crust.
This study presents a 3-D P-wave velocity model of the mid-crust beneath the MER and presents the first 3-D Vp/Vs model of the region. High P-wave velocities (6.5 km s−1) are observed beneath the axis of the rift at a depth of 12–25 km and the presence of high Vp/Vs ratios (1.81–1.84) in the same depth range, along with other geophysical evidence, suggests that they are cooled mafic intrusions. An additional high P-wave velocity and Vp/Vs ratio zone is imaged beneath the Debre Zeit volcanic chain located 30 km west of the current locus of rifting. This chain of Quaternary volcanoes may represent either a failed or underdeveloped (i.e. incipient) zone of strain. The high P-wave velocities and Vp/Vs ratios beneath the axis of the rift along with magnetotelluric and crustal seismic anisotropy results suggest that dyking is pervasive beneath the axis of the rift from mid-crustal depths to the surface.
Substantive helpful comments by Tom Brocher and an anonymous reviewer whose have improved greatly the quality of this manuscript. We acknowledge Laike Asfaw and technical staff in the Geophysical Observatory Addis Ababa University. Thanks also to Kathy Whaler, Peter Maguire, Christel Tiberi, Tyrone Rooney, Simon Klemperer, Katie Keranen and the EAGLE working group. Special thanks to Alex Brisbourne, SEIS-UK and Mark Longbottom for technical and computer support. Figures were created using GMT (Wessel & Smith 1998) and TOMO2GMT, VELEST and SIMULP codes were made available through Stephen Husen, who provided invaluable assistance at the early stages of this project. This work was supported by NERC grant NER/A/S/2000/01004 and NERC Studentship NER/S/A/2002/10547 and a NUI, Galway Millennium Research Fund grant.