The Main Ethiopian Rift (MER) is an active continental rift where magmatic intrusion is thought to play a key role by accommodating extension and thermally weakening the lithosphere (Kendall et al., 2005; Daniels et al., 2014). Since the start of the Ethiopia-Afar Geoscientific Lithospheric Experiment ( Bastow et al., 2011) in the early 2000s, consecutive and successful controlled and passive seismic deployments helped to delineate the seismic structure of the MER crust, especially the P-wave velocity (Vp) structure and crustal thickness (e.g., Ebinger et al., 2017). A major finding of previous P-wave images is that the Vp of the crust beneath the MER is faster than that of standard continental crust (Zandt and Ammon, 1995), a feature interpreted as caused by post-Miocene mafic intrusions that have accommodated extension (e.g., Keranen et al., 2004; Mackenzie et al., 2005).

More recently, however, the advent of ambient noise tomography at periods sufficiently short has facilitated imaging of shear-wave velocity (Vs) of the crust (e.g., Kim et al., 2012; Chambers et al., 2019). Results from these studies show that the MER crust has far slower Vs than the standard continental crust, with the absolute magnitude of the velocities in places interpreted to require the presence of partial melt (e.g., Chambers et al., 2019). The joint crustal seismic properties of relatively fast Vp and slow Vs are peculiar and poorly explored in previous literature. In addition, all previous constraints on the Vs structure of the MER come from models derived from surface wave imaging techniques, such as ambient noise tomography, with a lack of independent constraints provided by alternative methods.

In order to address this and provide additional and independent constraints on the Vs structure of the MER crust, we applied the receiver function (RF) techniques using open-source codes from Computer Programs for Seismology (CPS) (Herrmann and Ammon, 2004) to estimate the velocity of the crust and upper mantle. To this effect, we have used 17 new seismic stations deployed as part of the 2016–2017 RiftVolc project (Greenfield et al., 2019; Lavayssière et al., 2019) to improve our understanding of the spatial variations of the crustal Vs structure within the central MER (CMER) and adjacent Eastern Plateau. In addition, we use the RF technique to constrain the Vp and Vp/Vs ratio. In investigating the heterogeneous structure, we have chosen two vertical cross-sections to represent the area of our study (Figure 1). One profile (A–A′) is along the rift, and the other profile is across the rift (B–B′). This study improves on the previous velocity models and Moho depth estimates of the CMER and Eastern Plateau (Dugda et al., 2005; Keranen et al., 2009) by using a relatively large number of broadband seismic stations compared with the previous studies.

The CMER formed within the Oligocene Ethiopian flood basalt province and is thought to have initiated at between ∼20 and ∼10 Ma (e.g., Wolfenden et al., 2004; Bonini et al., 2005; Corti, 2009). The extension was initially localized to several ∼60-km-long, NE-SW striking, high-angle (>60°) border faults that bound the rift, such as the Boru-Toru and the Goba-Bonga structural lineament on the western side of the rift, and the Asela-Sire Border Fault on the eastern side of the rift (Bonini et al., 2005; Corti et al., 2020). Since the Quaternary, the locus of tectonic and magmatic activity within the CMER is thought to have become focused to a ∼20-km-wide zone of small offset faults, aligned cones, and active volcanic centres within the rift valley floor known as the Wonji Fault Belt (WFB), and also at a few rift marginal magmatic systems, such as the Silti-Debre Zeyit Fault Zone (SDFZ) towards the western side of the rift (Woldegabriel et al., 1990; Rooney et al., 2014; Chiasera et al., 2018).

Constraints on the crustal structure in the CMER come from several geophysical techniques including seismology, magnetotellurics (MT), and inversion of gravity data. Constraints on crustal thickness come from sparse RF measurements (Dugda et al., 2005; Keranen et al., 2009; Kibret et al., 2019), the wide-angle controlled-source along-rift EAGLE project line, and the intra-crustal Vs structure using ambient noise tomography (Kim et al., 2012; Chambers et al., 2019).

