The Impact of Complex Volcanic Plumbing on the Nature of Seismicity in the Developing Magmatic Natron Rift, Tanzania

Constraining the architecture of complex 3D volcanic plumbing systems within active rifts, and their impact on rift processes, is critical for examining the interplay between faulting, magmatism and magmatic fluids in developing rift segments. The Natron basin of the East African Rift System provides an ideal location to study these processes, owing to its recent magmatic-tectonic activity and ongoing active carbonatite volcanism at Oldoinyo Lengai. Here, we report seismicity and fault plane solutions from a 10 month-long temporary seismic network spanning Oldoinyo Lengai, Naibor Soito volcanic field and Gelai volcano. We locate 6,827 earthquakes with ML −0.85 to 3.6, which are related to previous and ongoing magmatic and volcanic activity in the region, as well as regional tectonic extension. We observe seismicity down to ∼17 km depth north and south of Oldoinyo Lengai and shallow seismicity (3–10 km) beneath Gelai, including two swarms. The deepest seismicity (∼down to 20 km) occurs above a previously imaged magma body below Naibor Soito. These seismicity patterns reveal a detailed image of a complex volcanic plumbing system, supporting potential lateral and vertical connections between shallow- and deep-seated magmas, where fluid and melt transport to the surface is facilitated by intrusion of dikes and sills. Focal mechanisms vary spatially. T-axis trends reveal dominantly WNW-ESE extension near Gelai, while strike-slip mechanisms and a radial trend in P-axes are observed in the vicinity of Oldoinyo Lengai. These data support local variations in the state of stress, resulting from a combination of volcanic edifice loading and magma-driven stress changes imposed on a regional extensional stress field. Our results indicate that the southern Natron basin is a segmented rift system, in which fluids preferentially percolate vertically and laterally in a region where strain transfers from a border fault to a developing magmatic rift segment.


