Seismicity and Crustal Structure of the Southern Main Ethiopian Rift: New Evidence From Lake Abaya

The Main Ethiopian Rift (MER) has developed during the 18 Ma‐Recent separation of the Nubian and Somalian plates. Extension in its central and northern sectors is associated with seismic activity and active magma intrusion, primarily within the rift, where shallow ( < 5 km) seismicity along magmatic centers is commonly caused by fluid flow through open fractures in hydrothermal systems. However, the extent to which similar magmatic rifting persists into the southern MER is unknown. Using data from a temporary network of five seismograph stations, we analyze patterns of seismicity and crustal structure in the Abaya region of the southern MER. Magnitudes range from 0.9 to 4.0; earthquake depths are 0–30 km. VP/VS ratios of ∼ 1.69, estimated from Wadati diagram analysis, corroborate bulk‐crustal VP/VS ratios determined via teleseismic P‐to‐S receiver function H‐ κ stacking and reveal a relative lack of mafic intrusion compared to the MER rift sectors to the north. There is a clear association of seismicity with the western border fault system of the MER everywhere in our study area, but earthquake depths are shallow near Duguna volcano, implying a shallowed geothermal gradient associated with rift valley silicic magmatism. This part of the MER is thus interpreted best as a young magmatic system that locally impacts the geothermal gradient but that has not yet significantly modified continental crustal composition via rift‐axial magmatic rifting.

in southern MER , likely suggesting the crust has been less thermally modified by magma intrusion (e.g., Daniels et al., 2014). Consistent with the expectation that it is a less magmatically mature rift sector, GPS data reveal that southern MER strain rates peak along the western rift margin in the Lake Abaya region (Figure 1: Kogan et al., 2012). However, geophysical constraints on crustal structure and patterns of seismicity are lacking in this region, meaning the magma-poor hypothesis for the southern MER remains untested.
To address these issues, we have deployed a new passive seismic network of five broadband seismograph stations in the Lake Abaya region to place new constraints on seismicity in the southern MER. To understand better how seismicity is distributed throughout the southern MER crust, we also perform a receiver function study to determine Moho depth and bulk crustal / P S V V ratio following a modified version of the Zhu and Kanamori (2000) H- stacking technique (Ogden et al., 2019).  Corti et al. (2013). Focal mechanisms are from the Harvard CMT catalog. Magenta squares: thermal springs. Magenta stars: fumaroles. White circles: hypocenters from the studies of Keir et al. (2006) and Lavayssière et al. (2019). Blue circles: earthquake locations from the NEIC catalog. Inset: Broad tectonic setting of the Lake Abaya region in Ethiopia. Blue dots: earthquake locations from the NEIC catalog.

