Body‐Wave Tomographic Imaging of the Turkana Depression: Implications for Rift Development and Plume‐Lithosphere Interactions

The Turkana Depression, a topographically subdued, broadly rifted zone between the elevated East African and Ethiopian plateaus, disrupts the N–S, fault‐bounded rift basin morphology that characterizes most of the East African Rift. The unusual breadth of the Turkana Depression leaves unanswered questions about the initiation and evolution of rifting between the Main Ethiopian Rift (MER) and Eastern Rift. Hypotheses explaining the unusually broad, low‐lying area include superposed Mesozoic and Cenozoic rifting and a lack of mantle lithospheric thinning and dynamic support. To address these issues, we have carried out the first body‐wave tomographic study of the Depression's upper mantle. Seismically derived temperatures at 100 km depth exceed petrological estimates, suggesting the presence of mantle melt, although not as voluminous as the MER, contributes to velocity anomalies. A NW–SE‐trending high wavespeed band in southern Ethiopia at < 200 km depth is interpreted as refractory Proterozoic lithosphere which has likely influenced the localization of both Mesozoic and Cenozoic rifting. At < 100 km depth below the central Depression, a single localized low wavespeed zone is lacking. Only in the northernmost Eastern Rift and southern Lake Turkana is there evidence for focused low wavespeeds resembling the MER, that bifurcate below the Depression and broaden approaching southern Ethiopia further north. These low wavespeeds may be attributed to melt‐intruded mantle lithosphere or ponded asthenospheric material below lithospheric thin‐spots induced by the region's multiple rifting phases. Low wavespeeds persist to the mantle transition zone suggesting the Depression may not lack mantle dynamic support in comparison to the two plateaus.

primary mechanism for strain accommodation throughout the crust (e.g., Cornwell et al., 2006;Mackenzie et al., 2005;Oliva et al., 2019) and mantle lithosphere (e.g., Bastow et al., 2010;Kendall et al., 2006;Tiberi et al., 2019). However, problematic in this picture is the low-lying, broadly rifted, and as-yet poorly studied Turkana Depression, separating the elevated Ethiopian and East African plateaus (Figure 1a). Precisely how Cenozoic strain has developed here is uncertain, not least because the Turkana Depression, herein referred to as Turkana, is host to a failed, NW-trending Mesozoic rift system (Figure 1a; e.g., Reeves et al., 1987), parts of which were below sea-level in Miocene times (Wichura et al., 2015). Geodetic data indicate 40%-100% of present-day strain is localized across faults and eruptive volcanic centers below Lake Turkana, meaning that southern Turkana may be in a state of incipient strain localization at crustal depths (Knappe et al., 2020). Whether this pattern is mirrored at mantle lithospheric depths, remains unresolved, including whether or not plate-scale magma-assisted extension is now underway. Strain patterns in northern Turkana and the southwestern Ethiopian plateau are still poorly resolved.
interactions (e.g., Corti et al., 2019;Ebinger et al., 2000). For example, the MER adopts its most distributed character in southern Ethiopia which also marks the northernmost limit of Mesozoic extension (Figures 1a and 1b). Whether this broadening is due to inherently different lithosphere in southern Ethiopia or whether complex MER-Eastern Rift linkage kinematics are at play, is unknown. Improved constraints on Turkana's lithospheric seismic structure are thus essential to our understanding of Cenozoic strain development.
Another research question exemplified by the East African Rift (EAR) concerns the cause of Cenozoic flood basalt magmatism (Figure 1a), specifically whether it is related to elevated mantle temperatures resulting from a single or multiple mantle plumes (e.g., George et al., 1998;Rogers et al., 2000;Rooney et al., 2012). The two-plume hypothesis would provide an elegant explanation for the existence of the elevated Ethiopian and East African plateaus and the interconnecting, low-lying nature of the Turkana Depression. However, Turkana's complex history of multiple rifting phases (Figure 1a; e.g., Purcell, 2018;Reeves et al., 1987), allows for the possibility that the plateaus are instead a continuously uplifted region attributed to a single superplume, with the low elevations explained by crust stretched during Mesozoic and Cenozoic times rather than a lack of dynamic support (e.g., Benoit, Nyblade, & Pasyanos, 2006;Mechie et al., 1994). Some geochemical and geodynamical modeling studies propose two plumes exist beneath East Africa (e.g., Furman et al., 2006;Lin et al., 2005;Rogers et al., 2000). Seismological interpretations, however, remain debated with one or more plumes suggested to impinge beneath East Africa (e.g., Ritsema & Allen, 2003), most favoring the view that the entire region is underlain by the African Superplume-a broad, steep-sided mantle upwelling rising from the core-mantle boundary (Ni et al., 2002;Ritsema et al., 2011). The plume head associated with the initial flood magmatism, 45-30 Ma, would have long since dissipated, but multiple remnant plume stems have been proposed (e.g., Emry et al., 2019).
To address these debates, we have produced the first teleseismic P-and S-wave relative arrival-time tomography models of the Turkana Depression, using data from the Turkana Rift Arrays Investigating Lithospheric Structure (TRAILS) project, a multi-institutional, US-UK-Ethiopian-Kenyan geophysical experiment designed to fill a critical data gap along the EAR. Body-wave tomography can help constrain upper-mantle structure and offers greater lateral resolution than surface waves, which to-date form the only existing constraints of Turkana's uppermost-mantle structure (e.g., Benoit, Nyblade, & Pasyanos, 2006). At lithospheric depths, we explore whether low wavespeed zones associated with ongoing rifting are focused or broadly distributed. Below the extending African plate, our models illuminate the architecture of mantle upwellings and thus inform East Africa's single-versus-multiple plume debate.
East Africa's pre-rift lithosphere comprises numerous Precambrian terranes, including Archean cratons (e.g., the Tanzania craton) and flanking Proterozoic mobile belts. The main collisional phase, the East African Orogen (700-550 Ma), consisted of a period of protracted island-arc and microcontinent accretion (e.g., Fritz et al., 2013). Subsequent deformation and metasomatism in Neoproterozoic mobile belts has meant the EAR developed within these regions (e.g., Purcell, 2018). However, given their accretionary mode of formation, the mobile belts can themselves encompass ancient lithospheric fragments, detectable by seismic tomography (e.g., Emry et al., 2019).
Turkana's geology is largely concealed by considerable Phanerozoic sedimentary and volcanic strata, meaning its Proterozoic substrate has yet to be explored. Southern Ethiopia is thought to contain a mixture of juvenile Neoproterozoic lithosphere and Precambrian continental crust, a combination typical of an island-arc accretionary setting (Teklay et al., 1998). In the same region, Proterozoic ophiolites and sutures mark locations of past subduction (e.g., Fritz et al., 2013). Geochemical (e.g., Rooney, 2019) and petrological studies (e.g., Bedini et al., 1997) at Dilo and Mega (Figure 1b), lack any consensus on southern Ethiopia's lithospheric properties, except to note the impact of ancient depletion events, mostly of Proterozoic age. However, older depletion estimates (Bianchini et al., 2014) and zircon studies (Stern et al., 2012) suggest the presence of Archean lithosphere, termed the "Southern Ethiopian Shield," although its presence is unsubstantiated geophysically.
Cenozoic rifting in southern Turkana developed along the eastern edge of the Archean Tanzania craton 25 Ma, which is stronger than surrounding Pan-African lithosphere. The craton boundary is marked by  32 Ma carbonatites associated with initial heating and metasomatism at the steep craton edge, localizing magmatism and strain during the initial rifting stages (Muirhead et al., 2020). Since then, the locus of active faulting migrated eastwards, manifesting in 3-4 sub-parallel NNE-striking half-graben basins (Figure 1b) with diachronous development, contributing to the unusual broadening of the Turkana Depression (e.g., Ebinger et al., 2000;Morley et al., 1992). Magmatism has also migrated eastward to the Lake Turkana rift basin, and to N-S-trending, largely unfaulted, shield complexes further east of Lake Turkana that have erupted large volumes over the past 3 Ma (Figure 1b; Guth, 2016). Morley (2020) interprets dike swarms along the western margin of Lake Turkana as evidence for extension via magma intrusion since 3 Ma. These patterns suggest the linkage between the MER and Eastern Rift may have changed over time as strain migrated eastward, inconsistent with the simultaneous development of strain predicted in numerical models (Brune et al., 2017). Extension and magmatism propagated southward to northern Tanzania between 20 and 5 Ma (e.g., Baker, 1986;Mana et al., 2015). In northern Turkana, magmatism and basin formation initiated 18 Ma in two sub-parallel basins, and propagated northward to form the modern MER (e.g., Ebinger et al., 2000). Corti et al. (2019) suggest the MER is now propagating southward into the Ririba Rift east of Lake Turkana, based on structural and stratigraphic analyses. In contrast to the rest of Turkana, broad shield complexes are absent across the southern MER (Figure 1b). Ebinger and Sleep (1998) proposed that buoyant plume material impinged on the base of the Turkana lithosphere at 45 Ma, where thinned, base-of-the-lithosphere topography channeled plume material susceptible to melting to surrounding regions. Ethiopian Plateau uplift, often attributed to mantle plumes, is estimated to have commenced 20-30 Ma (Pik et al., 2003) with some studies postulating more rapid Late-Miocene uplift resulting from lithospheric foundering in response to extensive heating of the lithosphere (Furman et al., 2016;Gani et al., 2007). However, more recent drainage analysis argues against rapid uplift relating to lithospheric delamination processes (Sembroni et al., 2016). East African Plateau uplift is likely to have proceeded Ethiopian Plateau uplift, with some estimates around 14 Ma (e.g., Wichura et al., 2010). The Turkana Depression, although low-lying in comparison, was below sea-level prior to the Neogene (e.g., Wichura et al., 2015), evidenced by marine whale fossils and sedimentary strata (e.g., Brown & Mcdougall, 2011), and is thought to have risen by 600 m concurrently with East African Plateau uplift (Wichura et al., 2015).