INTRODUCTION
Continental rifting is a geodynamic process in plate tectonics. However, the forces enabling the onset of rifting of comparatively strong, cold, thick continental lithosphere and the relative importance of magmatic and extensional processes remain poorly understood (e.g., Bialas et al., 2010;Brune et al., 2017).
A key question is how strain accommodation is partitioned between faulting and magmatism (e.g., Ebinger and Casey, 2001;Thybo and Nielsen, 2009;Marzen et al., 2020). Tectonic thinning is commonly accompanied by decompression melting and the rise of volatiles and magma through the lithosphere (White and McKenzie, 1989;Weinlich et al., 1999;Rooney, 2010;Lee et al., 2016). Therefore, intrusion of magma and volatile release are thought to play a pivotal, yet interacting, role in the continental rifting process (Lindenfeld et al., 2012a;Wright et al., 2012;Muirhead et al., 2016;Weinstein et al., 2017).
The magma-assisted rifting model by Buck (2004) describes how, if a steady supply of magma is available, extensional strain can be largely accommodated by dike intrusion at a fraction of the force required for tectonic faulting, which reduces activity on riftbounding border faults (also see Keir et al., 2006). Dike intrusions are enabled by the buoyancy of magma, which ascends from basaltic melt generation zones in the mantle toward a neutral buoyancy level at the base of the crust. Faulting and dike intrusion accommodate extension in the upper crust (Rubin and Pollard, 1988;Rowland et al., 2007;Baer et al., 2008;Trippanera et al., 2019), while released magmatic volatiles can weaken lithosphere through increased pore fluid pressures (Sibson, 2000;Reyners et al., 2007) and/or hydration mineral reactions (Moore and Rymer, 2007). All of these processes enhance the frequency of seismic activity and drive localization of fault-related strain (e.g., Muirhead et al., 2016;Chiodini et al., 2020).
Regions of active and evolving continental extension represent ideal locations to test the role of magma and magmatic volatiles in controlling rift processes Lee et al., 2016;Hunt et al., 2017;Oliva et al., 2019a;Kämpf et al., 2019). For example, the distribution of seismicity observed in tectonically active systems can provide unprecedented constraints on the 3D distribution of deformation associated with the magmatic, tectonic, and volcanic systems that interact to achieve continental extension (Geissler et al., 2005;Lambotte et al., 2014;Hurst et al., 2016;Nakai et al., 2017). These can be compared to geologically constrained models of state-of-stress from extensional stress and surface, internal, and subsurface density contrasts that load the plate (e.g., Karlstrom et al., 2010;Oliva et al., 2019b).
Geological, geodetic, geochemical and geophysical investigations of magma-rich rift systems in East Africa have revealed arguably the most dynamic magmatic-tectonic interactions observed in any such setting (e.g., Keir et al., 2006;Ebinger et al., 2010;Rooney et al., 2014;Wauthier et al., 2015). Combined geodetic and seismic observations reveal that magma systems situated within the center of rift segments can feed up to 45 km-long lateral dike injections (Wright et al., 2006). Stress changes associated with these dike intrusions drive shallow normal faulting (Rowland et al., 2007) and can trigger eruptions from magma chambers situated tens of kilometers from the initial point of magma injection (Wright et al., 2006;Ayele et al., 2009). Seismicity and InSAR observation of active intrusions in the Erta Ale rift segment (Ethiopia) reveal that ∼10 km-spaced magma chambers can synchronously feed eruptions through laterally connecting dikes . Observations in the Kenya Rift reveal transient uplift and subsidence from nearby sources at some volcanoes, which is interpreted to represent magma migration through laterally connected melt bodies (Biggs et al., 2009(Biggs et al., , 2016. This study utilizes a densely spaced seismic network in the southern Natron basin to provide unprecedented constraints on the 3D distribution of seismicity in an active magma-rich continental rift, which allows us to interpret the geometry and dynamics of the plumbing system. The early-stage Natron rift basin (∼3 Ma; Foster et al., 1997) contains Earth's only active natrocarbonatite volcano Oldoinyo Lengai and represents a region of recent dike and fault activity (Calais et al., 2008;Biggs et al., 2013). This volcano and a volcanic field are located in the interaction zone between a border fault-controlled rift segment and a developing rift magmatic segment (Muirhead et al., 2015Roecker et al., 2017;Weinstein et al., 2017), where mantle fluids drive hydraulic fracturing and associated fluid ascent along active fault systems Oliva et al., 2019a). Specifically, we report on ten months of seismicity and related focal mechanisms from the temporary SEISVOL network (Seismic and Infrasound Networks to study the volcano Oldoinyo Lengai) from 2019. These data are examined to discriminate zones of active intrusion and fluid migration and to understand how these processes control spatial variations in local stress states, rift kinematics and volcanic activity.

TECTONIC SETTING
Oldoinyo Lengai is located in the North Tanzanian Divergence, which is part of the East African Rift System (EARS, Figure 1). This is the longest continental rift worldwide with relatively slow extension rates, which are typically between 1-7 mm yr −1 in the present day in the North Tanzanian Divergence (Saria et al., 2014). In the North Tanzanian Divergence, magmatism initiated at ∼6 Ma (Mana et al., 2012), with evidence for faulting in the Natron basin region since 3 Ma (Foster et al., 1997), which in the present day is assisted by fluids rising through the plate as evidenced by earthquake swarms (Albaric et al., 2014;Lee et al., 2016;Weinstein et al., 2017). These interpretations are corroborated by soil CO 2 flux as well as carbon and helium isotope data from springs in the Natron basin, Tanzania, and Magadi basin, Kenya, revealing that magmatic CO 2 ascends to the surface along deeply penetrating faults in parts of the EARS Lee et al., 2017;Muirhead et al., 2020).
The North Tanzanian Divergence hosts the 200 km-wide Ngorongoro-Kilimanjaro Volcanic Belt, a diachronous belt of volcanoes that formed since ∼6 Ma, and shows an eastward younging trend (Le Gall et al., 2008;Mana et al., 2015;Muirhead et al., 2015, Figure 1). It consists of a large number of basaltic shields and composite cones that have a wide range of lava compositions including basanites, melilitites, basalts, nephelinites, phonolites, as well as carbonatites, implying multiple source depths and both asthenospheric and (dominantly) lithospheric mantle material (e.g., Dawson, 1992;Mana et al., 2015;Mollex et al., 2018). The unusual carbonatitic magmas are likely sourced from an enriched cratonic mantle lithosphere that experienced extensive metasomatism, with or without an asthenospheric melt component (Foley et al., 2012;Mollex et al., 2018;Muirhead et al., 2020). The magmas produced at these volcanoes exhibit high CO 2 contents, probably derived from small partial melting of a metasomatized mantle source (Fischer et al., 2009;Mana et al., 2015). Oldoinyo Lengai is situated at the southern end of the Natron basin and adjacent to its western border fault. In its immediate vicinity is Naibor Soito volcanic field, a monogenetic cone field with metasomatized mantle xenoliths, and the ∼1 Ma Gelai volcano (Dawson et al., 1995;Aulbach et al., 2011;Muirhead et al., 2016; Figure 1).
The 2007 events at Oldoinyo Lengai sparked new debate on the magmatic plumbing system and the volcano's connection with Gelai, the Naibor Soito volcanic field and surrounding rift faults (Biggs et al., 2013;Weinstein et al., 2017;Jones et al., 2019). Authors suggested that an undetected dike sourced from a deepseated magma chamber predated the initial fault slip detected during the 2007 unrest period (Baer et al., 2008;Biggs et al., 2013). From petrological observations, magma volume estimates and stress calculations, the 2007 eruptions could not have been sourced by a single shallow magmatic reservoir Biggs et al., 2013). Instead, the volcanic-tectonic event was most likely driven by a pressurized deep magma chamber, and magma was then distributed from there in the shallow crust, and fed into the Oldoinyo Lengai chamber, according to the local background stress (Baer et al., 2008;Biggs et al., 2013). Since then, a tomographic inversion by Roecker et al. (2017) imaged a low velocity, high Vp/Vs body interpreted as a single magma chamber between both volcanoes below a depth of 15 km. The combined tomography and seismicity data sets imaged two elliptical low velocity, high seismicity zones interpreted as sills that might have fed the 2007 dike intrusion, and triggered the eruption at Oldoinyo Lengai (Oliva et al., 2019a). Coulomb stress changes related to the volcanic-intrusive sequence may have also unclamped a section of the western Natron border fault, with the resulting slip detected on a local GPS network in early 2008 (Jones et al., 2019).

Data
The SEISVOL project (Seismic and Infrasound Networks to study the volcano Oldoinyo Lengai) is a temporary station deployment with 33 seismic and four infrasound stations. The deployed stations form a seismic network and array ( Figure 2) that were operative from February/March 2019 to June 2020. The seismic network and array comprised 10 Trillium Compact 120 and 20 Mark L4-3D seismometers, three 3D 4.5-Hz geophones and 33 Digos Cube 3-channel data loggers from the Geophysical Instrument Pool Potsdam (GIPP) and Goethe University Frankfurt. The infrasound stations consisted of four MB2017 sensors (ISTerre) and four Digos Cube 3-channel data loggers from Goethe University Frankfurt. All seismic stations recorded data at 100 Hz. Due to harsh field conditions, batteries as well as GPS antennas had to be replaced on a regular basis, and resulted in occasional downtimes for many stations ( Figure 3A). Here, we use data from March 1st to December 31st, 2019, as the remaining data were inaccessible due to the Covid-19 pandemic.

Earthquake Locations
We used the python software QuakeMigrate (Bacon et al., 2019) to detect, trigger, and locate earthquakes automatically. This software is based on the Coalescence Microseismic Mapping technique (Drew et al., 2013) and reports earthquake locations as well as P-and S-wave picks on the vertical and horizontal components, respectively. We analyzed the data using a bandpass filter of 1-25 Hz and the 1D velocity model from Albaric et al. (2010). We then filtered these results by the earthquakes' global covariance error in x, y and z directions and kept only earthquakes with less than 5 km horizontal and 7 km vertical uncertainties.
To obtain accurate hypocenter locations, we then relocated these chosen events with the probabilistic, nonlinear earthquake code NONLINLOC (Lomax et al., 2000) using the 3D-velocity model by Roecker et al. (2017) for the North Tanzanian Divergence. The final seismicity catalog consists only of earthquakes between −2.3°and −3.0°Latitude and from 35.7°t o 36.4°Longitude, and with location errors less than 5 km in the horizontal and 7 km in the vertical direction as calculated from final NONLINLOC location uncertainties and at least six picks. We report 6827 earthquakes with average location errors of 1.57 km in longitude, 1.59 km in latitude and 3.8 km in depth. We have ensured the location accuracy by comparing the automatically obtained earthquake locations with a subset of manually picked and located events, as well as comparing the differences between using a 1D and 3D velocity model. We can support the assertion by Weinstein et al. (2017) that differences in earthquake locations between both models are inconsequential and mostly within the horizontal and vertical errors.

Magnitudes
Magnitudes were calculated automatically in QuakeMigrate after Greenfield et al. (2018): after de-meaning and tapering the data, the instrument response is removed on each waveform. Then, waveforms are convolved with the response of a Wood-Anderson seismograph with the corrected gain of 2080 (Urhammer, 1982). For all network stations, the maximum zero-to-peak amplitude is picked on both horizontal components. Amplitudes are only picked if the maximum amplitude is greater than twice the standard deviation of a 10 s window of the data before the P-wave onset. Accordingly, only good recordings contribute to the calculated magnitude. These were then calculated on the local magnitude scale using the formula: where M i is the local magnitude of earthquake i, A ijk the amplitude recorded at each station j and component k, and the second terms denotes that we normed the magnitudes by a magnitude 3 earthquake in 17 km distance, using the hypocentral distance r ij and the attenuation parameters K and n (Hutton and Boore, 1987). K and n were derived for Tanzania by Langston et al. (1998) and are suitable given the small distance range for which they lead to equivalent magnitudes as compared to the parameters reported by Weinstein et al. (2017). Reported magnitudes are an average of all horizontal components per event (Greenfield et al., 2018).

Focal Mechanisms
We manually inspected the automatic earthquake catalog for suitable earthquakes for the calculation of focal mechanisms. We visually inspected waveforms for those earthquakes with more than 11 automatic P-and S-wave picks. We re-picked P-and Swave arrival times on unfiltered waveforms for those events, which show little noise contamination using the SEISAN software (Havskov and Ottemöller, 1999), and relocated these earthquakes with the in SEISAN implemented HYPOCENTER algorithm by Lienert and Havskov (1995). We used the 1D velocity model by Albaric et al. (2010) as SEISAN only permits 1D velocity models. Given our previous analysis, location differences to the 3D model are negligible. We continue to use only those earthquakes with azimuthal gaps smaller than 180°and at least eight P-wave polarity picks, which we picked on raw vertical traces. Fault plane solutions were then derived using the FOCMEC algorithm by Snoke et al. (1984), which is also implemented in SEISAN. This is a grid search algorithm to find the orientation of nodal planes using a double-couple solution. We started our analysis with a grid of 5°and tested whether a unique solution can be obtained. If this is the case, we refine the grid down to 1°w here possible. If a unique solution cannot be obtained, but solutions for a larger grid show a preferred mechanism, we manually picked SV/P amplitudes on the vertical component and included those in the inversion. We repeated our grid search for a grid of 5°and allowed for larger deviations in the amplitude ratios given large local attenuation effects. We only included those solutions that are unique for a grid of 5°, but refined the grid where possible. This procedure allowed us to constrain these fault-plane solutions accurately and enables us to report focal mechanisms for more earthquakes than using P-wave polarities alone (for categorization examples, see Supplementary Figure S1). In total, we report 281 fault plane solutions corresponding to three categories. 48% consist of solutions using only P-wave polarities, 48% are solutions using SV/P amplitude ratios additionally, and in 4% of cases we had to allow for one polarity error. We estimated the uncertainty of solutions by calculating a 95% confidence level for P/T plunges and trends using all solutions per event and found no systematic differences between small and large magnitude events.

Earthquake Locations
We located 6,827 earthquakes within and in the vicinity of our network ( Figure 3; Supplementary Table S1). On average, we report 22 earthquakes per day, but the number of running stations heavily biases the number of reported events ( Figure 3A). The locations of failed stations also impact the number of detected earthquakes. This effect is especially pronounced for stations at the network border, i.e. east of Gelai. During a period when the original setup of 20 stations was running, we detected 29 earthquakes per day on average. Figure 3A clearly shows two periods in time (April 13-15th, and October 10th) that greatly exceeded the average earthquake number per day. During these periods, two seismic swarms were observed, which will be discussed in (Swarms).
The observed seismicity is heavily clustered beneath volcanic centers ( Figure 3B). More than half of the earthquakes are located beneath the eastern flank of Gelai, which is also the location for the first earthquake swarm. There is a separate cluster beneath the Naibor Soito monogenetic cone field. Earthquakes also cluster north and south of Oldoinyo Lengai, but locations below the border fault are poorly constrained given the station setup. Scattered seismicity appears beneath Lake Natron, with two distinct clusters between the lake and adjacent to the western flank of Gelai. One of these clusters corresponds to a discrete seismic swarm recorded in October 2019. Most of the seismicity occurs at a depth of 5-10 km ( Figure 4A).

Magnitudes and B-Value
Local magnitudes M L range between −0.85 and 3.6. The largest magnitude earthquake is part of the first seismic swarm in April 2019 beneath the eastern flank of Gelai (see (Swarms), Figure 4). We sort the magnitudes into 0.1 M L bins and find that the four most populated bins lie between −0.2 and 0.2 M L , which corresponds to 40% of the data. We derive the magnitude of completeness (M C ) using the maximum curvature method and apply the 0.2 correction factor of Woessner and Wiemer (2005) to obtain an M C value of 0. We calculate the b-value by a robust, linear least squares fit to be ∼0.99 ( Figure 4B). We note that our b-value estimate changes slightly throughout the observational period, as the detection threshold for small magnitude events strongly depend on the number of running stations (see Supplementary Figure S2).
The calculated local magnitudes were used to estimate the cumulative seismic moment release throughout the 10 months observation period ( Figure 4C). Given the sparse number of calibrated earthquakes, we adopt the crude assumption that the moment magnitude M W M L . This assumption is supported by a comparison of local and teleseismically calculated magnitudes in the study region by Weinstein et al. (2017), who found that 2 M L 4.7 and 4.7 earthquakes on June 3, 2013 corresponded to body wave magnitude mb 4.5 and mb 4.6 events documented in the National Earthquake Information Center (NEIC) catalog. Indeed, mb is expected to closely resemble M w for earthquakes of mb ≤ 5 (Gasperini et al., 2013). Our results show that seismic moment release is nearly constant over time, but dominated by the 3.6 M L event on April 13, 2019.

Spatial Distribution
Here, we focused only on well-constrained events close to and within the network by removing all earthquakes at distances larger than 2.5 km from the seismic network boundaries. A total of 5,658 earthquakes occurred within the defined boundary, with horizontal and vertical errors of 1.5 and 3.7 km, respectively. The depth distribution of seismicity varies spatially within our focus site. To illustrate these depth variations, we project the seismicity along three profiles through all clusters ( Figure 5).
Profile X-X' transects the two clusters south and north of Oldoinyo Lengai ( Figure 5B). Seismicity in the south, close to the border fault, typically occurs between 8 and 13 km depth and gradually deepens northward, with earthquake depths ranging ∼13-17 km at the northern end of the profile. Similar earthquake depths are observed farther north of Oldoinyo Lengai below Lake Natron, with most events here ranging 12-20 km depth ( Figure 5B). Directly beneath Oldoinyo Lengai is a notable seismic gap. Magnitudes along this profile are mostly between −0.3 and −0.2 with one larger event at M L 3.27.
Profile Y-Y' runs southeast from the cluster north of Oldoinyo Lengai through the Naibor Soito volcanic field ( Figure 5C). In the center of the profile, the seismicity is located only in shallow depths of 5-10 km. Below the Naibor Soito volcanic field there are two clusters of seismicity, one at 5-10 km depth and another at 13-18 km depth. There is a noticeable gap between the deep cluster in the North (as seen in X-X') and the deep cluster beneath Naibor Soito. Toward the southeast end of profile Y-Y', seismicity is less abundant with scattered shallow (5-7 km) and deep (18-20 km) events. Overall, earthquake magnitudes along this profile are slightly smaller (∼0.1 M L ) than in profile X-X', with the biggest event M L 2.83. Profile Z-Z' runs from the southwest end of the Naibor Soito volcanic field northward along the eastern flank of Gelai ( Figure 5D). This profile reveals an area of intense seismicity. Earthquakes at the southwestern edge of the Naibor Soito volcanic field are generally shallow, with depths ranging from 3 to 10 km. As with profile B-B', shallow and deep seismicity clusters are present at 5-10 km and 13-18 km depth, respectively, below the Naibor Soito volcanic field. This is followed by an abrupt seismic gap at 13-16 km along the profile. Beneath the gap, seismicity is first between depths of 10-13 km, and then gradually becomes shallower toward the center of Gelai. Seismicity beneath Gelai is mostly between 3-10 km, and the largest event with M L 3.6 is visible at ∼ 28 km along the profile.