Tectonic Setting
The MER has developed over the past 18 Ma from the separation of the Nubian and Somalian plates (Wold-eGabriel et al., 1990;Wolfenden et al., 2004) (Figure 1). Extension began in the Miocene with rift basins defined initially by large offset, 60 km long border faults (Wolfenden et al., 2004). Through time, the degree of extension by magmatism is thought to have increased and manifested at the surface as the active volcanic centers, and systems of NNE-SSW trending, shorter length-scale, small-offset, intra-rift faults known as the WFB along the rift valley floor (Agostini et al., 2011;Ebinger & Casey, 2001;Mohr, 1967;Pizzi et al., 2006). The transition from border faulting to intra-rift localized extension is clearest in the northern MER, where the WFB defines the present-day rift axis (e.g., Ebinger & Casey, 2001). In the central MER, the WFB lies close to the eastern MER margin, while an additional volcanic belt, the Silti-Debre Zeyit Fault Zone (SDFZ) flanks the western margin (Rooney et al., 2007(Rooney et al., , 2011. While magmas rise quickly through the crust below the WFB before fractionating in the upper crust, the SDFZ magmatic plumbing system appears less evolved, with fractionation occurring throughout the crust (Rooney et al., 2007(Rooney et al., , 2011. South of 6.5-7.  5 N, in the southern MER, the morphological expression of the volcanic centers cut by young faults is less clear compared to the northern and central MER (e.g., . Significant slip on boundary faults is thought to have occurred during the Pliocene (e.g., Bonini et al., 2005), but axial deformation associated with the WFB remains near-negligible (e.g., Agostini et al., 2011;Corti et al., 2013;. Corroborating this structural geological observation, in contrast to the northern MER (Bilham et al., 1999;Ebinger & Casey, 2001), GPS data suggest that extension in the southern MER is focused along the Chencha escarpment of the western rift margin (Figure 1: Kogan et al., 2012). This suggests that strain is still localized to the border faults. Lake Abaya is positioned at the transition between central and southern MER deformation and is therefore an ideal study locale to understand early processes during the change from border faulting to more localized magma-assisted rifting.

Previous Geophysical Work
Current MER extension is approximately east-west oriented, with the rate of extension increasing from 4 mm/yr in the southern Ethiopian rift, to 6 mm in the northern Ethiopian rift (Kogan et al., 2012;Saria et al., 2013;Stamps et al., 2018). MER crustal structure has been constrained using a variety of geophysical techniques over the past two decades, particularly during the 2001-2003 Ethiopia Afar Geoscientific Lithospheric Experiment (EAGLE: Bastow et al., 2011;Maguire et al., 2003). Wide-angle controlled-source seismic profiling (e.g., Mackenzie et al., 2005;Maguire et al., 2006), supported by inversion of gravity data (Cornwell et al., 2006) and local earthquake tomography (Daly et al., 2008) demonstrate that dense, fast P-wavespeed, gabbroic intrusions are focused beneath Quaternary magmatic segments. According to the EAGLE along-axis controlled source profile, Moho depths vary from 26 km in the northernmost MER to 35 km at the southern end of the profile at Lake Hawassa in the Central MER . Moho depth constraints from teleseismic P-to-S receiver functions corroborate the wide-angle Moho depth observations and provide additional bulk-crustal / P S V V ratios. For example, at station E34 in Shashamene to the north of our study area (Figure 1), crustal thickness is 38.9 km . Further south in Arba Minch (ARBA: Figure 1), crustal thickness is 30 km (Dugda et al., 2005), with an associated bulk-crustal Poisson's ratio of 0.26 ( / P S V V ratio of 1.76) when crustal P V = 6.2 km/s: lower than the MER to the north, where / P S V V is generally 1.9-2.2 .
Although the EAGLE along-axis wide-angle seismic line reveals fast wavespeed, gabbroic intrusions along much of the Ethiopian rift, mean crustal P V reduces from 6.35 km/s in the north at ∼  10 N to only 6.1 km/s in the south approaching Lake Abaya , perhaps indicating a reduction in intrusive igneous material in the southern, less mature rift sectors. Corroborating this view, / P S V V ratios at mid and upper crustal levels constrained by local earthquake tomography (Daly et al., 2008) are lowest in the south near Lakes Langano and Shala (Figure 1), consistent with a more felsic upper crust than to the north. Shear-wave velocity constraints from ambient noise tomography (Chambers et al., 2019;Kim et al., 2012) are generally slow in regions where P V constraints are fast, likely indicating an increased sensitivity to melt over composition. With the caveat that our study area is on the edge of the resolved areas in the ambient noise studies, mid-to-lower crustal slow S V anomaly amplitudes are generally lower around Abaya than to the north of 7°N. In the central/northern MER, MT data (e.g., Whaler & Hautot, 2006) reveal enhanced OGDEN ET AL.

10.1029/2021GC009831
3 of 17 conductivity beneath the Boset magmatic segment at 20 km depth, interpreted as evidence for melt; in the uppermost upper crust (5 km depth), high conductivity material is more likely to be due to hydrothermal fluids than magmatism. InSAR and GPS reveal a number of persistently deforming volcanoes in the MER. For example, Albino and Biggs (2021); Biggs et al. (2011);Lloyd et al. (2018) document surface deformation in response to melt migration beneath the Corbetti volcano. Further south, below Aluto volcano at 7.  787 N, 38.  771 E, deep-seated high conductivity zones are lacking, implying ground deformation at the volcano is perhaps more likely the result of a shallow hydrothermal system (Samrock et al., 2015). Beneath the Tulu Moye geothermal prospect, Samrock et al. (2018) find evidence for a shallow crustal melt zone driving the hydrothermal system.
Little is known about the earthquake activity of Abaya and the nearby rift valley faults. The International Seismological Centre catalog reveals that since 1960, only 45 earthquakes have been reported within 200 km of the center of our network at 6.  75 N,  38.1 E with no obvious spatial patterns to note. Depths vary from 0 to 38.5 km; magnitude estimates range from 2.2 to 6.3. Of these, only six have associated focal mechanisms reported by the GCMT catalog ( Figure 1). This pattern of earthquakes is similar to previously compiled regional catalogs (Ayele & Kulhanek, 1997).
More is known about the microseismicity of the central and northern Ethiopian rifts, where dense deployments at Aluto (Wilks et al., 2017), Corbetti (Lavayssière et al., 2019), and Tulu Moye  show episodic swarms of densely distributed microseismicity in the uppermost crust (5 km depth) localized to the WFB and cross rift faults at these volcanic systems. This shallow microseismicity is interpreted as induced by localized hydrothermal fluid migration. This is set within a broader rift-wide pattern of seismicity with earthquakes mostly localized above 15 km depth on rift parallel Wonji and border fault systems (Ayele, 2000;Keir et al., 2006).
The mantle below the central and northern MER is known to be seismically slow (e.g., Bastow et al., 2008;Gallacher et al., 2016). Further, SKS splitting and surface-wave derived measurements of seismic anisotropy (e.g., Bastow et al., 2010;Kendall et al., 2005) are explained best by elongate melt inclusions in the MER lithosphere. The amplitude of uppermost mantle (75 km depth) shear wavespeed anomalies drops markedly south of ∼  7 N, while splitting delay times reduce from a maximum of  t = 3.12 s in the central MER to  t  1 s south of ∼  7 N (Kendall et al., 2006).

The Abaya Broadband Seismograph Network
In July 2019, we deployed five, 3-component Güralp CMG-6TD seismometers with a station spacing of 25-50 km, to the north of Lake Abaya ( Figure 1). Instruments recorded at 100 Hz until February 2020. Stations were deployed in secure health clinics and private settlements, resulting in an excellent data recovery rate of 98.8%. Analysis of power-spectral density plots (Figures S1 and S2) indicates that the network experiences generally high levels of noise in the period range 0.07-1.3 s where most body wave energy is anticipated: station YRBA, for example, plots particularly close to the Natural High Noise Model of Peterson (1993). For a typical local earthquake seismogram, we find that P-wave energy is found dominantly at 5 Hz; S-wave signals are typically lower frequency. We, therefore, opt to filter seismograms using a Butterworth bandpass filter with corner frequencies of 0.75 and 15 Hz to remove often-strong background noise at 1 Hz while retaining clear, un-distorted P-and S-wave signals.

Local Earthquake Hypocenter Determination
Continuous seismograms from the five Abaya stations were inspected manually for all possible earthquake signals, with suitable P-and S-wave phase arrivals picked accordingly on filtered seismograms. All picks were assigned a quality rating between 0 and 4 corresponding to visually assessed pick errors of 0.05, 0.2, 0.4, 1, and 9,999.9 s respectively. S-wave picks had generally poorer ratings, particularly for nearby earthquakes whose S-waves arrived within the P-wave coda. A quality rating of four was assigned in instances where P-or S-wave coda were identifiable, but first breaking energy could not be picked accurately. A total OGDEN ET AL.

10.1029/2021GC009831
4 of 17 of 313 accurate P-wave picks and 366 S-wave picks with quality ratings of 0-3 yielded 114 potential hypocenters (Table S1). An additional 42 earthquakes were detected but could not be located accurately.
A bespoke one-dimensional velocity model for our study area ( Figure 1) could not be calculated from our small aperture network. Instead, we used the P-wave velocity model derived by Lavayssière et al. (2019), whose focus was Corbetti volcano, 50 km to the north of our study region ( Figure 1). Shear velocities were determined assuming a / P S V V ratio determined following Wadati and Oki (1933): for earthquake-station pairs with both P-wave and S-wave phase arrival times, we plot the P-arrival time against the delay time between the P-wave and S-wave arrivals ( Figure 2), calculated from preliminary hypocenters estimated assuming  Figure 2. Earthquake hypocenters were located using the NonLinLoc method of Lomax et al. (2009). The non-global grid-search area sampled a 300  300 km area centered on northern Lake Abaya, extending to 50 km depth with 0.2 km node spacing in all directions. We used the weighted equal differential time likelihood function (Font et al., 2004) to estimate the posterior density function because it is less sensitive to outliers than the traditional least squares, L2 norm likelihood function (Lomax et al., 2009). NonLinLoc utilizes a three-dimensional, oct-tree approach to efficiently locate the hypocenter, with the maximum likelihood position in time and space chosen as the hypocenter. The best located hypocenters typically have spatial location and depth errors of 10 km. Two distinct geographical clusters of earthquakes are detected, identified herein as "Cluster 1" in the north of the study region near Duguna volcano, and "Cluster 2" in the south of the study region to the west of Lake Abaya.

Earthquake Magnitudes
Due to the relatively small number of located earthquakes in our study, we do not calculate a bespoke local magnitude scale for our earthquake catalog. Earthquakes located outwith our focused study region are biased inherently toward larger magnitudes, and sample different fault networks, potentially skewing the Gutenberg-Richter relationship of the study region. Hence these earthquakes are not included in our frequency-magnitude analysis. Accordingly, we calculate local magnitudes ( L M ) using only the 88 earthquakes (Table S2) located within the geographical region outlined in Figure 1. We use the scaling equation from Keir et al. (2006): where WA A is the maximum zero-to-peak amplitude recorded on the horizontal seismogram corrected to a Wood-Anderson response, r is the epicentral distance in kilometers and C is the individual station correction which is taken to be zero for all stations for this study. The reported L M for an individual earthquake is the average of the two horizontal component magnitudes. The magnitude of completeness cut-off for the catalog ( C M ), below which earthquakes are under-sampled owing to the network's insufficient station coverage, is calculated using the maximum curvature method (Wiemer & Wyss, 2000) to be C M = 1.4, a reasonable choice of method for a small, likely broad-peak frequency-magnitude distribution (Roberts et al., 2015). However, the maximum curvature method often underestimates where M is the mean local magnitude (for  L C M M ) and M is the magnitude binning width. We calculate an error ( b ) of 0.08 associated with this b-value following Shi and Bolt (1982): where n is the number of earthquakes in the catalog and i M is the local magnitude of the ith earthquake. We recognize that the small sample size of earthquakes is not ideal, which combined with the so-called broadpeak distribution (Roberts et al., 2015), means that we interpret our b-value and C M with caution.

Focal Mechanism Analysis
Double couple focal mechanism solutions were determined following Snoke (2003) using first motion polarities of P (picked on vertical component seismogram) and SH (picked on tangential component seismogram) phases. The Snoke (2003) method uses a grid-search approach around the focal sphere (at  2 intervals) to calculate acceptable nodal plane orientations. Earthquakes with nine or more picked phases were analyzed on raw, unfiltered, velocity seismograms to remove the effect of processing artifacts, with only confident P and SH arrivals picked for analysis. Take off angles and earthquake-station azimuths were determined using our hypocenter solutions. We allow zero polarity errors owing to the high level of confidence during polarity picking. Because our network is small and the relatively distant earthquakes in cluster 2 are constrained with a large azimuthal gap, our analysis yielded only four acceptable quality focal mechanism solutions. Table 1 describes the final fault plane solutions (see also Figures S6 and S7 for full focal mechanism solutions).

Receiver Function H- Stacking
Receiver functions capture P-to-S wave conversions at velocity contrasts in the crust and mantle, below the seismograph station, that are recorded in the P-wave coda of teleseismic earthquakes (e.g., Langston, 1977;Phinney, 1964). The travel times of the Moho P-to-S conversion (Ps) and subsequent crustal reverberations (PpPs,  PpSs PsPs) are then routinely studied to provide information on the crustal thickness (H) and / P S V V ratio () via the H- stacking procedure outlined by (Zhu & Kanamori, 2000   Note. Errors are derived from the standard deviation of the possible solutions for a given earthquake.

Table 1 Focal Mechanism Solutions (Probable Focal Plane Only) and Hypocenter Information for the Respective Earthquake
where N is the number of receiver functions; 1 w , 2 w , 3 w are stacking weights, ( ) j i r t are the amplitudes at the arrival times for each of the arrivals evaluated. 1 t , 2 t , and 3 t are the travel times of the Ps, PpPs, and  PsPs PpSs phases, respectively: OGDEN ET AL. Abbreviations: ACE, amplitude comparison estimate; CCC, cross correlation coefficient; SNR, signal-to-noise ratio.

Table 2
The 10 Criteria Used for Determining Result Quality Following (Ogden et al., 2019). A Criteria is Passed if the  where p is the ray parameter and P V and S V are bulk crustal P-and S-wave velocities, respectively. The optimal combination of H and  is determined where s(H,) is maximized. The choice of stacking input parameters including average crustal P-wave velocity ( P V ), the assigned phase weightings ( 1 w , 2 w , and 3 w ) and the choice of stacking strategy: linear versus phase-weighted (e.g., Schimmel & Paulssen, 1997), can each exert strong control on the final result. So too can data parameters such as receiver function frequency content and the subset of receiver functions stacked during H- analysis. To overcome this parameter sensitivity issue, we follow the modified H- stacking approach of Ogden et al. (2019) in which cluster analysis selects a final solution from 1,000 individual H- results, each calculated using randomly selected receiver functions and H- input parameters. Ten quality control criteria that variously assess the final numerical result, the receiver function data set, and the extent to which the results are tightly clustered, are used to assess the reliability of H- stacking at a station (Table 2; Figure 4). These criteria include: 1. The Amplitude Comparison Estimate (ACE): the amplitude ratio of the predicted Ps phase (defined by the final H- solution) and the root mean square (RMS) amplitude of the receiver function between  1 2 t s and  2 2 t s. A simple layer over half-space model with a sharp crust-mantle transition will exhibit an impulsive high amplitude Ps phase compared to the general receiver function signal. A gradational crust-mantle transition will produce a lower-amplitude Ps phase and thus a lower ACE. Complex intra-crustal structure will yield additional P-to-S conversions between  1 2 t s and  2 2 t s, which also lower the ACE and render H- analysis less reliable (Figure 4). 2. Signal-to-noise ratio (SNR): the amplitude ratio of the predicted Ps phase (defined by the final H- solution) and the RMS amplitude of 8 s of pre-P-arrival noise. A large Ps amplitude from a sharp crust-mantle transition yields a higher SNR than for a gradational Moho (e.g., due to lower crustal intrusions) (Figure 4). 3. The cross correlation coefficient (CCC): of all possible pairs of receiver functions calculated with the same frequency, for each different frequency of receiver function. The CCC tests the effect of noise and backazimuthal variations in structure at a seismograph station. Stations with a high CCC yield more stable H- estimates.
Seismograms for all earthquakes with M4+ (from the NEIC earthquake catalog) within  30 - 90 epicentral distance were analyzed. All seismograms were filtered using a Butterworth bandpass filter with corner frequencies of 0.04 and 3 Hz. Seismograms were visually checked for pre-P-arrival noise and discarded if unsuitable. Receiver functions were calculated for seismograms deemed acceptable, using the Extended Time Multitaper Receiver Function method of Helffrich (2006) for a range of frequencies (0.4-2.0 Hz at 0.1 Hz intervals). Receiver functions were then manually inspected and those with high pre-P-wave noise were discarded. Those remaining were analyzed with the modified H- stacking method of (Ogden et al., 2019) to obtain Moho depth (H) and bulk-crustal / P S V V ratio (). In addition to the five Abaya stations, station ARBA which operated between 2000 and 2001 as part of the Ethiopian Broadband Seismic Experiment (Nyblade & Langston, 2002) was also analyzed. A P-wave crustal velocity ( P V ) range of 6.0-6.6 km/s was used in this study, to span the 6.2 km/s average P V estimated by the EAGLE along-rift wide-angle seismic profile  that sampled crustal structure proximal to this study region in the central and northern MER. Stations ARGA, ABYA, BEDA, and YRBA have insufficient high quality receiver functions to obtain an H- result, likely due to a complex upper crustal velocity structure associated with sediments and volcanic related deposits: P-to-S converted energy from such interfaces is known to disrupt successful H- analysis (e.g., Ogden et al., 2019). However, HMBA and ARBA yielded good quality results when limiting the results to using lower frequency receiver functions: 1.7 and 0.9 Hz, respectively (Figures 5 and 6). Moho depth is 37  0.5 km with

A New Seismicity Catalog for Southern Ethiopia
The overall pattern of the 114 earthquakes in the full catalog ( Figure 7; Table S1) reveals that seismicity is mainly clustered to the northwest of Lake Abaya on the southern flank of Duguna volcano (cluster 1), and west of Lake Abaya (cluster 2). Both clusters are along the western border fault of the MER; magnitudes are in the range 0.8-4.0. Earthquakes are deeper in cluster 2 than cluster 1 (Figure 7 inset). Analysis of hypocenters of the 29 earthquakes with horizontal errors of 10 km (Figure 8; Table S3) indicates that all but four of the most accurately located earthquakes occur in these two clusters of earthquakes with low levels of seismicity outwith these two clusters. When viewed temporally in Figure 9, cluster 1 comprises two discrete earthquake swarms that occurred in July and November of 2019, seismicity in most of cluster 2 occurred in January 2020. Therefore, clusters 1 and 2 both display earthquake swarms. Elsewhere seismicity is scattered spatially and temporally.
To interpret the cluster geometry, in Figures 10 and 11 we focus on interpreting the best located earthquakes with horizontal errors of 10 km. Cluster 1 depths are generally 8-12 km depth below the southern flank of Duguna Volcano, in and around a zone of short length-scale (5-10 km) faults (Figures 8 and 10). Such a pattern is observed elsewhere in the Ethiopian rift where, for example, earthquakes below the Fentale-Dofen and Aluto magmatic segments are found at very similar, upper-crustal depths (e.g., Keir et al., 2006).
With the caveat that earthquakes in cluster 2 are outside the network and have a high azimuthal gap of   180 , we see hypocenters fall along or near to the 25 km long faults bounding the western margin of OGDEN ET AL.

10.1029/2021GC009831
11 of 17 12 of 17 Figure 8. Earthquake distribution for the 29 earthquakes hypocenters with horizontal location errors 10 km. The same two major clusters are identified as in Figure 8. Four reliable focal mechanisms are plotted and shaded according to depth. Black lines: faults from Corti et al. (2013). Top left inset: Earthquake distribution as a function of magnitude for these hypocenters, colored by cluster location. Bottom right inset: Earthquake depth distribution for these hypocenters, colored by cluster location. Figure 9. All 114 earthquakes plotted as a function of time and latitude. A low background level of seismicity occurred along the length of the western fault system throughout the 8 months of observations. Three temporal swarms of earthquakes are identified, two swarms are 4 months apart at ∼  6.92 N in cluster 1. The second of these swarms were initiated by a L M = 3.8 earthquake. The swarm in January 2020 at ∼  6.31 N is in the region of cluster 2. the MER; most earthquake depths are 16-24 km (Figure 11: C-C′) and are interpreted to be in the mid-tolower crust based on our receiver function-derived estimates of Moho depth (H = 37 km at station HMBA and H = 30 km at station ARBA: Figures 5 and 6, respectively). A similar depth range of earthquakes has been reported in southern Ethiopia in continent-scale seismicity studies that constrain earthquake depths using teleseismic and regional seismic data (e.g. Craig & Jackson, 2021). Focal mechanisms show normal faulting ( Figure 11), with NNE-SSW striking nodal planes parallel to mapped faults at the surface along the MER's western margin. The mid-to-lower crustal depths are robust, with consistently tight probability-density functions (e.g. Figures S3 and S4), and are demonstrably deeper than the earthquakes in cluster 1 (Figure 10).

Implications for Extension and Magmatism in the Southern MER
Cross sections through cluster 1 (Figure 10), which lies at the boundary between the central and southern MER (Figure 1), do not unambiguously delineate a simple, single fault plane associated with the western border fault system of the MER. Minissale et al. (2017) found that the chemistry of thermal springs and fumaroles discharging naturally in the area immediately to the north of Lake Abaya is dominated by a 2 CO -rich gas phase and discharges along the active faults bordering the western edge of the MER. Cluster 1 is therefore readily explained by deformation associated with upper crustal magmatic fluid migration, akin to the central and northern MER (e.g., Keir et al., 2006). South of ∼  6.9 N, however, within the southern MER, a striking observation in Figures 8 and 11 is the clear association of seismicity with the western border fault system of the MER to 15 km depth. This corroborates GPS observations that southern MER strain rates are OGDEN ET AL.

10.1029/2021GC009831
13 of 17 highest along the Chencha escarpment of the western rift margin (Figure 1: Kogan et al., 2012). According to our H- stacking-derived Moho depths at stations HMBA (37  0.5 m: Figure 5) and ARBA (30.3  1 km: Figure 6), southern MER seismicity persists to lower crustal depths, in line with estimates of effective elastic plate thickness ( t e ), which increase abruptly from t e = 8 km in the central MER to t e = 17  2 km in the Abaya region . Furthermore, continuous faults of length 25 km (Corti et al., 2013), along which deeper earthquakes may be expected to manifest, are mapped to the west of Lake Abaya in the vicinity of cluster 2 hypocenters. Accordingly, with the caveat that it is an observation based on a very short earthquake catalog, our b-value of 0.69  0.08 is consistent with a fault network dominated by longer length-scale, larger offset faults than in the central and northern MER. b-values in the central and northern MER are consistently larger in various regional and volcano-specific studies where the locus of extension has now shifted from mid-Miocene border faults, to rift-axial, small-offset, short length-scale faults associated with the WFB (Ebinger & Casey, 2001). For example, b = 1.13  0.05 (central and northern MER: Keir et al., 2006); b = 1.32  0.07 ; b = 1.40  0.14 (Wilks et al., 2017) (Figure 3). Our observations point strongly toward a hypothesis of continued brittle strain on border faults in this sector of the Ethiopian rift. However, while our earthquake depths and locations are consistent with them being of a purely tectonic origin, the swarm-like nature of much of our seismicity, coupled with independent measurements of 2 CO release along the fault systems (e.g., Minissale et al., 2017), indicates that fluid transport along the faults may also contribute to the observed seismicity.
Our observations join a growing body of evidence that earthquake depths in East Africa generally increase southwards, mirroring increases in plate thickness and strength (Craig et al., 2011). Most recently, abrupt changes in earthquake depths near major terrain boundaries have been cited as evidence for sharp changes in crustal composition (Craig & Jackson, 2021). To this end, that earthquake depths are the shallowest near Duguna volcano perhaps suggests that short length-scale variations in brittle layer thickness may result from variations in geothermal gradient associated with rift valley silicic magmatism. Compressible gas OGDEN ET AL.

10.1029/2021GC009831
14 of 17 Figure 11. Left panel-The southern cluster of seismicity to the west of Lake Abaya is aligned north-south proximal to the major western border fault, but is located in the mid-to-lower crust.Black lines: faults from Corti et al. (2013). Right panel: C-C′ fault-perpendicular cross-section. phases in the hydrothermal and volcanic systems, and basin fill of sediments, ash, and evolved volcanic rocks may also contribute to the low / P S V V (1.69) observed from our local earthquake Wadati analysis, which dominantly samples the upper crust. Such a conclusion is similar to other volcanic regions in the Ethiopian rift such as Corbetti and Tulu-Moye (e.g., Greenfield et al., 2019). However, receiver function-derived bulk-crustal / P S V V ratios (1.7) at stations HMBA and ARBA (Figures 5 and 6) imply the southern MER crust is relatively unaffected by voluminous mafic intrusion and/or present-day melt.
The combined constraints of seismicity distribution and depth, when reviewed in light of / P S V V ratios from our seismicity and receiver function analysis point toward the Abaya region being a young magmatic system that impacts the geothermal gradient locally but has not yet significantly modified continental crustal composition. Consistent with this, at greater depth, the lithospheric mantle below the central and northern MER is known to be seismically slow (e.g., Bastow et al., 2008;Gallacher et al., 2016). Further, SKS splitting and surface-wave derived measurements of seismic anisotropy (e.g., Bastow et al., 2010;Kendall et al., 2005) are explained best by elongate melt inclusions in the MER lithosphere. The amplitude of uppermost mantle (75 km depth) shear wavespeed anomalies drops markedly south of ∼  7 N into the Abaya region, while splitting delay times reduce from a maximum of  t = 3.12 s in the central MER to  t  1 s south of ∼  7 N (Kendall et al., 2006). Considered in light of these observations, our results point strongly toward continued border faulting during the initiation of in-rift magmatism. This, in turn, may suggest that the localization of magmatism guides the in-rift migration of faulting.

Conclusions
Using a new network of five seismometers, we have constrained bulk crustal structure (Moho depth and / P S V V ratio) and seismicity in the Lake Abaya region of the southern MER. We constrain 88 hypocenters, which provide clear evidence for earthquake clusters along the western margin of the southern MER, including below the southern flank of the Duguna volcano. In contrast to the more evolved northern MER, rift axial seismic activity is generally lacking in the southern MER. Magnitudes range from 0.9 to 4.0, while earthquake depths are 0-30 km. / P S V V ratios of 1.69, estimated from Wadati diagram analysis, corroborate bulk-crustal values determined via teleseismic P-to-S receiver function H- stacking analysis and reveal a relative lack of mafic intrusion compared to the MER rift sectors to the north. The spatial association between seismically active faults, Duguna volcano, and the numerous hot springs and fumaroles in the region, suggests some earthquakes may be induced by fluid migration; deeper events reveal the ongoing importance of border faulting. The Abaya region is thus interpreted best as a young magmatic system that locally impacts the geothermal gradient but that has not yet significantly modified continental crustal composition via rift axial magmatic extension.

Data Availability Statement
Data are archived with the IRIS Data Management Center under the network code 1Q. Seismograms were analyzed using the Seismic Analysis Code (Helffrich et al., 2013). Hypocenters were determined using the NonLinLoc method of Lomax et al. (2009), available from http://alomax.free.fr/nlloc/and pick information can be found in supporting information. Figures were prepared using the Generic Mapping Tools (Wessel et al., 2019).