Previous Geophysical Studies
Global tomographic studies have imaged a broad (500 km-wide) slow wavespeed structure originating at the core-mantle boundary beneath southern Africa, impinging on the lithosphere below East Africa (e.g., Li et al., 2008;Ritsema & Allen, 2003). Continent-scale ambient noise (Emry et al., 2019), surface-wave (e.g., Celli et al., 2020), body-wave tomography (e.g., Debayle et al., 2001;Hansen et al., 2012), and their joint inversion (Chang & van der Lee, 2011) recognize similar features, imaging Turkana as an anomalously slow wavespeed zone. Emry et al. (2019) image several, laterally distinct, slow wavespeed structures along the EAR, suggesting the presence of multiple plume stems linked to the broader African Superplume at depth. Other body-wave tomographic models find that the Superplume splits into two separate upper-mantle structures beneath both Afar and Tanzania (e.g., Chang & van der Lee, 2011;Hansen et al., 2012). More recently, however, Chang et al. (2020) and  find evidence for two separate whole-mantle plumes: a near-vertical plume beneath the Afar triple junction as well as the African Superplume. The morphology of the latter, as determined from gravity data (Simmons et al., 2007) and seismic waveform modeling (e.g., Ni et al., 2002), is characterized best as a sharp-sided thermochemical plume, whereas the plume below Afar is likely a purely thermal upwelling (e.g., Boyce & Cottar, 2021).
Regional relative arrival-time tomographic studies provide strong evidence for continuation of the Superplume into the Ethiopian upper mantle (e.g., Bastow et al., 2005Bastow et al., , 2008Benoit, Nyblade, & VanDecar, 2006;Civiero et al., 2015). Beneath the East African Plateau, tomography studies reveal a slow wavespeed anomaly in the upper mantle, interpreted as ascending plume material (e.g., Mulibo & Nyblade, 2013;Weeraratne et al., 2003). Tomography, in conjunction with receiver function studies of the mantle transition zone (e.g., Benoit, Nyblade, Owens, & Stuart, 2006;Nyblade et al., 2000;Thompson et al., 2015) that observe near-normal mantle transition zone thicknesses below the MER but thinned in southern Ethiopia, suggest the plume is likely to cross the mantle transition zone below the Turkana Depression, which has been a large data gap prior to this study.
At shallower depths ( 100 km), a focused (150 km-wide) low wavespeed rift-zone is imaged below the MER (Bastow et al., 2005(Bastow et al., , 2008 and Eastern Rift (Park & Nyblade, 2006;Tiberi et al., 2019). The Kenya Rift International Seismic Project (e.g., Prodehl et al., 1994) conducted wide-angle reflection studies to determine variations in crustal thickness and uppermost-mantle structure along the axis of the Eastern Rift, reaching as far north as Lake Turkana (e.g., Achauer & The KRISP Teleseismic Working Group, 1994;Mechie et al., 1994). Crustal thicknesses of 35 km and 20 km are revealed beneath the East African Plateau and Lake Turkana respectively, the latter ascribed to widespread Mesozoic rifting (Mechie et al., 1997). These authors find comparable results to the Rayleigh wave dispersion model of Benoit, Nyblade, and Pasyanos (2006), where shallow slow wavespeeds, spanning the whole Depression, suggest underlying lithospheric thinning with an average crustal thickness of 255 km. Overlapping phases of prolonged rifting have rendered Turkana the area of thinnest continental crust along the EAR excluding Afar, with crustal stretching factors of 1.2-1.6 (Ebinger & Ibrahim, 1994;Hendrie et al., 1994). Knappe et al. (2020) measure higher present-day strain rates in the southern Turkana Depression than the northern MER, but geodetic data are too sparse to evaluate the current localization of strain across the 300-km-wide broadly rifted zone in southern Ethiopia.

The TRAILS Seismograph Network and Teleseismic Dataset
Seismograph deployments in the Turkana Depression and adjoining areas on the MER and Eastern Rift were completely lacking compared to the rest of the EAR until 2019, when the NSF-NERC funded, US-UK-Ethiopian-Kenyan, TRAILS project began. The TRAILS broadband seismograph network ( Figure 2a) comprised 34 Güralp CMG-ESP and CMG-3T instruments, capable of responding to periods reaching 60 and 120 s, respectively (see Table S1 in Supporting Information for further station details). Additional data were sourced from permanent GEOFON seismic station LODK (Figure 2a during the period January 2019 to September 2020, were inspected visually for a good signal-to-noise ratio on at least 4 stations. Station-earthquake pairs with    30 were excluded to avoid non-unique raypaths caused by mantle-transition-zone-associated triplications. Ultimately, 607 earthquakes were analyzed for direct P-phases, together with 8 core-diffracted (Pdiff) and 70 core (PKP, PKIKP, and PKiKP) phases. Similarly, inspection yielded 205 direct S-and 105 SKS-phase earthquakes (Figures 2b and 2c). To provide greater resolution and improve the uniformity of the data coverage with respect to backazimuth and epicentral distance, further visual inspection of 16,738 lower magnitude earthquakes (4.0     m b 5.0) was conducted. This yielded an additional 294 earthquakes for P-wave analysis and 67 for S-wave analysis. Earthquake coverage is best from the east and south east (western Pacific subduction zones: 60-130; Figures   (b and c) The backazimuthal and distance distribution of earthquakes used in this study between January 2019 and September 2020. (b) m b  4.0 direct P-wave (circles) and core phase (crosses and red circles) event epicenters. (c) m b  4.0 direct S-wave (circles) and SKS (crosses and red circles) event epicenters. The concentric circles mark 30 intervals in epicentral distance from the center of the network at 4N, 37E (black star). Note that although some earthquakes are located   30 from the study area on the figure, no individual raypaths with   30 were included to avoid triplications associated with the mantle transition zone.

Method of Relative Arrival-Time Determination
Filtered waveforms, using a zero-phase, two-pole, Butterworth bandpass filter with corner frequencies of 0.4 and 2 Hz, and 0.04 and 0.2 Hz for P-and S-waves, respectively, were manually picked on the first identifiable coherent peak or trough of body-wave energy across the network. These bandwidths are comparable to previous studies on the EAR (e.g., Bastow et al., 2008), but with a higher S-wave low-pass cutoff than is used in some studies (0.1 Hz: Bastow et al., 2005;Benoit, Nyblade, & VanDecar, 2006). Low-pass filter frequency cut-offs were kept as high as possible to allow for the inclusion of high frequency signals, satisfying the infinite frequency approximation assumed in ray theory which we adopt in our inversions. Direct Pand S-phases were picked on the vertical and tangential component seismograms, respectively; SKS-phases were picked on the radial. Subsequent improvement of manual trace alignment was employed through Multi-Channel Cross-Correlation (MCCC; VanDecar & Crosson, 1990) allowing for accurate computation of relative arrival-time residuals.
The MCCC procedure windows manual phase picks (3 s for P-waves; 12 s for S-waves) to avoid interference from secondary phase arrivals, and searches for the cross-correlation maximum position,  max ij , between each pair of traces for a single earthquake. These are subsequently used to determine the most appropriate cross-correlation derived relative delay-time,  ij t , for each station pair for each earthquake: where P i t and P j t are preliminary-picked-arrival estimates for the ith and jth traces, respectively. Rays with a mean cross-correlation coefficient 0.85 were rejected. A least-squares minimization procedure is required to optimize relative arrival-time measurements for each station by forcing the mean over all stations to be zero. The standard error of each arrival is also quantified through the MCCC method: relative arrival-times for the P-and S-wave dataset have a mean standard deviation of 0.01 and 0.06 s, respectively. Given the trace sampling rate of 0.01 s these errors were considered optimistic, agreeing with comparable studies (e.g., Tilmann et al., 2001).
Relative arrival-time residuals ( rel t ) were then computed for each station: where i t is the relative arrival-time for each station i, e i t is the expected IASP91 (Kennett & Engdahl, 1991) arrival-time for station i, and e t is the mean IASP91 predicted arrival-time for the specific earthquake across the network. The final datasets comprise 14,420 P-wave and 4,813 S-wave relative arrival-times.

Analysis of Relative Arrival-Time Residuals
Relative arrival-time methods document seismic anomalies with respect to the regional, not global, mean (e.g., Bastow, 2012). To help us understand better the pedestal in which our relative arrival-times sit on a global scale, we computed absolute arrival-times for each station using the Absolute Arrival-Time Recovery Method (AARM; see Section S3 and Figure S1) of Boyce et al. (2017). In our P-wave dataset, a near-zero relative arrival-time residual exists at station QORK, where an absolute arrival-time residual of  t p abs = 1.11 s is observed, corresponding to an average of  V p  −1.4% slow upper mantle (0-670 km depth) relative to IASP91. Similarly, for S-waves at station BASK,  t s abs = 4.58 s, implying a V s  −3.3% slow upper mantle. Therefore, given our range of relative arrival-time residuals (Figures 3a and 3b), both positive and negative arrivals in this study are, on average, late compared to the global mean.
Both P-and S-wave mean relative arrival-times (Figures 3a and 3b) reveal a trend of late arrivals around Lake Turkana ( t p  1 s;  t s  3 s), and early arrivals toward the eastern edge of the network ( t p  −1 s;  t s  −3 s); a trend also evident in displaced residual histogram peaks when grouped by tectonic domain (Figures 3c and 3e). These data have been corrected for topography using correction velocities of 3 km/s and 1.7 km/s for P-and S-waves, respectively. Most stations across the network display a consistent residual variation with backazimuth and epicentral distance (Figures 3d, 3f and 4).
Stations in eastern Turkana consistently display the earliest arrival-time residuals (e.g., MOYE:  t p = −0.78  0.34 s;  t s = −3.31  0.73 s) with relatively early arrivals from all backazimuthal and epicentral distances, except for some rays arriving from the west (225-270 backazimuth; Figure 4). Eastern Turkana stations also display large peak-to-peak residual variations ( t p  2 s;  t s  5 s). The eastern part of the TRAILS network may thus encompass the eastern edge of the African Superplume. Stations BUBE ( p t = 0.52  0.21 s), TBIK ( p t = 0.88  0.12 s), and KERK ( p t = 0.5  0.18 s), located along the length of Lake Turkana, have small peak-to-peak variation ( t p  1 s;  t s  4 s) and the largest mean arrival-time residuals, mostly originating from short, eastern, epicentral distance earthquakes. These stations are located in 3 km-thick sedimentary basins (Sippel et al., 2017). Residual trends are more pronounced in the S-wave dataset (e.g., BUBE:  s t = 1.35  1.12 s), highlighting S-wave sensitivity to partial melt. Mean relative arrival-times drop in a NW-SE-trending band across southern Ethiopia, possibly suggesting a break in focused, rift-related low wavespeeds. Stations GAJO ( p t = -0.02  0.12 s) and QORK ( p t = 0.0  0.16 s), in the central Turkana Depression, have mean relative arrival-time residuals close to zero (i.e., near the regional background mean) and the smallest peak-to-peak variations ( t p  0.5 s;  t s  3 s).

Model Parameterization and Inversion Procedure
We adopt the regularized, least-squares inversion method of VanDecar et al. (1995) to invert relative arrival-time residuals for mantle velocity perturbations. Source terms and station statics were included to account for distant heterogeneity/event mislocations, and unresolvable near-surface structure, respectively. Similar inversion methods have been used by other regional studies along the EAR (e.g., Bastow et al., 2008;Mulibo & Nyblade, 2013). Model slowness is parameterized using B-splines under tension (Cline, 1981) over a dense grid of 64,800 knots, covering 27 knots between 0-1,000 km in depth, 50 knots between 26-47  E in longitude, and 48 knots between 5  S-14  N in latitude. The central, best resolved, portion of the model is sampled at 25 km intervals at  400 km depth, and 0.1  in latitude/longitude. We extend parameterization, at coarser spacing, beyond the area of interest to create zones where spurious structure can map without being incorporated into the interpreted parts of the model.
Our inversion method utilizes ray theory, which is an infinite-frequency approximation where arrival-times are only influenced by structures along an infinitesimally narrow region. Following Montelli et al. (2004), Fresnel zones in our study are of the order L for a wave of wavelength, , and ray length, L, where  is determined from the average frequency content of our data (1.2 Hz for P-waves and 0.12 Hz for S-waves) and L represents ray length which in this case we assume to be approximately equal to the depth. In the uppermost, best-resolved, parts of our model, we can thus justify resolution of features no smaller than 20 km and 50 km in the P-and S-wave models, respectively. KOUNOUDIS ET AL.  Because our inverse problem is under-determined (more unknowns than observations), we regularize through the minimization of a seven-point finite element approximation to the Laplacian operator to penalize model roughness and yield a smooth model that fits the data. A preferred model was chosen by investigating the trade-off between root-mean-square (RMS) residual reduction (data fit) and model roughness (see Figures S2 and S3 in Supporting Information), avoiding models achieving greater residual reduction than justifiable from MCCC-defined data noise levels which we consider, like others, optimistic (e.g., Villemaire et al., 2012). The chosen models fit 92.27% and 89.07% of the RMS estimated relative-delay times for the P-and S-wave models, respectively (see Tables S2 and S3 in Supporting Information). When subtracting the station static terms from the delay times, the RMS relative arrival-time residuals decrease from 0.33 to 0.25 s for P-waves and from 1.29 to 1.12 s for S-waves; the resulting residuals reflect more accurately the proportion of delay time anomalies incorporated into the region of the model where we make our interpretations.

Resolution
To assess resolution, a standard checkerboard analysis was conducted ( Figure 5 and Figure S4), using alternating V p s ,     5%, 50 km diameter spherical anomalies, defined by Gaussian functions across their diameter: structures 50 km are not warranted since they lie below the S-wave model Fresnel zone. The spherical anomalies were placed at 200 km depth intervals from 100 to 500 km depth, every 1 in latitude/longitude, then inverted using identical model parameterization and regularization to the observed dataset inversion. Travel-times are calculated assuming ray paths through 1D velocity model IASP91. Gaussian residual arrival-time error components with standard deviations 0.03 and 0.1 s were added to theoretical arrival-times for P-and S-waves, respectively. The retrieved spheres are laterally distinct, with an amplitude recovery of 50% at 300 km depth. However, anomalies experience increasing lateral smearing and reduced amplitude recovery (40%) with depth ( Figure 5). To assess vertical smearing further, we place a 150 km diameter sphere at 4N, 36E at 600 km depth, to simulate a large low-velocity mantle-transition-zone anomaly. Some vertical smearing is observed throughout the model, increasing toward the model edges ( Figure 5).
Our dataset's ability to retrieve geodynamically plausible structures is tested via several realistic synthetic tests. Models were created to simulate localized rifting, multiple segmented upwellings versus a laterally continuous superplume, and assess smearing between a mantle plume extending up to the mantle transition zone and shallow lateral flow ( Figure 6). Gaussian residual time errors were added and the data inverted as with the checkerboard tests. Amplitude recovery in the P-wave model is  60% (Figure 6), increasing with depth in central Turkana and decreasing toward the west. Long-wavelength anomalies are inevitably retrieved more easily than steep-gradient, short length-scale structure when regularizing via smoothing: hence the greater amplitude recovery compared to checkerboard tests. Lateral smearing is minimal in both models, however, S-wave model amplitude recovery (50%; Figure S5) is lower and smearing more pronounced, than the P-wave model. The multiple plume model smears upward 100 km (Figure 6; C-C'). Lateral smearing near the mantle transition zone is minimal. Vertical resolution is limited to  150-200 km, since features with smaller separation smear and merge, as shown in the synthetic combining shallow low velocities with a deep mantle-transition-zone plume anomaly (Figure 6; A-A'). Low wavespeeds along the length of Lake Turkana would be resolved at lithospheric depths, but smear downwards somewhat into deeper features that smear upwards.

Tomographic Results
The uppermost mantle beneath Turkana is marked by pronounced low wavespeed anomalies ( p V = −1.5%,  s V = −2.5%; Figures 7 and 8), widely distributed across the region that persist to mantle-transition-zone depths. At  150 km depth, low wavespeeds coincide with recent Pleistocene and Holocene volcanic centers which are themselves broadly distributed throughout Turkana (Figure 1b). Only south of Lake Turkana, along the Suguta Valley do we see evidence for low wavespeeds confined to a relatively narrow, 100-kmwide zone at shallow depths ( p V = −1.5%;  s V = −3%; Figures 7 and 8; B-B'). The north-south continuity of shallow low wavespeeds is broken at 5N by an elongated, narrow, high wavespeed anomaly ( p V = 1.5%), traversing NW-SE in southern Ethiopia to  200 km depth (Figure 7; 100 km; C-C', E-E'). This feature is less    pronounced in the S-wave model (Figure 8; 100 km), perhaps due to the lower resolving power of the noisier, longer wavelength S-wave dataset. High wavespeed anomalies are also present in the Lotikipi basin west of Lake Turkana to 100 km depth (Figures 7 and 8; 75-100 km). East of Lake Turkana, a pronounced, laterally robust, N-S-trending, high wavespeed ( p V = 1.5%) structure persists sub-vertically to 600 km depth, forming an eastern boundary to the low wavespeeds (Figure 7; A-A'). At 300 km depth, the entire region is underlain by broader low wavespeeds ( p V = −1.5%;  s V = −3%; Figure 7), particularly in the S-wave model (Figure 8; 500 km). Laterally distinct low wavespeed heterogeneity within the mantle transition zone exists in the P-wave model (Figure 7; 500-600 km), but resolution of such features is limited outside central Turkana ( Figure 5).

Causes of Seismic Heterogeneity
When interpreting relative arrival-time tomographic models, the  , p s V = 0% contour seldom matches the global zero mean (e.g., Bastow, 2012). Even peak-to-peak amplitude variations will not necessarily document how anomalous the mantle in a region truly is: if phase arrivals across the network are all consistently late or early, that information is lost completely during relative arrival-time calculation (Equation 2). Taking a 0-500 km depth average through global P-and S-wave models LLNL-G3Dv3 (Simmons et al., 2012) and SL2013sv (Schaeffer & Lebedev, 2013) below stations where  , t p s = 0 (stations QORK and BASK for P-and S-waves, respectively),  , p s V = 0% in Figures 7 and 8 could be considered to differ from the global mean by  p V = −1.24% and  s V = −3.13%. These estimates corroborate our absolute arrival-time observations of  t p abs = 1.11 s, and  t s abs = 4.58 s, which imply  V p  −1.4% and V s  −3.3% for the Turkana upper mantle below QORK and BASK, respectively. These values contrast the melt-rich MER to the north of our study area, where  p V is considered to be 5% slow compared to the global mean (e.g., Bastow et al., 2005;; V s  −11% (Bastow et al., 2010;Gallacher et al., 2016).
Various factors can affect uppermost-mantle wavespeeds, but temperature exerts greater influence than composition: for compositions lacking strongly depleted Mg-rich hazburgites, velocity variations are usually 1% (e.g., Cammarano et al., 2003;Goes et al., 2000). Anisotropy can also affect wavespeed variations but is unlikely to be a major contributor to the observed anomalies, as evidenced by our good backazimuthal earthquake coverage (Figure 2) and the largely consistent P-and S-wave models (Figures 7 and 8). We thus first explore the expected mantle thermal structure when assuming our seismic velocity anomalies are purely due to temperature variations.

Elevated Mantle Temperature and Melt
Smooth, isomorphic temperature derivatives for a pyrolitic composition along a 1300C adiabat, d s V /dT conversions taken from Styles et al. (2011) and d p V /dT based on mineral parameters from database stx08 (Xu et al., 2008) and composite attenuation model Qg (van Wijk et al., 2008), predict a 100C thermal anomaly associated with decreases in V p  1.5% and V s  2.5% at 100 km depth, to account for anelastic and anharmonic effects on seismic wavespeeds. Our observed peak-to-peak variations of V p  3% and V s  5.5%, at 100 km depth (Figures 7 and 8), yield lateral temperature variations of   200C and  220C for the P-and S-wave models, respectively, if cross-network heterogeneities are attributed to temperature alone. These are likely lower bound temperature estimates since anomaly amplitudes are invariably underestimated in regularized, under-determined tomographic inversions. Additionally, since the relative wavespeed anomalies we retrieve are on a pedestal that is slow ( V p  −1.4%;  V s  −3.3%), an additional thermal component of   130C would be required to explain Turkana's mantle wavespeeds. Such high-temperature estimates are petrologically implausible, meaning mantle melt is also required to explain the observations: lavas erupted in Turkana since 10 Ma reveal mantle potential temperatures of 1,450-1500C (Furman et al., 2006;Rooney et al., 2012); 100-150C above ambient mantle.
Assessing arrival-times rather than velocity perturbations for temperature estimation, circumnavigates amplitude recovery issues associated with smoothed, under-determined tomographic velocity models (e.g., differing resolutions, ray-paths, and regularization levels between the P-and S-wave models), and artifacts associated with the inversion procedure. Several studies have cited the ratio of S-and P-wave arrival-time ) for common earthquake-station pairs as a diagnostic tool for causes of seismic anomalies (e.g., Bastow et al., 2005;Bolton & Masters, 2001;Civiero et al., 2016). A . for absolute arrival-time residuals. The purple line is the slope in the partial melt affected Main Ethiopian Rift (slope of 10.00; Bastow et al., 2005). The red line is the slope in eastern Anatolia (slope of 4.01; Kounoudis et al., 2020), where anomalies are primarily due to temperature and small volumes of mantle melt.
are path-integrated values and S-waves are generally longer wavelength than P-waves, these ratios are likely to be underestimates (e.g., Schmandt & Humphreys, 2010).

High Velocity Anomaly in Southern Ethiopia: By-Product of Cenozoic Magmatism or Proterozoic-Age (Pan-African) Lithospheric Heterogeneity?
Perhaps the most distinctive feature at 0-200 km depth in the P-wave model (Figures 7 and 10a), is a NW-SE-trending high wavespeed band (V p = 1.5%) in the broadly rifted zone of southern Ethiopia, just south of where the MER broadens to three sub-parallel rift basins: Omo, Weyto-Chew Bahir, Segen (Figures 1b,  10a and 10b). Synthetic tests ( Figure 6) show that high amplitude, 40-km-wide structures can be recovered from the P-wave inversions. The feature has sharp, near-vertical, anomaly boundaries (Figure 7; C-C'), potentially a result of juxtaposed hot/partially molten mantle and a relatively thick, cold, and compositionally distinct lithospheric fragment; purely thermal anomalies are expected to be more laterally diffuse. The more subtle increase in  s V (Figure 8) observed coincident with the high P-wave band may imply a lithospheric composition to which S-waves have a lower sensitivity than P-waves, however, we cannot preclude the possibility that the model differences are due, at least in part, to the larger Fresnel zone (50 km) and reduced spatial resolution of the S-wave dataset.

Archean Lithosphere
Although most geochronological studies (e.g., Reisberg et al., 2004;Yeshanew et al., 2017) suggest the southern Ethiopian mantle lithosphere is Neoproterozoic in age, Stern et al. (2012) and Bianchini et al. (2014) attribute the presence of abundant 2.5 Ga zircons to Archean lithosphere at depth in the southern Ethiopian Shield. Xenoliths around Mega (Figure 10b) also preserve evidence of ancient Precambrian depletion events, unlike elsewhere in the northern sector of the EAR (e.g., Beccaluva et al., 2011). Iron-depleted, Archean lithosphere is characterized by fast wavespeeds (Griffin et al., 2009), affecting  s V more strongly than  p V . This contradicts the observation that the high wavespeed band is illuminated best in our P-wave model (Figures 7 and 8). With the caveat that the lower resolving capabilities and larger Fresnel zone of the S-wave dataset may at least partly explain the model differences, we suggest that the P-and S-wave amplitude recovery, when reviewed in light of the lack of contiguity with other Archean blocks in East Africa, render the Archean lithosphere hypothesis unlikely.

Cenozoic Processes
Alternatively, the high wavespeed band may be a by-product of Cenozoic magmatism. The earliest phase of basaltic magmatism in East Africa (40-45 Ma; e.g., Davidson & Rex, 1980) occurred in SW Ethiopia, potentially depleting the lithosphere of its most easily fusible elements. Melt depletion is generally characterized by increased wavespeeds (e.g., Lee, 2003;Matsukage et al., 2005), with most studies favoring greater increases in V p than V s (Saltzer & Humphreys, 1997;Schutt & Lesher, 2010), as per our results (Figures 7  and 8). However, melt depletion would have to occur through fractional melting to impact seismic wavespeeds discernibly (e.g., Afonso et al., 2008;Schutt & Lesher, 2006). Southern Ethiopia is considered to have experienced relatively low degrees of lithospheric mantle melting (Alemayehu, Zhang, & Sakyi, 2017;George & Rogers, 2002;Rooney, 2010). Petrological studies at Dilo/Mega (Figure 10b; e.g., Alemayehu, Zhang, & Seitz, 2017;Meshesha et al., 2011) instead attest to widespread Cenozoic metasomatism from intruded fertile plume material (e.g., Nelson et al., 2012); these would impart seismically slow phases in the lithospheric mantle (Eeken et al., 2018) whilst also increasing its temperature. Therefore, the style, location, and spatial extent of magmatism (in 10 Ma at least 30,000 3 km of magma was widely dispersed across SW Ethiopia; Figure 10b; Ebinger et al., 1993) does not lend itself to the generation of the linear melt depletion zone required to generate our arcuate, high velocity anomaly.
evidence for downward-thrusted ancient oceanic lithosphere. Teklay et al. (1998) also established, through zircon dating and isotope geochemistry, that southern Ethiopia comprises a mixture of juvenile Neoproterozoic lithosphere and small volumes of Precambrian continental crust, a combination typical of a Pan-African island-arc setting.
Quaternary lavas in Turkana show remarkable geochemical similarity (e.g., Furman et al., 2006;Rooney, 2019), except those in southern Ethiopia: Rooney (2019) notes that around Mega (Figure 10b), volcanism is unusually mafic. In geochemically distinct southern Ethiopia, xenolith studies (e.g., Bedini et al., 1997;Bianchini et al., 2014) lack any consensus concerning lithospheric composition. Mega lies at the edge of our high wavespeed band, hinting at a possible transition between two lithospheric domains as a cause for our observations: high olivine forsterite (Fo) content in Mega xenoliths is consistent with oceanic peridotite (Alemayehu, Zhang, & Sakyi, 2017), signifying compositional heterogeneity as a cause for the higher wavespeeds. An argument for continental lithosphere has also been made: the 45-35 Ma southern Ethiopian basalts have Na and Fe contents consistent with melting beneath thick continental lithosphere (e.g., George & Rogers, 2002). Either way, the East African Orogeny was certainly a period of protracted island-arc and microcontinent accretion involving juvenile Neoproterozoic oceanic crust that formed in and adjacent to the Mozambique ocean (e.g., Fritz et al., 2013). Our linear, elongate, high wavespeed band may thus be attributed to a fragment of oceanic or continental lithosphere in the presence of a high-grade ocean arc accretion zone. A relict Proterozoic subduction zone trapped in continental lithosphere potentially forms a rheological feature, strongly affecting the localization and development of subsequent continental tectonics. Such a scenario is not without modern-day analog. Oceanic lithosphere trapped in the uppermost mantle has been proposed in the NW United States: the 55 Ma accretion of the Siletzia microplate during the Cordilleran Orogeny is imaged seismically as a curtain of high wavespeeds to 200 km depth (Schmandt & Humphreys, 2011). Siletzia appears to be a competent crustal block with a diffuse layer of seismicity close to the Moho (Merrill et al., 2020). Zhang et al. (2009) also imaged a 100-km-wide high wavespeed band beneath the south-central part of the Gulf of California that is attributed to a 12 Ma slab fragment, trapped at lithospheric depths.

Implications for Mesozoic and Cenozoic Strain Localization
The lack of Quaternary eruptive centers and comparatively thin sequence of Eocene-Oligocene flood basalts coinciding with the Southern Ethiopia high wavespeed band (Figure 10b; Rooney, 2017) suggest the presence of relatively cold, refractory Proterozoic lithosphere. Intriguingly, the band delimits the northernmost extent of Mesozoic extension in Turkana, and marks an offset in the surface expression of the Mesozoic rift system (Figure 10a). Excluding the zone of continuation of the MER at the center of the high wavespeed band, a marked increase in topography (1.5 km) occurs south-to-north ( Figure 10c; A-F): abrupt in the west; more subtle in the east (Figure 10c; A-F). These topographic observations support a pre-Mesozoic age for the high wavespeed band, and potentially deem the band a refractory feature governing the northern limit of Mesozoic extension.
A zone of strong lithosphere should exert some control on EAR development too. Intriguingly, Cenozoic rifting lacks focus within and immediately to the north of the high wavespeed zone, where present-day seismicity (Figure 10a) is also diffuse, in contrast to the relatively narrow (80 km-wide) expression of the MER north of 6 5 . N.
Previous analog experiments and numerical models that have attempted to explain the linkage between the MER and Eastern Rift (Brune et al., 2017;Corti et al., 2019) assert the importance of previously thinned lithosphere below Turkana. However, the seismically fast band of refractory lithosphere illuminated in southern Ethiopia by our tomographic models, may have influenced both Mesozoic and present-day strain localization, likely including the complex, broadly rifted, transfer zone between the MER and Eastern Rift within Turkana. The initial conditions in the analog experiments and numerical models therefore require some revision: specifically, a strong zone in southern Ethiopia, not just a simple previously thinned lithosphere in the Turkana Depression between two plateaus, and a diachronous migration and propagation of several 100-km-wide rift zones.

The Central Turkana Uppermost Mantle: Ponded Asthenosphere or Melt-Infiltrated Lithosphere?
Below the 80-km-wide central and northern MER to the north of the TRAILS network, the uppermost mantle is characterized by large low wavespeed anomalies at 200 km depth (Bastow et al., 2005(Bastow et al., , 2008. Similarly focused mantle anomalies also characterize the Eastern Rift (Achauer & The KRISP Teleseismic Working Group, 1994;Mulibo & Nyblade, 2013;Park & Nyblade, 2006;Tiberi et al., 2019). Our results corroborate these observations in the south of our study area, beneath the narrow Suguta Valley (Figures 7  and 8; 75-100 km; B-B'). In the absence of high-resolution lithospheric thickness measurements, the low wavespeeds could signify concentrated melt intrusions in still-thick lithosphere or a lithospheric thin zone where laterally flowing plume material is ponding. However, the presence of extensive, Pliocene-Recent dikes (Morley, 2020) and high present-day strain rates (Knappe et al., 2020), mostly focused at the southern end of Lake Turkana, agree with the gradual focusing of magma-assisted rifting to narrower zones.
Absent from our models across most of Turkana is a single obvious N-S-trending, localized, low-velocity structure at <150 km depth (Figures 7 and 8). Instead, the focused low wavespeed zone at 100 km depth in Suguta bifurcates in central Turkana and merges into a laterally continuous feature in southern Ethiopia ( p V = −1.5%;  s V = −3%; Figure 11). The absence of a focused, narrow band of low velocity material is consistent with Turkana's missing rift-valley morphology, widespread faulting and sedimentary basins (Figure 11;Ebinger et al., 2000). The low wavespeeds instead mirror volcanic centers and patterns of seismicity, confined to the branches in central Turkana and become more broadly distributed south and north of the high wavespeed band in southern Ethiopia (Figure 11). Our S-to P-wave arrival-time ratio analysis (Figure 9) does not make a strong case for voluminous mantle melt. Extension-related decompression melting is also likely negligible due to the ultraslow extension rates that dominate Turkana (4 mm/yr; e.g., Birhanu et al., 2016;Knappe et al., 2020). Flanking relatively high wavespeed anomalies at 100 km depth, mostly coinciding with regions of failed Mesozoic and Paleogene rifting in the Lotikipi and Anza basins ( p V = 1%;  s V = 1%; Figures 7 and 8), perhaps signify areas of relatively colder and re-thickened lithosphere (since failed Mesozoic and Paleogene rifting) that GPS data reveal are not actively straining (Knappe et al., 2020). Thus, because we cannot completely preclude the possibility of melt at mantle lithospheric depths, whether KOUNOUDIS ET AL.  (Global Volcanism Program, 2013). Gray circles are earthquake epicenters occurring from January 2019 to September 2020 (Musila et al., 2020), increasing in depth with darker shades. The shaded blue region denotes the Cretaceous Anza rift, and the shaded white region shows the regions of Paleogene extension (e.g., Boone, 2018;Ebinger & Ibrahim, 1994).
Turkana's low wavespeed branches at 100 km depth are ponded asthenospheric material beneath broad zones of variably pre-thinned lithosphere, or heavily melt-intruded mantle lithosphere is equivocal. These low mantle wavespeeds are confined to the Eastern Rift, showing no evidence for connection to the Western Rift at the northern end of the East African Plateau.

The Turkana Sub-Lithospheric Mantle and Transition Zone
Our tomographic models reveal low velocity anomalies ( p V = −1.5%;  s V = −3%) that persist with depth, at least to the mantle transition zone (Figures 7 and 8; A-A', D-D'), however we cannot preclude the possibility of separate sub-lithospheric and lithospheric low wavespeed anomalies that are merged by vertical smearing (Figures 5 and 6). Absolute arrival-times of P-and S-waves recorded by the TRAILS network also attest to significant, deep-seated, low wavespeed mantle structure, particularly from southerly backazimuths, where  t s abs  14-17 s (Figure 9). Collectively, our observations support the view that African Superplume material dominates sub-lithospheric upper mantle heterogeneity below Turkana.
Despite a lack of coverage below Turkana, several studies (e.g., Simmons & Grand, 2002;Simmons et al., 2007) show thermal and compositional complexity in the East African mantle transition zone, consistent with the influence of mantle plumes. Temperature impacts  p V and  s V similarly, however, our S-wave model depicts more concentrated low wavespeeds at mantle-transition-zone-depths than the P-wave model (Figures 7 and 8; 500-600 km; A-A'). This discrepancy could be the result of compositional heterogeneity or a hydrated mantle transition zone (e.g., Cornwell et al., 2011;Meier et al., 2009;Thompson et al., 2015), consistent with the near ubiquitous view that the African Superplume is compositionally, not just thermally anomalous Ni et al., 2002;Simmons & Grand, 2002;Simmons et al., 2007).
A fundamental question concerning Turkana is whether or not a lack of dynamic upwelling or lithospheric structure are the main causes for its low-lying nature compared to the surrounding plateaus. Our observations of deep-seated, low velocity material provide no support for the hypothesis that Turkana's low elevation is due to lacking dynamic support. With the caveat that high-resolution crustal and lithospheric thickness measurements are lacking in Turkana, an isostatic response contribution from overlapping Mesozoic and Cenozoic phases of rifting (Purcell, 2018;Reeves et al., 1987), significantly thinning the lithosphere (Benoit, Nyblade, & Pasyanos, 2006;Mechie et al., 1994), thus instead likely govern Turkana's low elevations relative to the two plateaus it separates.

Conclusions
Using data from the recently deployed TRAILS seismograph network of 34 broadband stations, we present the first P-and S-wave tomographic models of upper-mantle seismic structure beneath the Turkana Depression.
Relative arrival-time residuals and velocity-constrained temperatures at 100 km depth exceed petrological mantle potential temperature estimates, suggesting the presence of mantle melt, albeit in smaller volumes than below the MER to the north. At 200 km depth, a NW-SE-trending, 50 km-wide, high wavespeed anomaly ( p V = 1.5%) is interpreted as relatively refractory Proterozoic lithosphere. The coincidence of this high velocity band with the northern limit of the Turkana Depression, and the southerly extent of MER-related fault scarps, implies the band has exerted some control on the spatial extent of Mesozoic and Cenozoic rifting. At lithospheric mantle depths (100 km), a single localized low wavespeed zone ( p V = −1.5%;  s V = −2%) is only present in the southernmost part of our study area where a focused zone of lithospheric extension in the Eastern Rift is illuminated. North of this, within central Turkana, the low wavespeed zone bifurcates to either lithospheric thin-spots that host ponded asthenospheric material or regions of melt-intruded mantle lithosphere, that merge into a broadly distributed low wavespeed zone in southern Ethiopia. Below mantle lithospheric depths, low velocity anomalies ( p V = −1.5%;  s V = −3%) provide evidence that dynamic support is continuous below East Africa, not just below the uplifted Ethiopian and East African plateaus. Figure 12 illustrates our main conclusions.