Focal Mechanisms
Following our procedure described in Data and Methods, we report 281 well-constrained focal mechanisms (see Supplementary Table   S2). We use the plunge of the P-and T-axis to categorize the solutions into the three main faulting types: strike-slip (blue symbols; P-and T-axis plunge 0-45°), normal (red symbols; Pand T-axis plunge 45-90°and 0-45, respectively) and reverse mechanisms (green symbols; P-and T-axis plunge between 0-45°and 45-90°; see Figure 6A). Fault plane solution are dominantly strike-slip (57%), with normal mechanisms making up 39% of the data, and rare occurrences of reverse mechanisms (4%). The mean direction of T-axes is NW-SE ( Figure 6B).
For ease of visualization, we show the best 50 solutions on a map only. We define those as being well constrained with a fine grid by using only P-wave polarizations, thus giving these focal mechanisms more importance than those constrained with amplitude ratios, and have at least 15 polarity picks (Figure 7, see Supplementary Figure S1 for categories). Focal mechanisms around Oldoinyo Lengai are purely strike-slip, with T-axis trends varying between N-S to E-W. Above the Naibor Soito volcanic field and at Gelai volcano, normal faulting dominates with minor strike-slip mechanisms. The northernmost focal mechanism at Gelai is part of the first swarm. T-axis directions are ∼ N-S or NW-SE. Three normal faulting events were constrained below the southern shore of Lake Natron. Owing to the network design, we cannot constrain focal mechanisms in other areas.
To study the direction of extension more closely, we plot the T-axis direction as a bar centered on the earthquake location and group the data into different depth bins ( Figure 8). Here, we omit all T-axes with confidence level widths above 40°while the remaining 82% have an uncertainty smaller than 10°. We constrain only a few mechanisms above 5 km depth of which most are located beneath Gelai in a NW-SE direction. This is more evident between 5-7 km depth and persists down to 9 km. T-axis directions below Naibor Soito volcanic field are very mixed throughout the 5-7 km, 9-11 km and below 13 km depth slices. North and East of Oldoinyo Lengai, T-axes can be divided into two groups: N-S and E-W directions. In 7-11 km depth, axes east of Oldoinyo are mostly N-S, while in 9-13 km, axes north of Oldoinyo are E-W.

