Seismicity of the Bora‐Tullu Moye Volcanic Field, 2016–2017

The Bora‐Tullu Moye (TM) volcanic field is a geothermal energy prospect in the central Main Ethiopian Rift, but little is known about the seismicity of the region. Here we document seismic activity between February 2016 and October 2017, locating more than a 1,000 seismic events of local magnitude 0 to 2.7. This provides new insights into fluid movement and deformation beneath what we only now realize is a complicated volcanic system. A discrete cluster of events lies beneath TM, but, surprisingly, most of the seismicity lies in two clusters that are beneath neither the Bora or TM edifices. In these regions, we use earthquake cluster orientations, fault plane solutions, and fast seismic shear‐wave orientations to show that seismicity is triggered by hydrothermal circulation of fluids along preexisting fractures. The fractures trend in multiple directions and are, in general, not parallel to rifting related structures. Instead, the fractures are parallel to structures created during previous caldera forming eruptions at both Bora and TM. Highly fractured regions such as this could be attractive targets for geothermal power generation. We estimate a minimum depth for a magmatic body beneath TM to be 6.5 km using the mapped brittle‐ductile transition. Frequency analysis of the earthquake waveforms reveal the bulk of the events to be volcano tectonic, but some low‐frequency seismicity is present at a depth of 5 km beneath the TM edifice triggered by high pore fluid pressures.


Introduction
The Main Ethiopian Rift (MER) stretches from the Afar triple junction in the north to the Turkana depression in the south. More than 31 volcanoes are located within the rift (Global Volcanism Program, 2013), and approximately 11% of the population of Ethiopia live within 30 km of a volcano; making it important to understand the hazard posed by these volcanoes (Aspinall et al., 2011). Additionally, volcanoes in the MER have the potential to generate large amounts of low-carbon geothermal power, a resource that is not well-exploited. There is only one operating geothermal power plant in Ethiopia, located on Aluto volcano (Teklemariam et al., 1996), although more are being developed. Hydrothermal systems are likely to be driven by subsurface magma replenishment (Hill et al., 1985), making it important to understand the magmatic plumbing systems of volcanoes in order to assess their geothermal power potential. However, while regional scale geophysical surveys have reported on the overall structure of the MER, few studies have explored the detailed structure and magmatism of individual volcanoes.
This study presents the first detailed picture of the magmatic, hydrothermal, and fracture systems at Bora-Tullu Moye (B-TM), as part of the Natural Environmental Research Council-funded RiftVOLC project. We used a seismic network deployed at B-TM between February 2016 and October 2017 in order to detect and locate seismicity from around the volcano, characterizing the seismic activity during this period. We use the catalog of earthquakes to generate a 1-D velocity model suitable for the local region, compute focal mechanisms, and image the anisotropic structure of the crust. In addition, we use the spectral content of the earthquakes to map spatial changes in the earthquake source processes. These data allow us to image the state of stress around the volcano as well as the active fracture sets facilitating fluid flow through the volcano. This study serves as a benchmark upon which further studies can be built.