Previous RF studies in the CMER using relatively sparse station spacing show that the crust is 27- to 40-km thick in the rift (Dugda et al., 2005; Keranen et al., 2009) and 33- to 44-km thick beneath the plateaus (Dugda et al., 2005; Kibret et al., 2019). The crustal structure beneath the Western Plateau is more heterogeneous (33- to 44-km thick) (Keranen et al., 2009) than that beneath the Eastern Plateau (38- to 41 km thick) (Dugda et al., 2005; Keranen et al., 2009; Kibret et al., 2019). Wide-angle controlled-source seismology provides further constraints in the CMER, and Western and Eastern Plateaus. The EAGLE across-rift line shows similar crustal thicknesses of 38–40 km beneath the CMER (Dugda et al., 2005; Maguire et al., 2006; Stuart et al., 2006), and 35–45 and 37–42 km beneath the Western and Eastern Plateaus, respectively (Dugda et al., 2005; Stuart et al., 2006; Cornwell et al., 2010; Hammond et al., 2011; Kibret et al., 2019).

The southern end of the EAGLE along-rift wide-angle controlled-source line is as far south as Lake Hawasa (near HAWA station) and shows a varied along-rift crustal structure in the CMER (Mackenzie et al., 2005; Maguire et al., 2006). The Vp structure, modelled by the wide-angle 2D profile studies, shows that the velocities of the upper crustal layers beneath the rift are 5%–10% higher than outside the rift, a feature interpreted to be caused by mafic intrusions associated with magmatic centres (Mackenzie et al., 2005; Maguire et al., 2006). Consistent with this, 3D controlled-source tomography of the upper crust by Keranen et al. (2004) imaged rift parallel high Vp (∼6.5–6.8 km/s) elongated bodies with a size of 20-km wide and 50-km long, and interpreted them as cooled mafic intrusions that are separated laterally from one another in a right-stepping en-echelon pattern, which corresponds with the surface segmentation of the WFB. These fast Vp regions correlate to a region of distinct positive Bouguer anomalies in gravity studies that are modelled as regions of dense rock (∼3,000 kg/m³) such a gabbro (e.g., Mahatsente et al., 1999; Cornwell et al., 2006).

In addition to the earliest studies revealing crustal structure in the MER based on Vp structure, later studies applied surface waves to render the Vs structure. Ambient noise tomography has been used to construct Rayleigh-wave group velocity maps covering the northern MER (NMER), CMER, and southern MER (SMER), and parts of the surrounding plateaus (Kim et al., 2012; Chambers et al., 2019). Chambers et al. (2019) also presented an absolute 3D Vs model of the crust and uppermost mantle of the region. An important feature of the Vs images is that the MER crust is mostly significantly slower than away from the rift, in contrast to the Vp, which is generally faster within the rift. The absolute Vs of less than 3.20 km/s + 0.03 in the lower crust are difficult to explain except with the presence of a fluid phase in the rock, such as partial melt (Chambers et al., 2019). In addition, slow Vs (<3.6 km/s) in the uppermost crust observed by Chambers et al. (2019) is consistent with the presence of sediments and/or partial melt (Diaferia and Cammarano, 2017).

Some studies reported that the anomalous high temperature is an important player on velocity structure in the case when it can trigger the transition of α–β quartz. In case of hydrated compositions (as one can presume about the current case study for the rift zone), the amphibole breakdown at increasing pressure and temperature produces a discontinuity that can be detected by RF or refraction studies (Guerri et al., 2015; Diaferia and Cammarano, 2017).

Similarly, several MT studies carried out in the CMER identify high conductivity anomalies associated with young surface volcanism (Whaler and Hautot, 2006). These conductive anomalies tend to be imaged in the uppermost crust at <1–2 km, in the upper crust at ∼3–6 km depth, and in the mid-lower crust at 20–25 km depth (Ebinger et al., 2017; Hübert et al., 2018). The shallowest anomaly is interpreted as being caused by hydrothermal fluids, whereas the other deeper high conductivity anomalies are interpreted to be caused by partial melt in the subvolcanic plumbing system (Ebinger et al., 2017; Hübert et al., 2018). Broadly speaking, there is a good correlation between the loci of slow Vs from seismology and high conductivities, giving additional remark to the idea that these anomalies are caused by partial melt (Chambers et al., 2019).