Swarms
We observed two 1-3 days-long seismic swarms within the 10 months of data acquisition. During these swarms, the study region experienced a 5-to 8-fold increase in the number of earthquakes in a single day ( Figure 3A), most of which were confined to a <10 km 2 area ( Figure 9). The first swarm occurred April 13th-15th 2019 on the eastern flank of Gelai. We repicked all swarm events manually, including amplitudes, in SEISAN and relocated them in NLL (see Supplementary  Table S3). The swarm consists of 266 well-located earthquakes with a spatial distribution resembling a ∼2 × 2 km-wide and ∼11 km-deep pipe-like structure between 2 and 9 km depth ( Figure 9A). Mean location errors are 0.8, 1.1 and 4 km in latitude, longitude and depth, respectively. Among the very first earthquakes was a 3.6 M L event, which was located at the lower depth bound of the swarm pattern. The swarm has a b-value of 1.01 (see Supplementary Figure S3 for the magnitude-frequency relation and distribution of magnitudes over time). There is no clear pattern in the distribution of seismicity in space or time, with events distributed throughout the interpreted pipe-like feature occurring the entire duration of the swarm. We are unable to constrain a fault plane solution for the largest event (M L 3.6), due to a small arrival before the P-wave. The second largest event with M L 2.2, which occurred on April 15th at 05:47, was a normal faulting event. We constrain seven mechanisms in total, with four normal, two strike-slip and one reverse mechanism ( Figure 9A left). All mechanisms, except the reverse faulting event, have a T-axis trend of NW-SE.
The second seismic swarm comprises 151 events that occurred on October 10th, between Lake Natron and the western Gelai flank ( Figure 9B). We did not manually repick this swarm as it is north of the network and the two closest stations were not running. The locations and depth are likely less well-constrained, although average location uncertainties are 1 km in longitude, 1.5 km in latitude and 4.5 km in depth. Earthquakes in this swarm were located between 6-15 km depth and their distribution resembles a 2 × 2 km-wide pipe-like structure. Given the distance to the closest seismic stations, we likely do not observe many earthquakes below M L 0. Our b-value estimate is 0.72 (see Supplementary Figure S4 for the magnitude-frequency relation and distribution of magnitudes over time), and we cannot constrain any focal mechanisms for this swarm.