Tectonic and Volcanic History
Spreading between the Nubian and Somalian plates is currently accommodated by the MER at~6 mm/yr in an approximately east-west direction (Bendick et  Supporting Information: • Supporting Information S1 • Data Set S1 • Data Set S10 Correspondence to: T. Greenfield, tmgreenfield@gmail.com Rifting initiated in the central part of the MER in the Miocene (Bonini et al., 2005;Woldegabriel et al., 1990;Wolfenden et al., 2004) and was bounded by large border faults that accommodated the bulk of the plate spreading. Since 2 Ma, there has been a shift in the style of rifting with a transition to magmatically facilitated rifting resulting in relatively short and small offset intrarift faults (Corti et al., 2013;Ebinger & Casey, 2001). These right-stepping en-echelon faults form the present-day Wonji Fault Belt (WFB; Pizzi et al., 2006). The WFB is the also locus of most of the Quaternary-Recent volcanism, including major volcanic centers near B-TM such as Boset, Gedemsa, and Aluto, and aligned monogenetic cones (Ebinger & Casey, 2001). South of Aluto however, the narrow zone of intrarift faulting of the WFB is less well defined or absent, though major volcanic centers are still present on the rift floor (e.g., Corbetti).
B-TM is composed of two volcanic centers: Bora and Tullu Moye (TM). TM is located within the WFB, while Bora is offset to the west and closer to the center of the rift. The eruptive histories of both centers are poorly constrained, with few dates available to constrain eruption rates. Ignimbrites exposed in the WFB to the south of TM give the ages of large, possibly caldera forming eruptions between 115 ka and 1.8 Ma (Bigazzi et al., 1993). Eruptive products are a mix of silicic pyroclastics, obsidian, and basalt lava flows. Geochemical analysis of the silicic products reveal significant differences between the pantellerites found around Bora and commendites around TM, suggesting that each center has its own magmatic plumbing system (Fontijn et al., 2018).
TM is the more recently active, with fresher looking trachytic flows emanating from rift parallel fissures. The freshest looking flows may be related to eruptive activity between 100 and 200 years ago (Di Paola, 1972;Gouin, 1979). Significant horst and graben structures are observed extending to the north and south of TM and indicate the location of the WFB in this area ( Figure 1). In contrast, Bora is much less heavily faulted.

Crustal Structure and Seismicity
Little previous work has been done on the detailed structure of the B-TM volcanic region. However, some regional geophysical studies have included this region of Ethiopia and have constrained seismically active locations and crustal structure (e.g., Daly et al., 2008;Keir, Ebinger, et al., 2006;Keranen et al., 2004;Maguire et al., 2006).   Keir, Ebinger, et al. (2006) are colored gray and scaled by magnitude.
The P-wave crustal structure of the MER was imaged by Maguire et al. (2006) using wide angle reflection and refraction seismology. B-TM lies within high-velocity (6.3 km/s at 10-km depth) regions defining the presentday rift axis. The velocity structure within the rift showed low velocity (3.3-4.8 km/s) sediments near the surface extending to 5-km depth and P-wave velocities slightly faster (0.5 km/s) than normal continental crust within the rift (6 km/s). The best estimate for crustal thickness in this area of the MER, constrained using Moho refracted and reflected waves in combination with receiver function analysis, was 35 km (Cornwell et al., 2010;Maguire et al., 2006).
Local seismicity, detected between 2001 and 2003 using a regional scale seismic network, was characterized by two clusters of earthquakes (Keir, Ebinger, et al., 2006). One cluster was located beneath the TM edifice at depths between 4 and 11 km; the other, more active cluster was located south of Bora and extended between the surface and 20-km depth. It is worth noting that the closest seismic station to B-TM was more than 15 km away meaning the depths reported by Keir, Ebinger, et al. (2006) may have large errors.
While no detailed seismic studies have been published on B-TM, a dense seismic network deployed between 2012 and 2014 instrumented Aluto, a volcano located 50 km to the south of the study area (Wilks et al., 2017). The paucity of local seismicity recorded in the MER makes Aluto the best volcano to compare our study of B-TM to. Within 15 km of Aluto, more than 1,300 events were clustered into a weakly linear orientation parallel to the rift and centered on Aluto. In depth, most (53%) earthquakes were shallower than sea level, and only 7% were located deeper than 9 km. The seismicity revealed a shallow geothermal reservoir, additionally constrained using magneto-telluric methods (Samrock et al., 2015), at depths shallower than sea level. The seismicity within the reservoir was seasonally modulated by loading of the nearby lakes . Wilks et al. (2017) proposed a heavily intruded layer between 2 and 9 km below sea level (bsl), with volatile-rich lenses. The top of this layer hosts the best-fitting deformation source imaged using satellitederived deformation measurements between 2004 and 2010 (Biggs et al., 2011;. Two cycles of uplift and subsidence were imaged and interpreted to be caused by a physical connection between the magmatic and hydrothermal systems. The intrusion of magmatic fluids at a depth of 5 km caused uplift which relaxes by volatile exsolution and degassing. Below 9 km bsl, the abrupt change to relatively little seismicity suggests that 9 km bsl is the brittle-ductile transition beneath Aluto. Below this, a mushy, aseismic magma reservoir is proposed (Wilks et al., 2017).
Eighty kilometers southwest of Aluto, the volcano Corbetti has also been the site of geophysical observations. Surface deformation measurements have revealed sustained uplift of the volcano since at least 2012 (when measurements began; . The depth of intrusion was calculated to be 6 km and coincided with the location of a conductive region (Gíslason et al., 2015). Seismicity measured using a fairly sparse network recorded only 224 earthquakes over a 2-year period during uplift. The seismicity was located around the high conductivity region within more resistive areas of the crust .
Little is known about the long-term surface deformation of the B-TM region. Between 2004 and 2011, two pulses of uplift and subsidence were observed using Interferometric Synthetic Aperture Radar (InSAR) (Biggs et al., 2011). Uplift of 5 cm was focused at the Bora volcanic center, deformation that was best fit by a shallow (<2.5 km) penny-shaped crack with a radius of 3-10 km. While the observed displacement was not as much as that observed around Aluto or Corbetti, the presence of deformation shows that Bora could currently be active, and magma could be currently intruding into a shallow magma reservoir.

Seismic Anisotropy
The dominant mechanisms for upper-mantle and crustal anisotropy in Ethiopia are summarized in Kendall et al. (2006) and Hammond et al. (2014). The preferential alignment of segregated melt at depth holds the signature of strain accommodation along the plate boundary (Holtzman & Kendall, 2010). There is likely to be strong anisotropy due to the layering of sediments and volcanics in the upper 5 km of the rift valley (Bastow et al., 2010), but vertically propagating body waves, as used in our study, will not be sensitive to this style of anisotropy.
More recently, Nowacki et al. (2018) studied the SWA around Aluto. The observed pattern of SWA was complex, with interactions between two fracture systems used to explain variations in ϕ near the edifice. The amount of SWA was found to increase toward the center of the volcano reaching a peak of 4%. The SWA was also found to be concentrated at depths less than 3 km bsl. Nowacki et al. (2018) concluded that high gas overpressure in a geothermal reservoir maintains fluid flow in multiple fracture sets. Lloyd, Biggs, Birhanu, et al. (2018) presented observation of seismic anisotropy near Corbetti volcano. They observed fast shear-wave polarizations that parallel preexisting, cross-rift structures. Wadge et al. (2016) discuss the broader implications for the state of stress in East African Rift volcanism, including constraints from such shear-wave splitting (SWS) observations. Based on a range of observations, including anisotropy, they conclude that the state of stress in the volcanic systems of the East African Rift is complicated and often affected by preexisting structures and temporal changes in magmatism.

Seismic Network
A network of eight 3-component, broadband seismometers was deployed around B-TM in February 2016. It was further augmented with five additional seismometers in February 2017 ( Figure 1). All instruments were either Guralp CMG-6TD (30 s corner frequency) or Guralp CMG-ESPCD (60 s corner frequency) and were buried at least 1 m deep in the ground. Accurate timing was maintained by syncing internal clocks to the timing signal provided by the Global Positioning System. The instruments were powered by a single 20 W solar panel and a 60 Ahr car battery. The risk of instrument theft or vandalism is high in Ethiopia, and so sites were chosen largely on the basis of good security. The majority of sites included fenced compounds already with a local guard, such as schools, medical centers, government offices, and military camps. Data were stored on internal hard drives at 50 samples per second and downloaded in the field approximately every 6 months. All instruments were retrieved in October 2017, resulting in 1.5 years of data.
The requirement to place instruments close to people means anthropogenic noise is very high during the day. Anthropogenic noise typically ranges between 2 and 25 Hz and overlaps the frequency content of the detected earthquakes (1-10 Hz). This significantly affects the ability to detect earthquakes during the day. As a result, 4-5 times more earthquakes were detected at night (see Figure S1 for an example). In contrast, microseism noise levels are low, close to the new-low-noise model (Peterson, 1993), presumably as this area of Ethiopia is far from the coast (see Figure S2).

Earthquake Catalog Generation
We detect and locate seismicity using the Coalescence Microseismic Mapping (CMM) technique (Drew et al., 2013). The raw seismic data are filtered between 2 and 16 Hz, and then a characteristic function (CF) is calculated using the average in a short window (0.1 s) divided by an average in a long window (10 s). The length of the windows is chosen to maximize earthquake detections while limiting triggering due to artifacts. Peaks in the CF below a signal-to-noise ratio of 3 are removed, while those above have a Gaussian function fitted to them. The CF from all the stations is then migrated through a precalculated 3-D travel-time look-up table to generate a coalescence function. Peaks in this coalescence function indicate the origin time and hypocentral coordinates of an earthquake. To generate our catalog, we keep only those events with an overall signal-tonoise ratio greater than 3.
The 3-D travel-time look-up table used to migrate the CF contains stations embedded within a Cartesian grid (i.e., including station elevations). Topography is not explicitly included allowing earthquakes to locate above the surface of the Earth. We allow these "airquakes" to increase the number of earthquake detections in the knowledge that subsequent relocation usually moves the offending events to below the surface. We populate the 3-D travel-time grid using the 1-D velocity model for Aluto from Wilks et al. (2017). We expect the true velocity model for B-TM to be laterally heterogeneous because of the numerous faults and lateral variability in magmatism and heat flow. However, 1-D velocity models are a simplification of the 3-D velocity structure and capture the bulk of the seismic velocity variation (i.e., with depth). Some improvements to the earthquake locations can be gained by using 3-D velocity models, but these would not substantially change the results presented in this study.
Using the velocity model of Wilks et al. (2017), we detect~1,200 events with a signal-to-noise ratio greater than 3. These earthquakes form an initial catalog which we can interrogate to gain further insight into the structure of the B-TM region. In particular, we first aim to improve the earthquake locations by generating a more appropriate 1-D velocity model and refining the arrival time picks.

1-D Velocity Model Calculation
We use the program VELEST (Kissling, 1988) to simultaneously invert for the earthquake locations, a "minimum 1-D-velocity model" (i.e., a well-suited velocity model for earthquake location), and station corrections. VELEST solves the coupled earthquake location-velocity model problem using full inversion of the damped least-squares matrix. Because of the nonlinear nature of the inversion, the solution is obtained iteratively.
An important part of the inversion for a minimum 1-D velocity model is to fully explore the parameter space. This relies on trialing a large range of starting models and then converging on the best-fitting, geologically reasonable model. Most 1-D velocity models have similar features: a high-velocity gradient in the upper crust, a sharp change in the gradient at some depth, and a low velocity gradient in the lower crust. This model has three unknowns: the velocity at the surface, the velocity at the gradient change, and the depth of the gradient change. To explore the parameter space, we generate a suite of initial velocity models by varying the three identified parameters using geologically reasonable ranges. The final model, including station corrections, is that which minimizes the root-mean-square misfit with geologically reasonable earthquake locations and seismic velocities.
Calculation of the velocity first proceeds by generating a catalog of events which have a large number of accurate P and S-wave arrival time picks and are located within, or close to, the footprint of the array. We select 300 events from the CMM catalog which have a high signal-to-noise ratios (>5) and manually refine the arrival times of each event. We then select those events which have more than six manually refined arrival time picks. This leaves us with a catalog of 251 well-located events to use in VELEST ( Figure S3). All of these events are located at depths less than 9 km bsl. As such, the derived velocity model is only for the upper crust.
The Vp/Vs ratio calculated from VELEST is not smooth because VELEST performs a separate inversion for P and S-wave velocity structure ( Figure 2). We chose to use a fixed Vp/Vs ratio to calculate the S-wave structure rather than that output from VELEST. To calculate this Vp/Vs ratio, we calculate the Vp/Vs recorded by each earthquake using the method of Wadati (1933). In this technique, the Vp/Vs ratio for each event is calculated by plotting the P-wave arrival time against the time delay between the P and the S-wave arrival time for each station that records the earthquake. The gradient of this line is proportional to the Vp/Vs ratio. The mean Vp/Vs ratio is 1.78 ( Figure 3), similar to the Iceland rift (Greenfield et al., 2016) and higher than observed in regional Ethiopian studies (1.75; Keir, Ebinger, et al., 2006). Earthquakes are commonly clustered in depth around sharp changes in velocity; to avoid this, we have linearly interpolated the velocity model output from VELEST to generate the final model (red lines, Figure 2).

Comparison to Other Areas of the MER
The local, 1-D, seismic velocity model has been calculated in this study for B-TM. The only previous model covers the entire rift (Daly et al., 2008). To locate earthquakes around Aluto, Wilks et al. (2017) used a priori information from boreholes to constrain the upper 2 km of the crust and a regional model (Daly et al., 2008) for the remaining part of the crust. The model for B-TM is slower in the upper 5 km of the crust than the regional model and is~1 km/s faster at greater depths. The rift-wide model includes areas which are influenced by volcanism and others with little volcanism but a large amount of sediments. As such, the model is an average and is possibly not the best to use in local networks beneath volcanoes where high-velocity intrusive bodies are present in the crust and low-velocity ash falls and lava flows are located near the surface. The active-source model from Maguire et al. (2006)

Refined Earthquake Catalog
We refine the~1,200 CMM detected and located events by manually picking the arrival times for every event detected by CMM. Weights are assigned to each pick by considering the quality (i.e., sharpness) of the arrival and converted into a time-pick error. We use a scheme in which the best quality arrival time picks (weight = 0) have an error equal to the sampling interval (0.02 s). There are three lower weights (1-3) with time pick errors of 0.05, 0.1, and 0.5 s, respectively. Phases with a weight of 4 are not used and are assigned very high errors (9,999 s). This step removes traces which do not record the detected earthquake and negatively contribute to the CMM generated location.
The program NonLinLoc (Lomax et al., 2009) is used to relocate the manually refined events with the minimum 1-D velocity model calculate using VELEST. NonLinLoc produces optimal earthquake locations by efficiently sampling the location probability density function (PDF) built using a likelihood function. We choose to use the Equal Differential Time likelihood function (Font et al., 2004) as it is more  Geochemistry, Geophysics, Geosystems robust to the presence of outliers than the traditional least-squares, L2-norm likelihood function. The PDF is sampled using an oct-tree approach in which 3-dimensional cells are recursively divided to create a cascade of sampled cells. The density of the sampled cells follows the PDF values of the cell. The advantage of the octtree approach is that complex PDFs can be imaged much faster than using exhaustive grid-search techniques. Beneath B-TM earthquake location, PDFs were mostly unimodal with depth errors of approximately 1.5 km. Some shallow earthquakes had bimodal earthquake locations PDFs. We chose to use the deeper location as they matched the location of the surrounding, better located earthquakes. By using NonLinLoc, manually refined arrival time picks, an optimized velocity model, and P-wave station corrections, we reduce the overall root-mean-square misfit by 25%. The final catalog contains 1,012 events with an average number of picks per event approximately equal to 8.

Earthquake Magnitudes
The magnitude of the detected earthquakes (M L ) is calculated using the equation from , where A WA is the maximum zero-to-peak amplitude corrected to the response of a Wood-Anderson seismograph, C is an individual station correction, and r is the hypocentral distance in kilometers. We attempted to invert for a new local magnitude scale using the earthquakes recorded around B-TM. This would highlight any significant differences between the attenuation structure of B-TM compared to the wider rift. However, the magnitudes calculated using the new equation were identical (within the expected error of at least ±0.1 magnitude units) to the original  equation. Instead, we calculate individual station corrections for the stations used in this study to reduce the error in the magnitude and account for near-surface structure. Local magnitudes vary between 0 and 2.7, although the magnitude of completion (M c ) is approximately 1.3. The Gutenberg-Richter b-value (b) and its accompanying error (σ b ) are calculated using the following equations (Aki, 1965), ( 2) where M is the mean magnitude in the catalog, ΔM is the magnitude binning width, M i is the magnitude of earthquake, i, and n is the total number of earthquakes in the catalog. The data show a good straight line fit at lower magnitudes (Figure 4), but a change in slope is observed at local magnitudes greater than 2.2. This is probably caused by the limited length of time our network was deployed not allowing us to see the full earthquake-size distribution. The overall b-value for the B-TM region is 1.32 ± 0.07 (Figures 4, S4, and S5), and there are no significant differences when our data are split into separate clusters ( Figure S6). The b-value is significantly higher than the global average of 1 (El-Isa & Eaton, 2014) but similar to that observed nearby at Aluto (1.40 ± 0.14; Wilks et al., 2017).

FPS
Fault plane solutions (FPS) are calculated using first motion P-wave polarities picked on this catalog of welllocated events. We invert for the best-fitting FPS by following a Bayesian approach outlined in Pugh et al. (2016) and implemented in the program MTFIT (Pugh & White, 2018). The method allows for the robust inclusion of measurement and location uncertainty estimates in the final PDF. This PDF is determined over the entire moment tensor space given the observed data which can be polarity picks, amplitudes, or amplitude ratios.
We use only P-wave polarity picks to calculate FPS because the amplitudes generated from this noisy data set with small magnitude earthquakes degraded the final solutions. We also constrain the solution to be double-couple because of the low number of polarity picks (the maximum possible is 12). Tests reveal that no event requires non-double-couple components in order to explain the observed polarity distribution, but we cannot preclude the possibility that some events may have non-double-couple components. FPS are calculated for those events with more than five P-wave polarity observations, although many of the events were too poorly constrained to display. From a total of 96 events with more than five

10.1029/2018GC007648
Geochemistry, Geophysics, Geosystems polarity observations, we calculate 15 well-constrained FPS ( Figures 5 and S7). Poorly defined FPS often have more than one solution with large variability in the strike, dip, and rake parameters.
Every FPS, except two, indicates normal fault failure. Eight contain fault planes which are subparallel to the WFB and record the predicted failure in this extensional environment. Five contain faults planes at significantly oblique angles to the predicted trend and indicate variability in the orientation of the ruptured faults.

Earthquake Location Analysis
The earthquake catalog highlights three persistently active regions (Figure 4). The most active region is located between the volcanic centers of Bora and TM (red-colored circles, Figure 4). This cluster (earthquakes closely spaced in space but not necessarily with time) is elongated east-west and contains 54% of the total number of detected earthquakes. These earthquakes occur in at least eight swarms (earthquakes closely spaced in time and space; red arrows, Figure 6). To the north of this cluster, a cluster of low magnitude events is located within a broader region of distributed seismicity (cyan-colored circles, Figure 4). Earthquakes within this region are outside the footprint of the array between February 2016 and January 2017, until additional stations are deployed in February 2017 (see Figure 1 for locations). The earthquakes within this cluster are associated with Bora and also occur in swarms (cyan arrows, Figure 6). Significant swarms occur in May and August 2016, and a number of smaller swarms occur in Geochemistry, Geophysics, Geosystems January 2017. The last significant cluster of seismicity occurs beneath TM (yellow-colored circles, Figure 4). The cluster contains 110 events (11% of the total). In contrast to the other clusters, earthquakes located here do not occur within swarms but are generated continuously during the observational period (yellow-colored circles, Figure 6).
In depth, the earthquake distribution is peaked at depths between 0.5 and 1 km bsl ( Figure 5). Only 70 (7%) earthquakes are located shallower than sea level, and only 23 (2%) earthquakes are located deeper than 5 km bsl. Of the earthquakes shallower than sea level, 24 are located above the surface. These indicate poor arrival time picks and/or heterogeneities in the near surface velocity structure, which are likely to vary significantly across the region. We include them in the catalog as the earthquakes are real and their epicentral locations are likely to be fairly accurate. The high proportion (35%) of shallow (depth < 0 km bsl) earthquakes which are airquakes suggest that the depths of many shallow earthquakes are underestimated and that most seismicity in the area is deeper than 0 km bsl. Most earthquakes deeper than 5 km bsl are located near the edge of the study area, away from the volcanic centers. This suggests that a brittle-ductile transition beneath B-TM exists at a depth of between 4 and 5 km bsl (5.6-6.6 km below the surface).
The detected seismicity occurs mainly in swarms. More than 20% of the detected seismicity occurs in the early half of April 2016 ( Figure 6). Many more earthquakes than we report here occurred in this swarm. However, they are only detected by the two nearest stations and thus cannot be well located. Swarms were analyzed to see whether the activity migrates over time, possibly indicating the movement of fluids along the fault plane. We do not observe any swarm migration in the data set, but with the reported errors being large (>1 km) and swarms being confined to relatively small areas, we cannot rule out any small-scale swarm migration.
The seismicity rate over the length of the deployment is 2.0 detected earthquakes per day. During the first 12 months of the deployment, the seismicity rate is higher (2.3 earthquakes per day); although during some swarms, the rate is significantly higher than this. The seismicity rate clearly falls for the final 6 months of the deployment (Figure 6) with, on average, only 1.1 earthquakes per day. The fall in seismicity rate cannot be related to changes in the network as the number of stations increases during this period ( Figure 6). Any seasonal changes in the seismicity rate, such as those seen around Aluto due to rainfall and lake loading , are not possible to observe in this short data set.

Refined Earthquake Catalog
For further interpretation, we refine the catalog by selecting only those events that are well located. We define a well-located event as having a maximum azimuthal gap between azimuthally adjacent stations <270°, RMS misfit <1 s, and horizontal errors of less than 4 km. A total of 505 events (50%) remain after this filtering process. These limits are more generous than those normally used in earthquake studies, but they increase the number of earthquakes available for interpretation at the cost of some precision. Figure 5 shows only these events located using NonLinLoc. The removal of relatively poorly located events increases the visibility of the clusters, especially between Bora and TM. Many events were lost from the cyan cluster ( Figure 4) as they had very high maximum azimuthal gaps.
Clusters of earthquakes located between Bora and TM illuminate an intensely faulted region with at least two right-stepping, NE-SW oriented, en-echelon clusters. Focusing on this region (Figure 7) reveals significant complexity between multiple fault orientations illuminated by different swarms. The along rift cross section from A to C cuts the most seismically active region in the study area. Coloring the earthquakes in different swarms in different colors reveals discrete faults which rupture during the time period. The earthquakes indicated by the orange circles occur during August 2016 and are located beneath a volcanic ridge defining part of the Bora caldera (Figure 7). One well-constrained FPS is calculated for an earthquake during the swarm and is a thrust mechanism with a fault plane oriented subparallel to the trend of the WFB and the volcanic ridge.
On further inspection, the most active cluster in the region is clearly split into two adjacent clusters. One is oriented subparallel to the WFB and is defined by the earthquake swarm in March 2017 (light-blue circles, Figure 7). The other cluster is active during multiple swarms, which are all detected in April 2016 (dark-blue, green, and pink circles, Figure 7). While initially appearing to be trending subparallel to the WFB, the cross section reveals that this swarm defines a fault plane dipping to the northeast at an angle of~45°. Unfortunately, the events during April 2016 do not have enough observations to calculate a FPS. The April 2016 swarm is located shallower than the March 2017 swarm, and it would therefore appear that the two faults do not cross at depth.
Further north, a well-defined fault is clearly imaged in the cross section from B to C (Figure 7). The fault dips to the southeast at a low angle (~20°). By coloring earthquakes which define this plane purple, it is clear this fault is active throughout the experiment but is especially active during May 2016, September 2016, and in a number of smaller swarms during January 2017. These swarms are active for a couple of weeks, in contrast to swarms active to the south which are typically active for fewer than 3 days. No events within this cluster have good FPS, although deeper events located in the same area have extensional mechanisms with a range of

10.1029/2018GC007648
Geochemistry, Geophysics, Geosystems fault plane strikes. The surface above this cluster is heavily altered due to hydrothermal processes (orange regions, Figure 7) suggesting that seismicity may be reflecting the movement of hot fluids in this area.

Frequency Content
The frequency content of earthquakes can be useful in interpreting the mechanisms driving seismicity in an area and is influenced by the earthquake source mechanism and the properties of the area imaged by the earthquake. Around volcanoes, a large variation in the frequency content of earthquakes can be observed (see McNutt, 2005, for an overview). Typically, these are split into discrete categories: volcano-tectonic (VT), long-period, and hybrid events, but in reality, the variation is a continuum. VT events are traditional earthquakes which contain high frequencies and typically have impulsive onsets. In contrast, long-period events contain predominantly low frequencies (LFs; <5 Hz) and often have emergent arrivals which can be difficult to pick accurately making them difficult events to locate. Hybrid events often have high frequency, impulsive onsets but have a long, LF trailing coda.
We analyze the frequency content of earthquakes within our NonLinLoc catalog by calculating a Frequency Index (FI) for each event (Buurman & West, 2010). The FI expresses the ratio between the amount of energy in a LF band (0.6-1.2 Hz) and that in a high frequency band (6-12 Hz). A logarithm (base 10) is used to rescale the difference in FI and aid plotting and interpretation, where A h and A l are the mean amplitudes in the high frequency and LF windows, respectively. The FI will vary with distance from the source depending on the amount of attenuation between source and receiver. By considering the effect of attenuation on (3), it can be seen that where FI 0 is the FI at the source (zero offset), Q is the seismic quality factor, r is the distance, and v is the propagation velocity. ω h and ω l are the mean angular frequencies in the higher and lower frequency windows, respectively. Assuming Q to be constant, (4) shows that, to first order, FI will decrease monotonically with increasing hypocentral distance.
The top panel of Figure 8 shows how the FI for each observation varies with raypath length. A clear trend is observed, but it is not a straight line (as predicted by equation (4)). To compare the FI 0 between events, the dependence of FI with distance must be removed. We remove this dependence by considering the change in FI through a simple 1-D attenuation model, where the symbols are as in (4) except that the FI for event, i, calculated at station, k, includes the effect of attenuation through layer j. We also include a station correction term, C k . Equation (5) is of the form Ax = b

10.1029/2018GC007648
Geochemistry, Geophysics, Geosystems and can be solved for the FI 0 of every event, station corrections, and a depth-dependent Q model. We use a least-squares approach to solve (5) and include a requirement that the Q model is smooth by adding a regularization term, where L is a regularization matrix of the form, and λ is a parameter (Lagrange multiplier) which can be adjusted. L is formed using a second-order finite-difference operator in the bottom corner, so only the Q model is constrained to be smooth. We also require Q to be positive and that the sum of the station corrections is equal to zero. This avoids any numerical instabilities during the inversion.
A range of layer depths for the Q model are trialed, calculating the distance each ray spends in each Q layer in addition to the average shear-wave velocity. A Q model with 1-km-thick layers to a depth of 5 km and a half- Figure 10. Examples of a volcano-tectonic (left panels) and low frequency (right panels) events. Top panels show waveforms recorded on the North component at station ANOL (see Figure 1 for location). The amplitude spectra plotted on a log-log plot and log-linear plot are displayed in the middle and bottom panels, respectively. FI = Frequency Index.

10.1029/2018GC007648
Geochemistry, Geophysics, Geosystems space below this produced the best fit to the data. The final Q model is shown in Figure 8. Q decreases from a surface value of 280 to a minimum of 115 at a depth of 2 km bsl.
After a correction for Q, most events have a FI between 1 and 2 ( Figure 9). There are, however, a small number of events with an FI as low as À4. These contain much more LF energy. This feature is clear in the seismic waveforms of the seismicity ( Figure 10) and is characterized by peak frequencies between 0.6 and 1 Hz ( Figure 10). The LF nature of this seismicity cannot be a path effect as VT earthquakes with high FI values are located in similar locations and the LF earthquakes have a similar magnitude to the VT seismicity. The arrivals of P and S-waves are impulsive (Figures 10 and 11), making the seismicity relatively easy to locate with small errors similar to VT seismicity. In fact, because the LF energy is at longer periods than the anthropogenic noise, the LF events are detected at stations farther away than the VT seismicity. Most of the LF events locate beneath TM (red circles, Figure 5), suggesting the source of this unusual seismicity is associated with TM.

Seismic Anisotropy
We investigate the seismic anisotropy of the B-TM volcanic complex by evaluating SWS observed in events from the refined NonLinLoc catalog. In an anisotropic medium, an incident shear wave is split into two orthogonal components, with one component propagating faster than the other. SWS measurements are made by calculating the polarization of the faster component (ϕ) and the delay time between the fast and slow component (δt).
Seismic anisotropy reflects structural order within an elastic body and leads to variations in seismic velocities with direction. SWS is arguably the most unique indicator of seismic anisotropy, but there are many mechanisms for anisotropy. Deformation can lead to the lattice-or crystal-preferred alignment of minerals (e.g., olivine in the upper mantle; Mainprice et al., 2013) or phyllosilicates in sedimentary rocks (e.g., Valcke et al., 2006). On a larger scale, the rhythmical layering of structures with contrasting elastic properties also leads to directional variations in seismic velocities (e.g., Backus, 1962). In the shallow crust, the alignment of cracks or fractures parallel to the direction of maximum horizontal compressive stress is thought to be a primary mechanism for anisotropy. Here the fast shear wave is polarized parallel to the orientation of the cracks or fractures (e.g., . The style of anisotropy is sensitive to the degree of fracture connectivity and nature of the fluids (see, e.g., Baird et al., 2018). Regardless of cause, the aligned structures must have lengthscales much less than the seismic wavelength. To calculate δt and ϕ from the raw waveforms, we follow the technique of Silver and Chan (1991). The shear-wave arrival is windowed, and the particle motion of the two horizontal components is analyzed. If anisotropy is present, the S-wave particle motion will be elliptical or cruciform (although at these near offsets a cruciform particle motion is unlikely). A grid search is performed over all possible values of δt and ϕ, where the two horizontal waveforms are rotated by ϕ and shifted by δt. The best-fitting values are those that minimize the second eigenvalue of the corrected particle motion covariance matrix. This signifies that the particle motion is approximately linear. This analysis also determines the source polarization of the shear wave before propagating through the anisotropic medium. Errors in δt and ϕ are calculated using an F-test. This technique is highly dependent on the window chosen, as well as other In each panel, the upper-left plot shows the filtered waveforms recorded on the east (e), north (n), and vertical (z) components. The S-wave arrival time pick is indicated by the thick black line, and the window used is shown by the gray box. In the upper-right plot, the upper two waveforms show the two horizontal components rotated into the incoming polarization direction (p) and its perpendicular value (p⊥). The lower two waveforms show the horizontal components rotated into the same reference frame but after correction by δt. The lower-left plot shows the two rotated horizontal components (upper plots) within the final window chosen and the particle motion (lower plots). Before (left) and after (right) correction is displayed. In the lower-right plot, the contours of the smallest eigenvalue of the covariance matrix are shown with the final chosen shear-wave splitting measurement indicated by the blue cross.

10.1029/2018GC007648
Geochemistry, Geophysics, Geosystems parameters such as filtering applied to the data. To obtain the highest quality results, we use a range of window lengths and pick the most stable results using cluster analysis (Teanby et al., 2004). In addition, a range of filters is used, and the "best" result picked (Savage et al., 2010).
We use our catalog of manually picked events as input for SWS analysis. This contains the high-quality S-wave picks required for analyzing SWS. We choose only those earthquakes which lie within a shear-wave window , thus avoiding effects caused by S-to-P conversions at the free surface. We choose an angle of 55°from vertical, larger than the critical angle that would typically be used (sin À1 V S /V P ≈ 35 ∘ ). However, because of low, near-surface velocities, rays are steeper at the station than the straight-line path used to calculate the shear-wave window.
We manually grade the SWS observations into one of three categories: A (best), B (acceptable), and X (not used). An example of each grading is shown in Figure 12. We grade observations based on error in δt and ϕ, linearity of the particle motion after correction, stability of the clustered solutions for δt and ϕ, and minimal energy on the component transverse to the initial source polarization after correction. From a total of 1,170 observations, 125 are graded A (10%), 214 are graded B (18%), and the remaining observations are graded X. It is worth noting that within the X-graded observations are those where no splitting is observed. Such null observations could be due to wave propagation through an isotropic medium, but in this environment, this is more likely due to the initial source polarization being aligned with the fast or slow shear-wave polarization. In reality though, distinguishing nulls due to alignment or due to noise is difficult, and we do not interpret the nulls in this study.
Results from the SWS analysis are displayed in Figure 13 as bars plotted at the midpoint of the raypath between source and receiver. The length of the bar is proportional to the amount of anisotropy, while the direction of the bar indicates the orientation of the fast direction. Most measurements extend in a band between TM and Bora in a southeast-northwest orientation. Another cluster of measurements is made to the north of this band. We normalize δt by the distance traveled along the raypath using where t slow and t fast are the S-wave travel times of the slow and fast shear waves, respectively, and s, a, and v are the distance along the raypath, fractional anisotropy (varying between 0 and 1), and the average shearwave velocity, respectively (Thomas & Kendall, 2002). Rearranging (6), where only the positive solution is valid for anisotropic strength to vary between 0 and 1.
We use the hypocentral distance as the length of the raypath, which better approximates the length of the raypath length than the epicentral distance for these earthquakes. A shear-wave velocity of 2.24 km/s is used which is the S-wave velocity in the earthquake nucleation region. We observe up to 10% anisotropy, although the bulk of the measurements lies between 2% and 6% ( Figure 13). Any depth dependence of the anisotropy is difficult to analyze because of the limited depth distribution of the seismicity and high uncertainty in earthquake depth. Nowacki et al. (2018) observed that the bulk of the observed anisotropy around Aluto is concentrated into the top 3 km. Initially, the results of this study appear similar. In Nowacki et al. (2018), the apparent strength of anisotropy decreases progressively at depths greater than the actual anisotropic layer. This is because the splitting magnitude is normalized over the entire raypath, even if it transits deeper isotropic regions. In contrast, Figure 13 shows little correlation between anisotropy magnitude and depth. A few low anisotropy observations are from the deepest events, but some of the largest are also from the deepest events. A broader depth distribution of seismicity beneath B-TM would be required to analyze any variations in the strength of anisotropy with depth. Geochemistry, Geophysics, Geosystems

Results and Interpretation
The overall distribution of SWS fast directions has a broad peak between 015°and 045°similar to the distribution of fracture trends observed around B-TM of 015-030°( Figure 14). However, observations at individual stations show significantly more variation (Figure 13). We can interpret the SWS fast-direction orientations by considering the cause of this variability. Within the upper crust, most seismic anisotropy is caused by the alignment of microcracks at wavelengths smaller than the seismic wavelength . This orientation of these cracks is the result of an interplay between stresses due to regional tectonics, local VT stresses, and any preexisting structures. Around B-TM, distinguishing between these mechanisms is difficult without a larger data set.
The four stations located in the southwestern corner of Figure 13 show a high degree of similarity in the observed ϕ distributions. The dominant direction at each station is NE-SW, similar to the orientation of the WFB. Stations in the east of the study area show a more N-S-oriented ϕ distribution. This is possibly the result of a subtle change in the orientation of the faults within the WFB around TM (Figure 13). This is most clear at the station ANOL (see Figure 13 for location) where individual measurements (mostly graded B) are parallel to the nearby faults exposed at the surface. The change in the WFB orientation is likely caused by the local stresses around TM.
The one station in the north of the area (GEBI, Figure 13) contains a very different ϕ distribution, with a dominant E-W direction. The earthquakes from which these observations are derived are located within the northern cluster of earthquakes (Figures 5 and 13) and are different to the earthquakes which are used at the other stations in the area. This could be the reason why few E-W ϕ are observed at any other station. It is interesting to note though that the few observations from these earthquakes recorded at the station TULL (circled red, Figure 13) are oriented NW-SE, significantly different to the dominant ϕ at TULL (NE-SW). This suggests that ϕ at TULL is dependent on back-azimuth, and therefore, the anisotropy is spatially heterogenous.
SWS observations around B-TM reveal the complex stress pattern induced around volcanoes due to local and regional stresses. The dominant orientation of ϕ is subparallel to the rift indicating that the regional stress exert a strong control on ϕ. However, the variability in ϕ (particularly at TULL and GEBI) shows that local stresses around the volcano significantly modify ϕ. This is similar to that observed at Aluto where Nowacki et al. (2018) matched the observed ϕ on the Aluto edifice to two competing fracture orientations. In contrast, ϕ beneath Corbetti is dominated by preexisting fault structures (Lloyd, Biggs, Birhanu, et al., 2018), although the seismic network and earthquakes used to generate these observations limit the spatial coverage of ϕ observations. Around B-TM, we suffer from similar difficulties, making quantitative analysis of the contribution of different fracture orientations impossible.

B-TM Seismicity
Seismicity around B-TM is closely associated with the expression of volcanism, being centered around TM, rather than the tectonically generated features within the WFB. This result is to be expected as strain rates measured around volcanoes (cm/yr) are an order of magnitude higher than tectonic extension (5-6 mm/yr; Bendick et al., 2006). The regional extensional stress regime has the effect of biasing the style of faulting to be normal, but local stresses induced by the volcano causes some the earthquakes to exhibit compressional style faulting.
We use the depth distribution of the seismicity to constrain the brittle-ductile transition. Beneath B-TM, we observe few earthquakes deeper than 5 km bsl. If this depth represents a brittle-ductile transition, then a high geothermal gradient (~60°C/km) is required, indicating the presence of a heat source at depth. Given that the seismicity is predominantly located near TM, we suggest that the heat source is located beneath TM rather than Bora. Such a heat source is required in order to drive the hydrothermal circulation observed around B-TM. The fluids involved in the hydrothermal circulation could be sourced directly from the magmatic body degassing or indirectly via meteoric water interacting with the elevated geothermal gradient. Geochemical evidence from other volcanoes in the MER shows that the bulk of the water found in hydrothermal fluids is sourced from rainwater (Darling et al., 1996), and we suggest the same is likely here.
Fluid circulation can cause locally high-pore-fluid pressures which act to lower the normal stress. This causes the triggering of earthquake swarms with high b-values (Farrell et al., 2009;Sholtz, 1968;Wyss, 1973) and at low angles to the maximum compressive stress (Collettini & Holdsworth, 2004). Around B-TM, seismicity predominantly occurs within swarms, and the overall b-value is much higher than the average global value. In addition, if we assume the purple cluster in Figure 7 to be an extensional fault, then high-pore-fluid pressures could be the reason for the unusually low dip. Fluid circulation therefore plays an important role in determining where seismicity occurs beneath B-TM and thus the seismicity can be used to image fluid pathways through a network of fractures in the crust. These seismically imaged fluid pathways could be targets for possible geothermal power generation.
Given the low amount of seismic moment released during the swarms of earthquakes beneath B-TM, we suggest that rather than new fractures being formed, preexisting fractures are utilized. We might expect the structures related to the WFB to dominate the fracture network within the crust. Indeed, the fast shear-wave orientations are broadly parallel to the WFB trend. However, the cluster orientation (Figure 7), FPS ( Figure 5), and details of the SWS observations ( Figure 13) reveal a significant amount of variability in the orientation of fractures at depth. Such fractures are either not dense enough to significantly alter the fast shear-wave orientation or reveal that SWS is caused by stresses only. This is different to Aluto, where interacting NE-SW WFB fractures and E-W trending caldera structures cause fast shear-wave directions to be~060° .
A possible origin for the variability in the preexisting fault network orientation is structures associated with the Precambrian crust underlying the MER. N-S rifting occurred during the Mesozoic (Corti, 2009), and the cross-rift structures produced may have been reactivated as recently as the Pleistocene (Corti et al., 2018). Cross-rift structures are particularly evident beneath Corbetti (but also elsewhere in the East African Rift system; e.g., Robertson et al., 2016) where a significant preexisting E-W-oriented structure causes: east-west elongation of the caldera, differences in uplift between the north and south, prominent changes in the resistivity between north and south and, east-west oriented fast shear-wave directions (Lloyd, Biggs, Birhanu, et al., 2018). Around B-TM, there are no mapped cross-rift structures (Corti, 2009;Corti et al., 2018), suggesting that any regional ancient fabrics are either not present or not strong enough to cause the local variability fracture orientation.
Instead, we suggest that the variability is caused by structures related to the formation of a number of calderas in the region. At least two but possibly three calderas are mapped around B-TM (Figure 1). These are undated but could be related to a series of large caldera forming eruptions across this region of the MER between 180 and 320 ka (Hutchison, Fusillo, et al., 2016). Both Bora and TM host calderas which overlap in the region where most seismicity is recorded. As mapped in Figure 1, the TM caldera is the most recent as it cuts the Bora caldera and is less infilled. The southern cluster of earthquakes in Figure 7 could be related to hydrothermal fluid flow inducing failure on caldera faults at this location. The swarms, active during April 2016 (dark-blue-, green-, and pink-colored circles, Figure 7), dip to the northeast, parallel to the TM caldera at this location. In contrast, the other highlighted swarms (indicated by the orange and light-blue circles, Figure 7, and active August 2016 and April 2017, respectively) are parallel to structures associated with Bora.

Comparison to Other Volcanoes in the MER
Only two other volcanoes in the MER have been studied using dense seismic networks: Corbetti (Lloyd, Biggs, Birhanu, et al., 2018) and Aluto (Wilks et al., 2017). All three volcanoes have different seismicity patterns and eruptive activity (Fontijn et al., 2018). This may reflect structural differences between the volcanoes or, given the short observation periods, observations of the volcanoes during different periods of their respective eruptive cycles.
Spatially, there are large differences between the seismicity distributions of Aluto, Corbetti, and B-TM. Seismicity around B-TM is clustered onto faults away from the center of the volcanoes, while around Aluto, the seismicity is not clustered but is at its densest at the center of the volcano. The distribution around Corbetti is difficult to interpret because of the high magnitude of completion but seems to be clustered in the center of the caldera at depths above the inflating deformation source (Lloyd, Biggs, Birhanu, et al., 2018). Aluto, with shallow, distributed seismicity, has a shallow hydrothermal reservoir and exploits a distributed fracture network. In contrast, the hydrothermal reservoir at B-TM is deeper and exploits larger fluid pathways caused by preexisting faults. Corbetti has SWS dominated by preexisting E-W-oriented structures (Lloyd, Biggs, Birhanu, et al., 2018). This structure is illuminated by deep seismicity, indicating that it acts as a highway for fluids to migrate through the crust. Such deep, preexisting, cross-rift structures are not seen beneath Aluto or B-TM and must be less important for controlling the overall structure of the rift.
The seismicity rate is difficult to compare between B-TM, Aluto, and Corbetti because of large differences between the density and number of seismic stations deployed. In particular, the deployment around Corbetti (Lloyd, Biggs, Birhanu, et al., 2018;Wilks, 2016) only consisted of a maximum of five stations deployed at any one time. We use the parameter a, calculated from the fit to the Gutenberg-Richter law (section 4.4) to compare the overall seismicity rate. a indicates the number of earthquakes we would expect to find in the region above a magnitude of zero per year. Aluto, Corbetti, and B-TM have a-values of 4.41, 4.35, and 4.18, respectively (Wilks, 2016;Wilks et al., 2017) indicating that the three volcanoes have similar levels of seismicity.

LF Earthquakes
LF events around volcanoes in the MER have not been observed before but are frequently observed around many volcanoes around the world (Chouet & Matoza, 2013;McNutt, 2005). LF events are characterized by LF waveforms (typically <2 Hz), long, "ringing" codas, and emergent phase arrivals (McNutt, 2005). They typically occur at shallow depths (often <1 km below the surface) and within swarms where event waveforms are highly similar (e.g., Bell et al., 2017). Models to explain LF seismicity include resonance of a fluid-filled cavity following a brittle-failure "kick" (e.g., Chouet, 1996;Neuberg et al., 2006) and seismogenic failure of unconsolidated sediments near the surface and wave propagation through a shallow, low-velocity layer (Bean et al., 2013). LF earthquakes around TM are significantly deeper than 1 km below the surface, so neither of these mechanisms can explain the observed seismicity.
Deeper LF events, often with impulsive phase arrivals, can be found beneath some volcanoes round the world (e.g., Greenfield & White, 2015;Koyanagi et al., 1987;Nakamichi et al., 2003). These are often deeper than 10 km below the surface and are associated with the movement of magma through the crust. Around TM, where most of the detected LF events are located in this study, there is little evidence for magma movement within the crust. We observe no migrating swarms of earthquakes associated with dike or sill intrusions near the surface, and there is no surface deformation observed at TM and therefore conclude that the LF earthquakes are not caused directly by magma movement.
There is, however, lots of evidence for the movement of other fluids (e.g., water) beneath TM. In laboratory experiments, rapid movement of fluids have been shown to emit acoustic waves similar to LF events (Benson et al., 2008). However, the LF events generated in these experiments are either emergent or hybrid style events with high-frequency onsets. Such events do not match our observations. The movement of large volumes of fluids also occurs on subduction zone interfaces and in the accretionary prism. In these settings LF earthquakes and tremor are observed and have been proposed to be caused by brittle failure on faults with low confining pressure, probably due to high (near lithostatic) pore fluid pressures and low rupture velocities (Peng & Gomberg, 2010). A similar model, with hydrothermal fluids locally causing high pore fluid pressures, better explains the origin of the LF events beneath TM.

Conclusions
We have detected and located more than 1,000 earthquakes in the B-TM volcanic region ranging in local magnitude between À0.08 and 2.72. The magnitude of completion for our deployed network is 1.3, and the Gutenberg-Richter b-value is 1.32 ± 0.07, indicative of low effective stresses which we interpret to be caused by high pore fluid pressures. Earthquakes around the B-TM volcanic center in the MER are highly clustered and occur within swarms. Swarms are likely triggered by the movement of hydrothermal fluids which increase the pore fluid pressure and cause brittle failure. Hydrothermal fluids are likely heated by a shallow heat source probably located at depths greater than 5 km bsl beneath TM.
The obvious structural orientation in the area, the WFB, provides only weak controls on the orientation of the seismogenic fractures. Instead, we propose that preexisting fractures, most likely related to a number of volcanic calderas, form the fracture network exploited by hydrothermal fluids. Such structures (or the stress that causes them) locally alter the SWS orientation and the strike of fault planes away from the expected rift parallel orientation. Around B-TM, a region where multiple calderas overlap is the most seismically active area because of significant hydrothermal fluid circulation. This region would be an ideal target for geothermal power generation.
Analysis of the frequency content of the earthquakes reveals the presence of LF events. These events occur in swarms and are well-located due to impulsive and easily identifiable P and S wave arrivals. The LF earthquakes colocate with VT events beneath TM. We propose that exceptionally high pore-fluid pressures cause slow rupture speeds and the observed LF events. would not have succeeded. Seismic data are stored at SEIS-UK and although currently embargoed, will be available, open-access, on IRIS (https://www.iris. edu/hq/) after October 2020. Data sets generated from this study, such as earthquake hypocenters, can be downloaded from the supporting information. Generic Mapping Tools (GMT; Wessel et al., 2013) was used to make some of the figures, and Obspy (Beyreuther et al., 2010) was used extensively for data analysis. Earthquake spectra were calculated using a multitaper spectrum analysis library (Prieto et al., 2009). The use of CMM and computing resources at the Bullard Laboratories was generously supplied by Bob White. We thank two anonymous reviewers and an anonymous associate editor for their thoughtful comments, which improved the manuscript.