We applied an RF technique using time series teleseismic earthquakes to provide constraints on the local velocity structure of the crustal and upper mantle (Langston, 1979; Ammon et al., 1990). To extract the RF for each event, we first window the three-component seismograms starting from 10 s before and 40 s after the predicted P arrival. Selected teleseismic seismograms are rotated to radial (R), tangential (T), and vertical (Z) components from east–west, north–south, and vertical components, respectively. Each pair of horizontal-component signals (i.e., north–south and east–west components) was rotated to their corresponding radial and transverse directions.

A straightforward frequency domain deconvolution can be unstable due to spectral holes in the vertical component, and stabilization of this process can be obtained by either “pre-whitening” (Roninson, 1982; Yilmaz, 2001); or “water-level” algorithms. The former adds a small component of random noise to the vertical component, while the latter sets a lower bound on the magnitude of the denominator terms (the vertical seismogram spectral elements) in a frequency domain spectral division. In this study, converted phases are isolated by iterative, time-domain spiking deconvolution (Gurrola et al., 1995; Ligorria and Ammon, 1999) with pre-whitening to stabilize the filtering. Iterative time domain deconvolution works well even with complex signals. However, regardless of a deconvolution algorithm, the response at the receiver depends on the complexity of structures. Simple structures generally lead to better RF images (Ligorria and Ammon, 1999). After deconvolving the vertical from the radial component, we removed the signature of source, travel path, and instrumental response effects (Langston, 1979; Ammon et al., 1990; Dugda et al., 2005; Kibret et al., 2019) employing the signals coming from four different back azimuths (Figure 3).

The RF technique is a time series when the radial component trace is deconvolved from its vertical component seismogram, where the timing and amplitude of the RF phases are sensitive to the near receiver local Earth structure beneath the seismic station (Langston, 1979). The dominant signal in the first few seconds of the RF is the Ps conversion from the Moho and/or intracrustal velocity contrast followed by reverberated phases within the crust (e.g., Last et al., 1997; Hammond, 2014). In case of using a relatively dense array, RFs can show the fine crustal heterogeneity, anisotropy, and dipping structures (Eckhardt and Rabbel, 2011; Liu and Niu, 2012; Niu and James, 2002; Thybo et al., 2019; Youssof et al., 2013, Youssof et al., 2015).

Each RF was deconvolved for 20 iterations with a limiting error of 0.001 by applying three different Gaussian width parameters of 0.5, 1.0, and 2.5. We applied an iterative deconvolution algorithm (Kiknchi and Kanamori, 1982), which is calculated by the division of the denominator from the numerator (Herrmann and Ammon, 2004). Also, in each case, we allowed iteration to continue until the change in misfit resulting from the addition of a spike was 0.01% (Ligorria and Ammon, 1999). The degree of fit between the synthetic and observed RFs is calculated from the three Gaussian width parameters. A sample of two RFs is selected from 17 stations based on their percent of fit to demonstrate the overall results throughout each step (Figures 4, 5).

The study applies the ak135 velocity model (Kennett et al., 1995) as the initial velocity model to calculate the best fit velocity structure. Finally, we identified the level of the model fit of the observed and synthetic models by using both visual inspection and the calculated percentage of signal power fit. When the synthetic signals show a high degree of a misfit from the calculated RFs, both RFs and the synthetic models are automatically discarded.

The observed (red colour) and synthetic (blue colour) RFs (Figure 5) as well as the initial and the final velocity models (Figure 6) are calculated by using programs from Herrmann and Ammon (2004). The final velocity models are calculated from the global velocity model ak135. The calculated absolute velocity values at every 2-km depth are obtained from the inversions of the RFs. The uncertainties of the calculated RFs are estimated from the percentage of fit between the observed and the calculated RF. Subsequently, well-constrained Vs structures of the crust and upper mantle are provided in the 2D profiles.