DISCUSSION
Here we discuss our observations in relation to previous studies. Seismicity observations from the southern Natron basin are sparse, and are only available for part of the 2007-2008 eruption (Calais et al., 2008;Albaric et al., 2010) and a 13months deployment in 2013-2014 (Weinstein et al., 2017;Oliva et al., 2019a), which marked a time of relatively little activity at Oldoinyo Lengai as inferred from satellite data (Coppola et al., 2016;Mirova, 2020). As discussed below, our results show that some seismicity patterns observed during these previous deployments are consistent over decadal timescales, while previously unrecognized patterns reveal shorter-term variability in local magmatic-tectonic processes.
Generally, seismicity clusters below the volcanic centers, and as such we define five areas of seismic activity: 1) Oldoinyo Lengai, 2) Naibor Soito volcanic field, 3) Gelai, 4) Crater Highlands, and 5) Lake Natron ( Figure 10A). Given our network design, the detailed 3D seismicity patterns are best constrained and show the greatest detail below areas (1), (2) and (3), which will be the focus of the following sections. The networks lacks resolution for earthquakes below the Crater Highlands (4), which is located at the SE edge of the network, and thus we largely omit these data from further discussion. Comparatively little seismic activity is observed below Lake Natron (5) associated with the east-dipping western border fault. This is consistent with the findings of Weinstein et al. (2017), and suggests that recent extensional strain accommodation has been facilitated through magmatic activity (e.g., dike, sill intrusions) rather than slip along border faults. The only major seismic swarm observed below Lake Natron correlates with a west-dipping fault at the western edge of Gelai volcano. This fault is characterized by sets of fault-parallel springs exhibiting high 3 He/ 4 He ratios and diffuse soil CO 2 discharges, indicative of tectonic degassing of deep mantle volatiles Lee et al., 2017).

Tectonic vs Fluid-Driven Seismicity in the Natron Basin
The southern Natron basin appears to have experienced persistent earthquake activity for the last 13 years, with likely more than 1,000 events of >0 M L occurring every year (Calais et al., 2008;Albaric et al., 2010;Weinstein et al., 2017; this study).
Our results indicate an increase in b-value from 2013-2014 of 0.87 (Weinstein et al., 2017) to 0.99, which is consistent with the 10-15 years eruptive cycle of Oldoinyo Lengai and increasing pressurization on the system. This persistent seismicity likely reflects the influence of a crustal magma plumbing system on the rifting process. For example, in many continental rift settings, Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 8 | Article 609805 crustal strain is accommodated by a combination of magmatism and faulting (e.g. Buck, 2004). Melting via decompression may be generated by large lateral gradients in topography at the lithosphere-asthenosphere boundary (LAB) with or without the influence of an anomalously hot upper mantle (e.g. Ebinger and Sleep, 1998). Basaltic magma is less dense than the mantle lithosphere, and rises to the base of the crust of equivalent density (Watanabe et al., 1999). Here, these magmas will temporarily pond as melt bodies that, depending on magma flux and the efficiency of heat transfer through the lithosphere, may induce further heating and melting of surrounding rock before eventually freezing as an underplate (e.g., Thybo and Artemieva, 2013;Karakas and Dufek, 2015;Rooney et al., 2017). Heating and volatile release in these regions serve to localize magmatism, and this heating may reduce crustal strength and further localizes strain (Bialas et al., 2010). Crustal thinning and intrusion of magma cause changes in crustal density, which, in addition to rift basin topography and volcanic construction, generates loads that moderate the tectonic state-of-stress (Karlstrom et al., 2009;Maccaferri et al., 2014). Flexural isostatic compensation for the loads leads to plate bending, and depth-dependent state of stress. Intrusions Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 8 | Article 609805 11 moderate the stress field in the crust and in large part determine the orientation and location of subsequent intrusions (e.g., Gudmundsson, 1995;Buck et al., 2006;Karlstrom et al., 2009).
Crustal seismicity in magmatic rift zones can therefore reflect a variety of processes, from tectonic-to fluid-driven faulting, magma transport and storage, and pre-, syn-and posteruption dynamics. Broadly, observed seismicity can roughly be differentiated into 1) fluid processes (e.g., magma ascent in conduits, gas pressurization in conduits, bubble formation) that lead to low frequency, long-period events; or 2) brittle shear failure or shear failure with dilatation, which produce highfrequency events (e.g. McNutt, 2005). Given the comparatively high frequency content of the earthquakes reported here, we assume most earthquakes are associated with brittle deformation, although fluid (over)pressures likely play an important role driving much of the observed seismicity (i.e., Lindenfeld et al., 2012b). Here, we examine seismicity patterns and focal mechanisms for 1) Oldoinyo Lengai, 2) Naibor Soito, and 3) Gelai shield volcano to characterize the extent to which both faulting and magmatic processes are active between eruptions.  Figures 5, 10) and there is little seismicity in the vicinity of the second, late 2007 dike intrusion inferred from surface deformation patterns (Biggs et al., 2013). In Iceland, deep crustal seismicity (e.g., below the brittle-ductile transition down to 25 km) has been interpreted to represent brittle failure resulting from melt migration (e.g. Hudson et al., 2017;Greenfield et al., 2018). Weinstein et al. (2017) assume a weakly hydrated mafic granulite lower crustal composition based on crustal xenoliths from the study region, and estimate the brittle-ductile transition at ∼16 km depth. Seismicity could also be caused by fluid-driven fracturing due to magmatic degassing (e.g. Lindenfeld et al., 2012b) or be related to the down-dip extension of the border fault below the volcano, along which deeply-sourced fluids such as CO 2 abundantly escape Muirhead et al., 2020) and may cause fault zone weakening (Reyners et al., 2007). To conclude, seismicity at 15-20 km could thus be indicative of volatile or magmatic fluid movement. Focal mechanisms at Oldoinyo Lengai are purely strike-slip. Similar strike-slip kinematics have been observed at other volcano-tectonics settings (e.g., Ágústsdóttir et al., 2016;Haug et al., 2018) and may be explained through a variety of mechanisms. In Afar, Pagli et al. (2019) show how linking rift segments can produce strike-slip tectonics. Dering et al. (2019) confirm the recent geophysical evidence for dike-parallel shear failure that can occur in the host rocks surrounding intruding dykes. Considering magmatic interactions, Hill (1977) suggests that observed strike-slip mechanisms result from the joining of offset opening cracks where hypocenters are located at the edges of interacting dike segments. Similarly, moment-tensor solutions from the 2014 Bárðarbunga-Holuhraun dike intrusion support strike-slip motions parallel to the advancing dike tip (Ágústsdóttir et al., 2016). Roman (2005) support a model where strike-slip earthquakes occur in the walls of inflating dikes. In this example, dike inflation causes a 90°r otation of P-axes from that predicted by the regional tectonic stress field, which allows for strike-slip motions. In all these examples, the regional extensional stress deviates locally during periods of magma influx, leading to strike-slip mechanisms as the predominant slip mode. Support for a local magma-driven stress field around Oldoinyo Lengai comes from the spatial pattern of P-axes for strike-slip mechanisms around the volcano. These exhibit a subradial pattern, with the trends of the P-axes converging near the center of the volcanic edifice ( Figure 11) in accordance with vent lineament data associated with the Naibor Soito volcanic field (Muirhead et al., 2015). This local stress pattern would also be enhanced by loading associated with both Oldoinyo Lengai and the Crater Highlands (Oliva et al., 2019b). Notably, the strike of the late 2007 dike intrusion beneath Oldoinyo Lengai also matches a predicted radial pattern, which may be best explained by a combination of the effects of magma overpressures in the volcanic feeder system and edifice loading (Biggs et al., 2013;Roman and Jaupart, 2014;Muirhead et al., 2015). Numerical studies investigating these effects (e.g. Karlstrom et al., 2009;Oliva et al., 2019b) reveal that edifice loading exerts a greater influence on the local stress state in the shallow upper crust, whereas magmatic overpressures exert the dominant stress control in vicinity of magma chambers (1) yellow for Oldoinyo Lengai, (2) blue for Naibor Soito, (3), cyan for Gelai, (4) magenta for Crater Highlands, (5) green for Lake Natron. The two swarms are shown in orange and gold. Inferred dikes are shown as red lines (Calais et al., 2008;Biggs et al., 2013) and roughly dimensioned magma chamber in shaded-red inferred from Roecker et al. (2017). Purple arrow shows view angle in (B). (B) Vertical Cross-section of (A) looking −20°N. Circles show earthquakes and are colored by group and scaled by their magnitude.
Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 8 | Article 609805 at depth. These processes play a critical role in driving local stress rotations and controlling fracture and dike orientations (e.g., Kervyn et al., 2009). In our study region, the depth of the strike-slip events (<9 km) and their proximity to magma pathways below the volcano support a role of magmatic pressures in controlling local stress deviations at depth, thereby capturing the effects of the Oldoinyo Lengai magmatic system on the surrounding stress state. Of note, our interpretation is somewhat limited by the fact we cannot utilize surface deformation patterns (e.g. InSAR data) as an additional independent dataset to support our arguments above, as the source depths of the analyzed focal mechanisms are not conducive to any detectable surface deformation.

