Modelling S-Wave Velocity Structure Beneath the Central Main Ethiopian Rift Using Receiver Functions

We applied the receiver function (RF) technique on high-quality teleseismic earthquake data recorded by the RiftVolc broadband network from February 2016 to October 2017. We calculate RFs at 17 stations, which are inverted to estimate Vs, and Vp/Vs structure beneath the Central Main Ethiopian Rift and the Eastern plateau. The observed slow S-wave velocity (Vs) in the uppermost crust (<6 km depth) is interpreted as sedimentary and/or volcanic layers. Beneath the rift valley, crustal Vs is heterogeneous both laterally and with depth. In particular, slow Vs (∼2–3 km/s) is localised beneath volcanic centres in the upper-mid crust but ubiquitously slow in the lower crust with Vs as low as ∼3.5 km/s common. The slow lower crust is associated with high Vp/Vs ratios of ∼1.9–2.0. The Vs and Vp are consistent with the observed seismic velocities, and interpreted the presence of the small fraction (<5%) of partial melt from previous seismic imaging studies of the lower crust. In addition, the velocity contrast is small between the lower crust and upper mantle. The results suggest that partial melt in the lower crust beneath magmatically active rifts might be more widespread than previously thought and an important component of the magma plumbing system. In contrast, Vs is far more homogeneous and faster beneath the Eastern Plateau, with a distinct velocity contrast between the crust and upper mantle suggesting less crustal deformation than what is observed beneath the central rift zone.


INTRODUCTION
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  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 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 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).
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 enechelon 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 3 ) 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).

Data
The data were acquired from the RiftVolc temporary network project that was conducted from February 2016 to October 2017 and recorded by three-component broadband Guralp CMG-6TD and Guralp CMG-ESPCD seismometers with a 50-Hz sampling rate. We downloaded the teleseismic waveform data and instrument responses of the RiftVolc project data archived at the Incorporated Research Institutions for Seismology (IRIS) Data Management Center (DMC).
To constrain the Vs and Vp/Vs structure beneath 17 stations, which are deployed along and across the CMER 60 teleseismic earthquakes with magnitudes, Mw ≥ 6 and source-to-receiver epicentral distances between 30°and 90°( Figure 2) were chosen. However, after calculating the RFs, only 8-38 signals were selected per station based on the percentage of signal power fit (Table 1).
Data were processed in SAC format. We applied a cosine taper function for the P-waveform signal for a length of 50 s (10 s before and 40 s after the onset of the P-wave arrival) before computing the RFs. To reduce the influence of low-frequency noise on the RFs, all the signals were filtered with a Butterworth bandpass filter of between 0.01 and 5 Hz to ensure the stability of the RFs and to avoid aliasing when decimating the data. Finally, each three-component signal was reviewed to remove signals that contained low signal-to-noise ratios and/or when any of the three components were not recorded properly due to instrument malfunction.

METHODS
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 "prewhitening" (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 wideangle 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-kmdepth 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.

RESULTS
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).  Mahatsente et al. (1999) FIGURE 3 | This is a sampled receiver function (RF), which is calculated from the deconvolution of the horizontal from the vertical component. The deconvolution is processed from the data collected by the ASSE station from an earthquake coming from a~178°azimuth. The diagram shows the direct P and the converted Ps and the PpPs multiples.
Frontiers in Earth Science | www.frontiersin.org March 2022 | Volume 10 | Article 773783 From the calculated 1D Vs models shown in Figure 6, the redcoloured 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-tomiddle 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.

DISCUSSION
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. The ASSE station located on the eastern plateau, which is the best fit RF of this study, whereas station OGOL is located in the central MER (CMER) rift margin near the eastern plateau, which has the largest misfit of all the RFs. Blue-coloured RFs in the background are synthetic, whereas the red-coloured RFs on top are the observed RFs. P represents the direct primary wave, and Ps is the converted phase at the Moho. The numbers on the left show the station name, Gaussian width parameter, percentage of fit, and the applied ray parameters for the specified RF. The numbers to the right of the RFs are the occurrence time of the earthquakes. In (A), two signals and in (B) six signals are calculated twice for the Gaussian width parameters of 0.5 and 1.0.