We applied the Delaunay triangulation interpolation method to estimate unknown velocities based on several known calculated velocities (Ping et al., 2009). The method uses three velocities at a time by assuming no points inside the circumference of any triangle. We applied this interpolation method as implemented in the GMT plotting software (Wessel et al., 2019) by triangulating and contouring the calculated velocity values to image the 2D velocity versus depth plots.

Crustal thickness and Vp/Vs ratio are estimated from the a priori known Vp value obtained from two-dimensional wide-angle seismic modelling from the EAGLE controlled-source survey (Mackenzie et al., 2005; Maguire et al., 2006) in the region. During our inversion, we calculated Vs values at 2-km-depth intervals. Again, we employed the mathematical model by Last et al. (1997) and Zhu and Kanamori, (2000) to get the Moho depth (H) at each station, where t_Ps − t_P is the time interval between the arrival of the direct P wave and the Moho Ps converted phase, and p is the average ray parameter calculated from RFs.

We computed observed and synthetic RFs at 17 stations where the degree of fit is between 77%–92% (at station OGOL and ASSE, respectively), as shown in Table 1. The range of degree of fit between observed and synthetic seismograms is similar to what is previously reported (70%–90%) in Ethiopia and Kenya by Dugda et al. (2005).

The current RFs are obtained with two clusters of range of back-azimuths of 30°–110^o and 185°–260^o (Figure 4). The first arrival spike is the direct incident P wave at the surface; however, the subsequent arrivals correspond to the partition of converted and reverberated phases (Figures 3 and 5).

For the two examples of observed RF in Figure 4, we present the RF of each event with the corresponding synthetic RF in Figure 5. In this model, the red-coloured RFs are the observed signals, whereas the blue colour shows the synthetic ones. The observed and synthetic RFs (Figure 5) show a high degree of fit for the Gaussian width parameters of α = 0.5 and α = 1.0.

Figure 6 indicates the 1D velocity models for the chosen two stations. These velocity models are calculated from the blue-coloured synthetic RFs shown in Figure 5. They are calculated in the depth range of 2–100 km. The blue-coloured nearly vertical line is the initial velocity model, which is assumed as a homogeneous half space with a Vs of ∼4.48 km/s, which is the value of most of the lithosphere in the ak135 velocity model (Kennett et al., 1995).

From the calculated 1D Vs models shown in Figure 6, the red-coloured 1D velocity value is the final and best fit calculated Vs model. From the models, the ASSE station, which is located on the Eastern Plateau shows very small heterogeneity in the upper and lower crust. However, station OGOL is located on the CMER floor and shows a heterogeneous velocity structure with a relatively high velocity of up to ∼4.6 ± 0.1 km/s in the upper crust and a relatively low velocity of as low as 3.4 ± 0.1 km/s in the lower crust.

For the remainder of the stations, we have shown the results in the form of the along- and across-rift profiles (Figure 7). Broadly speaking, the velocity models show a distinctive reduction in Vs in the mid to lower crust similar to that observed at OGOL (Figure 6), or a more regular increase in Vs with depth as observed at ASSE (Figure 6). Closer inspection for the stations along the rift shows that the velocity model varies considerably spatially with both styles of velocity structure observed in different places within the rift. In contrast, the across-rift profile shows that the stations on the Eastern Plateau have a velocity structure more similar to ASSE.

The upper to mid crustal high-velocity material (∼4–4.5 ± 0.1 km/s) observed in OGOL is also observed in the rift beneath YIRG, SHAS, KADO, and BESH stations for the depth range of 4–25 km (Figure 7). At these stations and beneath the observed high-velocity upper to mid crust, there is a relatively slow Vs (∼3.5 ± 0.1 km/s) for the depth range of ∼24–45 km. The slow velocity deep crust is commonly beneath normal upper-to-middle crust (4–4.3 km/s) such as beneath the JIMA station. In contrast, beneath JIRE, ODAS, ASSE, and SAGU stations (Figure 7), crustal Vs are relatively homogeneous.