Naibor Soito
Beneath Naibor Soito volcanic field, seismicity is divided into a shallower (5-10 km depth) and deeper cluster (12-18 km depth), with comparatively little seismicity occurring at depths of 10-12 km, except for a fine structure seemingly connecting both (Figures 5, 10B). Both clusters were observed in 2012-2013 (Weinstein et al., 2017). The seismicity cluster between 12 and 18 km occurs on top of a S-wave low velocity and high Vp/Vs zone, which has been interpreted previously as a lower crustal magma body . In agreement with previous authors, we assume that the seismicity detected adjacent to and above this low-velocity zone represents fluidtriggered fracturing, resulting from one or a combination of devolatization of magmas and the migration of melt from lower crustal reservoirs Weinstein et al., 2017). The shallow, disc-shaped seismicity pattern was interpreted as a complex of pressurized, hot sills likely fed by the deeper magma source (Weinstein et al., 2017). These shallow sills in-turn may have sourced the 2007 dike intrusion (Calais et al., 2008;Weinstein et al., 2017). The recurring non-double couple source mechanisms suggest that numerous intrusions from the same complex have occurred since 1995 (Oliva et al., 2019a). Thus, our and previous observations support a central, deep magmatic source which drives the volcanism locally. Below the Naibor Soito volcanic field, our high-resolution observations allow us to distinguish a complex pattern of focal mechanisms and T/P-axis directions. Above the magma chamber (below 13 km) and in the pressurized sill complex (5-7 km; 9-11 km), there is no dominant mechanism or stress direction. Between 7 and 9 km depth, T-axes are ∼ N-S (or nearly E-W P-axes, Figure 11), which represents a rotation of local stress away from the regional ESE-WNW extension direction (Saria et al., 2014). Both observations are best explained by magmatic processes associated with melt transported from the deeper magmatic body to the shallow sill complex. While the pressurized magmatic bodies exert a changing stress field depending on their state (i.e. inflation or deflation), earthquakes in the walls of the inflating intrusions show a local horizontal stress change according to Roman (2005). This is also supported by the findings of Oliva et al. (2019b), who predict stress rotation between Oldoinyo Lengai, Naibor Soito and Gelai driven by magmatic intrusions and magma bodies at depth. These fluid driven stress-changes as well as strong shallow and deep seismicity patterns are among the most persistent, clearly imaged features of our dataset through our observational period.