S-wave velocity structure within the rift
The slow velocity (2-3 km/s) imaged at 2-6 km depth is similar to the proposed Vs of~1.9-2.8 km/s typical of layered sediments (Benoit et al. (2006). This is also in good agreement with the work of Chambers et al. (2019), which interprets a similarly low velocity at the topmost upper crust as sedimentary and/or FIGURE 6 | Panels (A) and (B) are the two representative velocity models for the ASSE and OGOL stations. The nearly vertical start. mod is an initial half space velocity model derived from the ak135 global velocity model (Kennett et al., 1995) and the end. mod is the final and best-fit velocity model. The tmpmod96. xxx are the calculated velocity models from the relatively less fit RFs during an inversion. Frontiers in Earth Science | www.frontiersin.org March 2022 | Volume 10 | Article 773783 8 volcanic layers. This result agree with the interpretation of Cornwell et al. (2006), which interprets the existence of an upper crustal low-density (2,380 kg/m 3 ) layer that represents interspersed volcaniclastics, lava flows, and lacustrine sediments within the rift valley (Wolfenden et al., 2004). In support of this interpretation, the low Vs of the uppermost crust extends to greatest depths within the rift valley than outside of it ( Figure 9).
Both profiles shown in Figure 7 represent significant variations in the 1D velocity models. In particular, a number of the stations show elevated seismic velocity at 6-25 km, while others are less fast. When stations are organised spatially from NW to NE in Figure 8A1, the spatial variability of this is clearer. Typically, along rift, in-between the magmatic centres (such as beneath station YIRG, SHAS, KADO, BESH, and OGOL), the high Vs (~4-4.5 km/s) is present in the upper-to-mid crust. The across-rift profile in Figure 9 shows that these regions of higher Vs in the upper/mid crust are localised beneath the Wonji Fault Belt. The high seismic velocities coupled with their Wonji Fault belt position favours an interpretation of their origin being a solidified mafic intrusion, an interpretation in line with previous seismic imaging (Chambers et al., 2019), and spatially match high positive Bouguer anomalies constrained in gravity studies (Mahatsente et al., 1999;Tiberi et al., 2005;Cornwell et al., 2006). The average slow-velocity (~3.5 ± 0.1 km/s) regions at 24-45 km depth may represent a less mafic modification of a normal continental crust of Vp/Vs <1.85 (Zandt and Ammon, 1995), or a more complex modification from felsic intrusion, and/ or presence of partial melt with a Vp/Vs value of >1.9.
Beneath 25-km depth in the lower crust, the 1D models show that the majority of seismic stations show a reduction in Vs in the lower crust to 3.1-3.7 ± 0.1 km/s (Figure 7). Figure 8 shows that this feature is spatially ubiquitous. There is some spatial variability in the magnitude of the velocity inversion ( Figures  8A1,A2), with a hint that the most pronounced slow Vs regions in the lower crust are beneath the volcanic centres, such as Aluto, although this pattern is not particularly clear elsewhere. These regions of slow Vs correlate to high Vp/Vs of~1.9-2.1. The observation of slow Vs and high Vp/Vs in the lower crust in the rift valley is consistent with the ambient noise tomography by Chambers et al. (2019), which shows that the slowest velocities for all depths within the MER range from 3.28 ± 0.01 km/s at 10-km depth to 3.83 ± 0.01 km/s at 40-km depth. The magnitude of the slow Vs at this depth range, combined with the high Vp/Vs, is consistent with previous deep crustal imaging studies, which combined interpret between 0.5% and 5% partial melt (e.g., Chambers et al., 2019). More tightly constraining melt fraction from the seismic velocities alone is difficult since Vp and Vs measurements are potentially explainable by either lower melt fractions aligned vertically as dikes or more elevated melt fractions aligned as sills (Paulatto et al., 2010;Paulatto et al., 2012;Paulatto et al., 2019;Dvorkin, 2020). However, dominance of horizontal sill-like melt alignment is favoured by inversions for radial anisotropy derived from surface waves for the MER (e.g.,   Chambers et al., 2021), and petrological models of the deep crustal magma plumbing system globally (e.g., Annen et al., 2006). The interpretation of partial melt in the lower crust is also supported by high conductivities in the lower crust observed in crustal-scale MT studies at comparable depths (Whaler and Hautot, 2006).
Profile AA′ in Figure 8 shows variations in the Vs structure, which provides insights into the crustal-scale magma plumbing system. The ubiquitous slow Vs suggests a diffuse interconnected melt-rich lower crust beneath most of the rift valley, with potentially higher melt concentration beneath the volcanic centres. In contrast, in the upper half of the crust, slower Vs beneath the volcanic centres, with anomalously fast Vs in between the volcanic centres, is consistent with volcanic segment-centred melt supply, in which subvolcanic melt reservoirs focus and store melt, which is delivered episodically mafic intrusion along the rift axis. Such an upper crustal plumbing system has been proposed in Afar on the basis of episodic segment-centred fed dyke intrusions observed with InSAR and seismicity (Keir et al., 2009;Barnie et al., 2016). Here in the MER, a similar subvolcanic plumbing system is consistent with the seismic velocity structure of the upper/mid crust. In addition, in the MER, observations of such rifting episodes are lacking, with geodetic observations of magma-related ground deformation being restricted to volcanic centres such as Aluto and Tulu Moye (Biggs et al., 2011;Albino and Biggs, 2021). Similarly, localised subvolcanic pockets of melt beneath the volcanic centres (Gleeson et al., 2017) suggest localised longer-lived magma bodies in the shallow crust of the volcanic centres. However, our seismic imaging of the deeper crust suggests that the distribution of melt in the lower crust might well be widespread and enable significant melt transport along rift.

S-wave velocity structure of the Eastern Plateau
In contrast to the rift valley floor, the Vs structure beneath the Eastern Plateau is remarkably homogeneous ( Figure 9B1). In addition, the distinctive increase in Vs at~45 km depth, is remarkably similar to constraints on the Moho depth computed in our study, consistent with a previous wide-angle active source, and passive source RF studies (Mackenzie et al., 2005;Maguire et al., 2006), and adds support observations that the Moho beneath the Eastern Plateau is a sharp and distinctive seismological boundary (e.g., Ogden et al., 2019). This profile shows a smooth transition toward the shoulder compared with the western plateau margin in which sharp lateral contrast between plateau and rift is observed (Chambers et al., 2019). Limited heterogeneity of the crustal and mantle structure beneath the Eastern Plateau is typical of regions of stable continental crust with limited history of deformation and modification by magmatism (e.g., Thompson et al., 2010;Youssof et al., 2013, Youssof et al., 2015. The strong contrast in velocity structure from the Eastern Plateau into the rift (Figure 9) is in sharp contrast to the conjugate side of the rift valley, with the Western Plateau showing evidence for significant magmatic modification (e.g., Mackenzie et al., 2005;Chambers et al., 2019), indicating strong asymmetry to the rifting process. The lack of evidence for magmatic modification of the crust beneath the Eastern Plateau also favours a model of dynamic uplift from a deep-seated asthenospheric anomaly (e.g., Sembroni et al., 2016), as opposed to uplift being compensated by crustal magmatic additions (e.g., Keranen et al., 2009;Chambers et al., 2019).
Our study reveals new important insights regarding the variability in crustal structure and melt fraction on a local scale beneath the volcanic regions of the MER. The results demonstrate the continued need for more future efforts to understand crustal structure and distribution of partial melt in the wider sense beneath and near the East African rift. We would like to point out the need to have more international collaboration-although we would imagine that long-term and sustainable research in Ethiopia really needs local scientists to lead the way.

CONCLUSION
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.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: https://ds.iris.edu/wilber3/find_event.

AUTHOR CONTRIBUTIONS
BK customized the Computer Program for Seismology software and wrote some essential scripts. All the three authors developed the concept of this paper. BK selected and processed the teleseismic signals, conducted the modelling of all datasets, and led the writing of the paper. All authors contributed to the write up, discussion, and interpretation of the result of the paper.

FUNDING
The research work is sponsored by Addis Ababa University. The project was funded by the Natural Environment Research Council under NERC Grant NE/L013932/1. The publication charges were covered by the University of Southampton, and the laptop used to conduct the analysis was purchased using GCRF-UKRI funding from Ian Bastow, Imperial Collage London.