Figure 8B1 shows the 2D Vs structure and Figure 8B2 the corresponding Vp/Vs ratio of the along-rift profile in the CMER obtained from the Delaunay triangulation interpolation method. Throughout the crust, the depth to particular velocity contours generally deepens with proximity to the major volcanic centres. This is especially pronounced in the 5- to 20-km depth range where the Vs increase significantly in regions in between the major volcanic centres. For example, beneath the two high topographic peaks (marked as Aluto and Tulumoye) observed in Figures 8A1,A2, there are slow velocity (<3.8 km/s) and high Vp/Vs ratio zones in the upper-mid crust. A similar slow velocity zone in the upper-mid crust is also observed beneath the Wondo-Genet remnant Mega caldera rim. Vs is generally slow (∼3.1–3.7 ± 0.1 km/s) in the lower crust beneath the CMER, with less spatial variation in velocities compared with that observed in the upper-mid crust. Generally, our findings are consistent with previous ambient noise tomography results showing the presence of slow S-velocity shallow crust beneath mega calderas, such as beneath Aluto and Tulu Moye, and slow Vs found more ubiquitously in the lower crust (Chambers et al., 2019).

Figures 9B1,B2 show the variations in Vs and Vp/Vs structure across the rift, respectively. In a similar fashion to the along-rift profile, the topmost ∼5 km of the upper crust of the across-rift profile is a very low seismic velocity (2.0–3.2 km/s) material. The border fault of the eastern side of the CMER is marked by a topographic step from ∼1,700 m in the rift to ∼2,700 m on the rift margin (Figure 9). Outside of the rift on the rift flank, we observe a fairly homogeneous crustal structure with a distinct lack of slow velocities in the lower crust. Instead, the seismic velocity mostly increases with depth. In addition, there is a sharp increase in seismic velocity at ∼45 km depth, where previous studies based on different methods showed this change as Moho discontinuity (e.g., Mahatsente et al., 1999; Mackenzie et al., 2005; Cornwell et al., 2006), as shown in Table 2. However, similar to the along-rift profile, within the rift on the across rift profile, we see a more heterogenous velocity structure. At 20- to 35-km depths, particularly slow Vs and high Vp/Vs ratios are found beneath the eastern part of the across-rift profile beneath the JIMA and OGOL stations. In this depth interval, the lowest velocities are found beneath the eastern side of the CMER spatially associated with the surface position of the WFB volcanic centres.

We discuss here the Vs and Vp/Vs structure of the rift obtained from our data analysis in the context of magmatic and tectonic extensional processes, and with the aid of a priori constraints of Vp ∼6.8 km/s (e.g., Dugda et al., 2005; Mackenzie et al., 2005; Maguire et al., 2006). We also compare our findings with constraints inferred from density and conductivity analysis conducted in the area. We use both one- and two-dimensional Vs profiles to interpret velocity variations in the lithosphere to answer basic questions about the nature of the crust and upper mantle when rifting modifies the lithosphere.

We use RF to delineate the Vs structure of the lithosphere beneath 17 stations in the CMER, which are arranged in two profiles along and across the rift valley. The observed low Vs (∼2–3 km/s) uppermost crust (<6-km depth) is interpreted as sedimentary and/or volcanic layers. Beneath the rift valley crust, Vs is heterogeneous laterally and with depth. In particular, slow Vs and high Vp/Vs ratio is localised beneath volcanic centres in the upper-mid crust but ubiquitously slow in the lower crust. The Vs and Vp are consistent with the presence of the small fraction (<5%) partial melt interpreted in previous seismic imaging studies of the lower crust. In addition, the velocity contrast is small between the lower crust and upper mantle in the rift. The results suggest that partial melt in the lower crust beneath magmatically active rifts might be more widespread than previously thought and is an important component of the magma plumbing system. In contrast, Vs is more homogeneous and faster beneath the Eastern Plateau, with a distinct and sharp velocity contrast observed between the crust and upper mantle at Moho, jointly indicative of very little crustal modification from magmatism.