Gelai
All temporary seismic networks have recorded persistent seismic activity beneath Gelai volcano since 2007 (Calais et al., 2008;Albaric et al., 2010;Weinstein et al., 2017). However, the closer station spacing in our network allows examination of seismic activity in greater detail than previous studies (Figures 5, 10). Locations reported by Calais et al. (2008) spread widely, which may largely reflect location uncertainties associated with locating events outside of the network. Albaric et al. (2010) found that seismic activity was primarily focused beneath the southern flank of Gelai, while Weinstein et al. (2017) reported activity beneath the southern and western flanks of Gelai. Overall, these previous results combined with our newly located earthquake data suggest that activity may have migrated northwards below Gelai over time, as the recent activity now reaches the northern end of Gelai ( Figure 5). Beneath Gelai, earthquakes show predominantly normal faulting mechanisms (Figure 7) with T-axes subparallel to the regional extension direction (Saria et al., 2014;Muirhead et al., 2015;Weinstein et al., 2017), which is consistent with faulting in extensional rifts elsewhere (e.g. Keir et al., 2006;Lindenfeld et al., 2012a). These mechanisms contrast those at Oldoinyo Lengai and the Naibor Soito volcanic field. The strike of faults, dikes and previously determined focal mechanisms at Gelai typically indicate WNW-ESE extension (Calais et al., 2008;Albaric et al., 2010;Muirhead et al., 2015;Weinstein et al., 2017), which is consistent with geodetic data (Saria et al., 2014).
Gelai volcano is also associated with the April 2019 seismic swarm, which exhibits a very narrow, pipe-like shape. This swarm is in line with, and located ∼10 km north of, the 2007 dike intrusion as imaged by InSAR (Calais et al., 2008), which seems to have initiated the recent recorded Gelai seismic activity (Weinstein et al., 2017). However, seismicity associated with dike emplacement and propagation usually follows along a plane and migrates spatially through time in the direction of magma movement (e.g. Ebinger et al., 2010;Ágústsdóttir et al., 2016), and thus the 2019 swarm activity is not interpreted to relate directly to dike emplacement. Though not exhibiting a clear main shock after shock or clear temporal sequence, the 3.6 M L event is the first and second deepest event which seems to start the sequence. Afterward, other earthquakes rupture along the pipe, which is then active along the entire segment. Focal mechanisms show normal faulting representative of the regional stress field. Similar pipe-like seismic swarms were observed by Lindenfeld et al. (2012b) in the Rwenzori region of the Albertine Rift (Uganda) of the western branch of EARS, and were interpreted to originate from fluids and gases, sourced from a magmatic body, rising through the plate. Accordingly, this swarm may be caused by degassing along a currently hidden fault and is feasible given the high amount of CO 2 degassing in the entire Natron rift , Lee et al., 2017Muirhead et al., 2020). The shape and duration of the April 2019 swarm resembles the October swarm located between Gelai and Lake Natron which, given the nature of the seismicity and associated mantle volatile release at the surface (i.e., Lee et al., 2017), we attribute to ascent of deep mantle volatiles along a pre-existing fault. As crustal fluid speeds based on migrating seismic swarms are typically 10-100 m/day (Ibs-von Seht et al., 2001;Hariri et al., 2010), we assume that these earthquakes represent a mechanical response to subtle pressure changes in the fluid system due to changes in fault permeability (e.g. due mineral precipitation), in the local and/or regional stress state, and/or in the deep fluid flux (Sibson, 2000). The 3D seismicity pattern beneath Gelai cannot be explained by slip on a single fault plane. Given our observations, the strong microseismicity in the last ∼13 years and repeating earthquake swarms, indicate that fluids and volatile release must play an important role in accommodating the extensional strain, as is evident in this and other rift settings (Lindenfeld et al., 2012b;Wright et al., 2012;Muirhead et al., 2016;Weinstein et al., 2017). Diffuse, high magmatic CO 2 degassing from lower crustal intrusions throughout the entire Natron rift could cause significant earthquake activity Weinstein et al., 2017). In fact, the latter suggested that the recent locus of strain accommodation has been focused in the central rift valley below Gelai, and that stress field rotation and magmatism below Oldoinyo Lengai and the Naibor Soito volcanic field have facilitated stress transfer from the western Manyara-Natron border fault to a newly developing rift segment beneath Gelai.

Plumbing System
Our detailed 3D seismicity patterns and focal mechanisms allow us to map the magmatic plumbing system associated with Oldoinyo Lengai volcano in unprecedented detail. The persistent pipe-like aseismic structures down to depths of 12 km, as well as well resolved areas of high seismicity, are particularly striking ( Figure 10B) and may have implications about the long-debated plumbing system. We have summarized our findings in a sketch (Figure 12), which highlights our conceptual interpretation.
There is a distinct aseismic zone beneath the Oldoinyo Lengai edifice, which extends to a depth of 10 km and becomes successively narrower with depth ( Figures 10B, 12). This seismic gap may represent the region of persistent partial melt (i.e., Hudson et al., 2017) that connects the volcano to its deeper plumbing system, and/or has heated nearby country rock resulting in ductile rock behavior for extended time periods in otherwise brittle crust (Castaldo et al., 2019;Parisio et al., 2019). Indeed, renewed and persistent volcanic activity at Oldoinyo Lengai since late 2018, as inferred from satellite data (Massimetti et al., 2020), require open conduits transporting volatiles and melt to the surface. From tomography and complementary receiver function studies (Plasman et al., 2017), the volcano is underlain by a low S-wave anomaly and high Vp/Vs ratios at 5 km depth, which continue well under the rift shoulder  and may represent zones of partial melt below the volcano. Alternatively, the spatially varying high and low Vp/Vs and aseismic pipes in this region may, at least in part, reflect the sustained hydrothermal system interpreted for this volcano (Mollex et al., 2018), which is closely associated with its underlying magma plumbing system.
The exact spatial extent and overall geometry of partial melt bodies and hydrothermal fluids below Oldoinyo Lengai volcano are challenging to resolve. Surface deformation observed in InSAR data in late 2007 were modeled as a deflating shallow magma chamber at 3 km depth beneath Oldoinyo Lengai (Biggs et al., 2013), and chemistry from prior eruptions has been previously considered to reflect a two-level system, with magma storage at 0.5 and 3.3 km depth (Petibon et al., 1998). In all, some models of nephelinite and natrocarbonatite magma generation at Oldoinyo Lengai require only one open magma system (de Moor et al., 2013), while others necessitate at least three distinct zones of melt accumulation and storage (Berkesi et al., 2020). Regardless of these uncertainties, our high-resolution observations of earthquake activity below the volcano are best explained by a high temperature zone of partial melt and hydrothermal circulation underlying the edifice, where this zone is too ductile to produce any earthquakes (e.g. Thurber, 1984;Bryan et al., 1999;Castaldo et al., 2019, Figure 12). Additionally, our seismicity data and other geophysical constraints (i.e., Roecker et al., 2017) suggest that the inferred mid-crustal magma chamber for the volcano (i.e., Mollex et al., 2018) might either be located in the aseismic zone beneath Oldoinyo Lengai, or just above the deep-seated seismicity beneath Naibor Soto, where we image a circular aseismic feature (Figures 10, 12).
A second, pipe-like seismic gap separates seismicity patterns between Naibor Soito and Gelai (Figures 10, 12). It is particularly interesting, as there are no detectable surface displacements accompanying volumetric expansion in the upper crust since 2007, and the ∼2.4 m-wide 2007 dike (Biggs et al., 2009) would have long since cooled and solidified. This aseismic zone may instead mark a region of 1) hydrothermally altered rock following the 2007 dike intrusion, 2) a region of currently recognized weak rock or 3) part of the hydrothermal system (e.g., Mollex et al., 2018;Parisio et al., 2019). Considering the large amount of earthquakes beneath Gelai, which may be driven in large parts by fluids or degassing, this aseismic zone may also facilitate fluid transfer from the deeper parts of the plumbing system to the shallow crust beneath Gelai.
In conclusion, we cannot differentiate between different geochemical models of the natrocarbonatite genesis, but our results constrain the magmatic plumbing system of the Natron rift by interpreting areas of high and low seismicity as fluid or melt transport and locally heated, ductile zones, respectively. Focal mechanisms show the influence of magmatic intrusions and magma bodies, which together with seismicity may point to different localities for magma differentiation.

Interactions Between Rift Segments
Our observations are consistent with conceptual models of magmatic rifting revealing enhanced vertical permeability and associated fluid migration within accommodation zones between rift segments (e.g. Corti, 2004;Rowland et al., 2007;Muirhead et al., 2015; Figure 13). Uniquely, the North Tanzania Divergence offers the opportunity to explore strain accommodation between border fault-controlled extension in the Manyara basin, to largely magma intrusion-controlled extension within the central Natron basin, where a new, magmatically controlled segmentation is developing (e.g., Ebinger and Casey, 2001;Muirhead et al., 2015;Oliva et al., 2019a). Stress distributions in border fault-controlled segments are largely controlled by the shape of the basin and flank topography, and the Moho geometry, whereas segments with crustal magma chambers and volcanoes are controlled by crustal thinning, magma bodies, and volcanic loading (Karlstrom et al., 2009;Muirhead et al., 2015;Oliva et al., 2019b).
We image how the complex stress field in the inter-segment transfer zone drives the development of radial and rift-oblique dikes, faults and fracture systems, and how interconnected zones of magmatism feeding into these transfer zones facilitate strain transfer between the rift segments ( Figure 13). Our seismicity and fault plane solutions reveal that the high flux of gas-rich magma into these regions produces overpressured magma reservoirs that locally perturbed stress fields. This is supported by vent alignment (Muirhead et al., 2015), InSAR data (Biggs et al., 2013;Wauthier et al., 2013) and FIGURE 12 | Plumbing system of Oldoinyo Lengai with seismicity. The orange shaded areas are seismic gaps interpreted as part of the mushy plumbing system, overlain with small reddish ellipses denoting sill complexes. Black lines denote deep plumbing system, which connect the hydrothermal system to the deep magma chamber inferred from Roecker et al. (2017 temporally consistent earthquake source mechanisms (Oliva et al., 2019a). Furthermore, the shape of crustal thinning, the rift basin and flank topography guide zones of magma upwelling and storage, which promotes the transfer of strain to central volcanic systems (Oliva et al., 2019b). Extension-normal dikes inject into the rift to accommodate regional tectonic extension, while also focusing strain into newly developing rift segments (i.e., below Gelai volcano; Weinstein et al., 2017).
FIGURE 13 | Idealized conceptual illustration of magmatic-fluid-tectonic interactions associated with the volcanic plumbing system of the south Natron basin. The structure of the plumbing systems is based on data and interpretations presented in the current study, as well as previous seismicity and geophysical studies by Roecker et al. (2017) and Weinstein et al. (2017). Crustal thicknesses are broadly based on Plasman et al. (2017) and surface topography is from the 90 m SRTM DEM dataset. (A) Surface topography and general crustal structure of the south Natron basin, highlighting the position of the key volcanic features. Letters refer to the presented crustal sections in (B-D). (B) General distribution of seismicity at the northern end of the seismic network. Deep seismicity occurs below Lake Natron, with seismic swarm behavior associated with fluid overpressures and fluid migration through the crust. Helium isotope data of dissolved spring volatiles and carbon isotopes of diffusely degassing CO 2 support mantle fluid flow through these fault systems , Lee et al., 2017Muirhead et al., 2020). (C) General distribution of seismicity through central regions of the seismic network. Seismicity patterns generally follow those in B and support fluid-driven earthquakes, although no specific seismic swarms have been identified. Cooled dike intrusions represent regions of past magmatic rifting, although no significant melt or active diking is inferred below Gelai volcano currently. (D) General distribution of seismicity throughout the southern end of the seismic network, with deformation patterns strongly influenced by the presence of melt and fluids. The infrequent occurrence of seismicity in the lower crust below the Naibor Soto volcanic field support the presence of lower crustal magma bodies, consistent with the crustal velocity structure presented by Roecker et al. (2017). Distinct seismic gaps below these volcanic systems represent vertical zones of partial melt and/or thermally altered regions, and represent possible connections between partial melt bodies at different crustal depths (Weinstein et al., 2017). Melt bodies and topographic loads impose local stress perturbations (Oliva et al., 2019b) indicated by variations in fault plane solutions observed in this system. Persistent low magnitude seismicity throughout the region supports a model of earthquake activity driven by the percolation of melt and magmatic volatiles, consistent with full moment tensor analyses by Oliva et al. (2019a) who found various non-double couple earthquake source mechanisms.

CONCLUSION
We report seismicity and focal mechanisms from a ten-month deployment in 2019 in the southern Natron basin, which encapsulates recent volcanic activity at the only active natrocarbonatite volcano Oldoinyo Lengai worldwide and repeated seismic swarms consistent with volatile release and dike intrusion. Seismicity patterns beneath Oldoinyo Lengai and Naibor Soito reveal a complex magmatic plumbing system in which melt is transported from a deep-seated magma chamber to the surface, and focal mechanisms reveal the influence of the edifice load, magmatic intrusions at depth and magma-driven stress changes. Seismicity at Gelai is consistent with fluid release and degassing of CO 2 and focal mechanisms point to strain accommodation within the regional extensional stress field. Our results indicate that the southern Natron basin is a segmented rift system in which strain is transferred from the western border fault to a developing rift segment. The rift basin functions as a magmatic accommodation zone where fluids preferentially percolate vertically and thus prescribe the location of magma intrusion, storage and eruptive centers. These findings provide a framework for understanding and modeling complex interactions between tectonic and magmatic processes during continental rift initiation and development.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found at https://dx. doi.org/10.6084/m9. The dataset is archived at GEOFON Data Center, https://dx.doi.org/10.14470/4W7564850022.

AUTHOR CONTRIBUTIONS
MCR is the PI of the SEISVOL project, processed the data, composed Figures 2-12, and wrote the first manuscript. JDM composed Figures 1, 13. MCR, ASL, FL, and GR acquired the field data, with substantial support by EOK, MCR, and GR conceived and designed the seismic and infrasound networks. All authors contributed to the interpretation of the results and writing of the paper.

FUNDING
MCR is funded by the DFG (German Research Council), grant number RE 4321/1-1. JDM received funding through an NSF EAR GeoPRISMS Program grant 1654518. Additional support for fieldwork and instrumentation was provided through funds from Goethe University to GR.