Rifting of the Kalahari Craton Through Botswana? New Seismic Evidence

Botswana is a country with relatively low seismic activity that experienced an unexpected Mw 6.5 earthquake on 3 April 2017. Using data from the first countrywide seismic network, we established a Botswana earthquake catalog for the period January 2014 to February 2018. Two areas of elevated seismic activity were detected. The first one is the Okavango Rift Zone in northern Botswana, an area that is known to be active. The other one is associated with the 2017 mainshock and its aftershocks in central Botswana; it follows the Paleoproterozoic suture between the Limpopo Belt and the Kaapvaal Craton. Double‐difference relocation of these events revealed a reactivated fault system with a northwesterly strike with the aftershocks occurring at shallower depths than the mainshock at 29 km. The focal mechanisms of the mainshock and selected aftershocks are of normal faulting type with strikes similar to the orientation of the fault system. The unidirectional rupture of the mainshock in the lower crust combined with the upward migration of the aftershocks along the Moiyabana Fault and a thin low velocity anomaly the uppermost mantle are consistent with the events being produced by eclogitization of a dry metastable granulite facies rock by fluid intrusion from the mantle in an extensional stress regime. The Okavango Rift Zone is generally interpreted as the terminus of the southwestern branch of the East African Rift System. We suggest that the recent earthquakes in central Botswana may be considered as the southern continuation of this branch.

3 of 17 the Kaapvaal, Zimbabwe and Congo cratons down to depths of 200-250 km, whereas low upper mantle shear velocities (4.4-4.5 km/s) were obtained beneath the northeastern part of the Okavango Rift Zone which are likely associated with the southwestern branch of the East African Rift System White-Gaynor et al., 2021).
In this study we investigate the Botswana seismicity by analysis of the continuous data of the NARS-Botswana network. We then improve the (relative) locations of the 2017 mainshock and its aftershocks using the double-difference method and determine the focal mechanisms of the seven largest events as well as the rupture propagation of the mainshock to better understand the nature of this intraplate seismicity.

Event Detection and Localization
We used the NARS-Botswana data for the period from 1 January 2014 to 28 February 2018 to obtain a Botswana event catalog that improves on other catalogs (which are mainly based on data from stations outside the country) both in terms of detection as well as localization of the epicenters. Because there are various mines in the country (Figure 2), not only tectonic earthquakes will be detected but also mining related events.
We employed the SeisComP3 software package (GFZ, 2008) with automatic and manually picked P-wave arrival times to detect and locate seismic events. The procedure is summarized as follows. First, a simple STA/LTA (short-term average/long-term average) detector was applied to the 0.7-2.0 Hz filtered data to obtain rough estimations of potential first P-wave arrivals. After event identification, that is, recognizing that multiple P-wave arrivals are associated to a single event, an AIC P-phase picker, based on the the Akaike Information Criterion (St-Onge, 2011), was applied to the 0.5-5 Hz filtered data to achieve more accurate arrival times. Finally, all automatically determined P-wave arrivals were manually reviewed and modified if necessary. To ensure that no seismic events were missed, the detected events by SeisComP3 were compared to the catalog of unreviewed events published by the ISC (ISC, 2018). Events that were not identified by SeisComP3 but that were listed by the ISC were used as potential additional events for which the data were analyzed. If sufficiently clear, new P-wave arrivals were manually picked for those events as well. In the last step, the picked P-wave onset times were used to locate the seismic events using the reference model iasp91 (Kennett & Engdahl, 1991). Events were located if at least six P-wave picks were available with a maximum individual residual of 7.0 s and a final root mean square (RMS) residual less than 3.5 s. The procedure, including all parameter settings, is described in detail in Micallef (2019). The final catalog consists of 376 Botswana events and includes magnitude estimates that are derived from the high-frequency P-wave amplitudes. The magnitude of completeness was estimated from the Gutenberg-Richter relationship for events prior to the 2017 mainshock (Micallef, 2019), giving a value of 3.1. The catalog is provided in Table T1 of Supporting Information S1.
A map of the events is presented in Figure 2. It shows two main regions of natural seismic activity: one in the Okavango Rift Zone in northern Botswana with magnitudes ranging from 3.0 to 4.9, and the other in central Botswana which includes the 3 April 2017, mainshock and its aftershocks. Note that many other seismic events are related to mining.
Within the Okavango Rift Zone, seismic activity was highest in the eastern part of the region for the period of investigation. This is also observed in the seismotectonic map by Meghraoui and the IGCP-601 Working Group (2016) who re-evaluated historical and instrumental earthquake catalogs of the African continent. The correspondence indicates that the elevated seismicity in the eastern Okavango Rift Zone is a real feature, not biased by the limited time span of our data. It is broadly believed that seismicity in the Okavango Rift Zone is caused by the development of an active rift, that is, to the continuation of the southwestern branch of the East African Rift System which enters Botswana in the northeast (Kinabo et al., 2007;Modisi et al., 2000;Scholz et al., 1976). Moreover, it has been suggested that the earthquakes in the Okavango Rift Zone are facilitated by ascending fluids from a hot mantle along lower-crustal faults. This is inferred from the thermal structure derived from aeromagnetic and gravity data (Leseane et al., 2015) as well as from the seismic structure obtained from surface wave tomography . It is likely that the process of rifting and upward fluid transport from the mantle is further developed in the northeastern and eastern part of the Okavango Rift Zone (i.e., the terminus of the southwestern branch of the East African Rift System), which would explain the dominant seismicity at the eastern parts of the reactivated faults.
The other region of natural seismicity includes the 3 April 2017, Mw 6.5 earthquake in central Botswana at 120 km distance from the village of Moiyabana (Kolawole et al., 2017). The epicenter, located by our data at 22.636°S, 25.206°E, matches the location given by the International Seismological Center (ISC) within the uncertainties. It was the largest earthquake recorded by the network and was completely unexpected because it occurred in a region without major tectonic activity in the last ∼2 Ga (Kolawole et al., 2017). Previous studies, not using NARS-Botswana data, showed that it was a normal faulting event along a northwest striking fault (Gardonio et al., 2018;Materna et al., 2019;Midzi et al., 2018), although the studies do not agree on the direction of fault dip.
The main event was followed by 216 aftershocks (up to 4 February 2018) which are clustered along a NW-SE trending zone (Figure 2). This matches the aftershock distribution for the month of April 2017 measured by a temporary network installed immediately after the mainshock (Midzi et al., 2018). Yet, our longer-term earthquake distribution extends further to the northwest, and clearly shows that the events cluster around the boundary between the Kaapvaal Craton and the Limpopo Belt (Figures 2 and 4b).
We carefully checked our catalog for the presence of potential foreshocks because the main event occurred in a remote area without any seismicity that was reliably detected. Two events, with magnitudes between 2 and 3, occurred at distances of ∼40 km northwest of the main event (within the same zone as the aftershock seismicity) on 27 October 2014, and 10 March 2015. Although they are not true foreshocks in the sense that they occurred just before the main event, they do provide evidence for seismic activity along the boundary between the Kaapvaal Craton and the Limpopo Belt prior to the main event. Notably, even after careful data inspection, we Soda ash mine did not detect any of the 20 foreshocks with M > 3 identified by Gardonio et al. (2018), who applied a template matching technique to data recorded at distances larger than 1200 km from the main event.
On 12 August 2017, an earthquake with an estimated magnitude of 5.6 occurred at 23.626°S, 25.678°E. It is located within the Kaapvaal Craton, 120 km southeast of the main event, at the junction of the WSW-ENE trending Zoetfontein Fault with a southern branch and a northwesterly branch extending to central Botswana. The area is known to be seismically active, but mostly with magnitudes smaller than 5 (Meghraoui & the IGCP-601 Working Group, 2016). Several small events in our catalog can be associated to these known faults in the Zoetfontein Fault region of Kaapvaal Craton ( Figure 2).
In addition to the tectonic events, a few clusters of mining events were detected ( Figure 2). Most of them are from explosions in diamond mining areas (around 21.3°S, 25.5°E and 24.5°S, 24.7°E) and a few are related to coal mining (22.5°S, 27.0°E). Albano et al. (2017) verified that the 3 April 2017, main event was not caused by mining.

Relocation of the Aftershocks of the 3 April 2017, Earthquake
Aftershock distributions are useful to identify the fault plane and to determine the evolution of aftershock related fault slip. In an analysis that was carried out independently from the localization study presented in the previous section, we relocated the aftershocks of the 3 April 2017, earthquake with the double-difference method by Waldhauser and Ellsworth (2000). The double-difference method is based on the notion that the ray paths of two nearby events to the same station are similar along most of their ray paths if the hypocentral distance between the two events is small compared to (a) the event-station distance and (b) the scale length of velocity heterogeneity. If pairs of events meet these conditions, then their travel time differences observed at a common station can be attributed to the spatial offset between the events. This allows very accurate relocation of events with respect to each other if travel time differences of multiple stations with a good azimuthal distribution around the earthquake are used.
We relocated the 3 April 2017, mainshock together with 79 aftershocks with M ≥ 2.5 (until 9 November 2017) using the original locations from the ISC and USGS catalogs. P-wave arrival times were manually picked on seismograms of the NARS-Botwana network as well as of station LBTB of the Global Telemetered Seismograph Network. Relocation was done using the velocity structure at the location of the mainshock obtained from the tomographic model by Fadel et al. (2020). Tests were carried out on subsets of the data using various parameter settings and solving the system by either singular value decomposition (SVD) or the conjugate gradient algorithm LSQR (Paige & Saunders, 1982). Because of the number of events, the inversion of the entire data set was performed using LSQR and a full description of the parameter settings is given in Bouwman (2019).
Of the 80 events in total, 57 events were relocated in the area around the mainshock. Two events, the event in the Kaapvaal Craton on 12 August 2017, and an aftershock on the same day, were relocated as a separate cluster 113 km southeast from the mainshock. The other events could not be relocated because their initial locations were too far away from these two clusters. The relocated events were found to have mean two-sigma standard errors of ∼4.5 km in the three orthogonal directions (2Δx = 4.4 km, 2Δy = 4.6 km, 2Δz = 4.4 km) with depth errors (Δz) ranging between 0.4 and 1.2 km. More accurate relative locations might be obtained by including S-wave arrivals, although they are harder to pick.
The results, presented in Figure 3 and as Table T2 in Supporting Information S1, show that the dispersed distribution of hypocenters around the mainshock becomes a 24-km long zone of aligned aftershocks. Notably, the aftershocks are all located northwest of the mainshock and decrease in depth in the northwest direction. The cross-section of Figure 3 shows that the hypocenters are relocated along a plane with a dip of ∼70° toward the northeast. We observe that the main event is relocated at significantly shallower depth (18 km) compared to the original USGS location (29 km). This is caused by the fact that the double-difference method is a relative relocation method, keeping the center of the event distribution in place. With most of the aftershocks originally being located at the default depth of 10 km, whereas their actual depth is probably deeper (see next section), the relocated events are likely to be shifted upward compared to reality.

Focal Mechanisms and Hypocenter Depths
The previous localization and relocalization studies only used first arrival times. On the other hand, P and S waveforms (including surface reflected phases such as pP and sS) constrain focal mechanism, depth and magnitude of earthquakes. We modeled the P and S waveforms with the "cut and paste" (CAP) method of Zhao and Helmberger (1994) (see also W. Chen et al., 2015;Zhu & Helmberger, 1996). Green's functions were calculated with the Thompson-Haskell propagator method, and the 1-D model to calculate the synthetics is based on the average Botswana mantle model by Fadel et al. (2020) combined with the crustal structure at the location of the mainshock from the same study. A grid search is performed to find the optimum source parameters (focal depth, moment magnitude and strike, dip, rake for a double-couple point source) using the L2 misfit between data and synthetics. The method allows for travel time shifts caused by 3-D heterogeneity that are not accounted for by the 1-D velocity model. The P waveforms were filtered between 0.05 and 0.3 Hz and the S waveforms between 0.05 and 0.2 Hz to minimize the influence of anthropogenic and microseismic noise.
We applied the method to events with M > 4 which had sufficient signal-to-noise ratio: the 3 April 2017, main event, five aftershocks in its immediate vicinity and the 12 August 2017, earthquake at the Zoetfontein Fault. Figure 5 illustrates the waveform fits for the 4 July aftershock. Figure 6 presents all results and they are also listed in Table 1.
The 3 April 2017, mainshock and its aftershocks are all normal faulting events with a NW-SE strike, a strike that is consistent with the aftershock distribution (Sections 3 and 4) and the orientation of the boundary between the Kaapvaal Craton and Limpopo Belt. The event at the Zoetfontein Fault has a similar mechanism, suggestive of extension perpendicular to NW-SE oriented faults in the entire area. Furthermore, we found that all events occurred at depths between 20 and 30 km, with the main event at 29 km being the deepest. The focal mechanism and event depth of the mainshock are similar to previous studies (Gardonio et    Stations are listed with increasing epicentral distance. P waves at distances larger than 300 km were excluded because of insufficient signal-to-noise ratio. Red crosses in the focal mechanism indicate azimuth and take-off angle of the P waves on the focal sphere.  Materna et al., 2019) and suggest an initial stress release during the mainshock at around 30 km followed by triggered aftershocks that occur at shallower depths.

Rupture Propagation of Mainshock
We found that the waveforms of the mainshock were less well fit than those of the aftershocks for which modeling by a simple point source was adequate. This suggests that there are effects of rupture propagation in the waveform data of the main event which may be retrieved with the empirical Green's function approach (e.g., Hutchings & Viegas, 2012). In this method a nearby event is selected with a similar moment tensor, but with a size that is at least one magnitude lower. This lower magnitude event has similar radiation and propagation effects, but is better approximated by a point source. A recording of this smaller event is referred to as the empirical Green's function (EGF) for the station considered. Deconvolution of the main event by its EGF yields an estimate of the main event's source-time function (STF) for that station. Similar to the Doppler effect, the STF is shortest when the rupture propagates in the direction of the receiver and is longest when it propagates away from it. From the azimuthal distribution of the STF durations, the rupture direction, its extent in the horizontal plane and the rupture velocity can be estimated (Savage, 1965).
The aftershock on 4 July (2017-07-04, 11:37:06) was selected as the EGF event. A stabilized deconvolution of the vertical component of mainshock by its EGF was performed for a 320 s time window starting from the origin time. After frequency filtering (0.02-0.2 Hz) and rejection of stations with low-or high-frequency artifacts, results of 11 stations could be used to estimate their STF duration (9 NARS-Botswana stations plus two stations operated by the U.S. Geological Survey). Figure 7a shows the deconvolution results. The STF is interpreted as where L r is the (horizontal) rupture length, c r the rupture velocity and θ r the rupture propagation direction. Given that our deconvolutions are filtered between 0.02 and 0.2 Hz, an additional term W is incorporated to account for pulse broadening, resulting in shifts of the start and stop times. We estimate the pulse width W caused by the band pass filter to be 4.9 s (Figure 7a). With a lower crustal shear wave speed c of 3.9 km/s at the source depth of 29 km , the following parameters were fitted: L r = 12.3 km, c r = 2.27 km/s, θ r = 314° (wih an RMS error of 0.295 s). We note that the values of L r and c r are strongly dependent on c but that θ r is tightly constrained by the good azimuthal station coverage. The rupture propagation direction (to the northwest) agrees with the boundary between the Limpopo Belt and the Kaapvaal Craton, the distribution of the aftershocks (all northwest of the mainshock), as well as the strike of the nodal planes of the main event and the aftershocks.

Discussion
Events of the 2014-2018 NARS-Botswana catalog ( Figure 2) demonstrate that Botswana seismicity does not predominantly occur in the Okavango Rift Zone, but appears to be wider spread. In particular, the 3 April 2017, mainshock and its aftershocks indicate the presence of a seismically active zone around the southern boundary of the Limpopo Belt with the Kaapvaal Craton, a region that was unknown to be active before 2017.   The natural seismicity of our Botswana catalog is added in Figure 4b (events located in mining regions were discarded from our catalog). Although most of the events have magnitudes smaller than 4, their distribution seems to fill the region between the seismicity of the Okavango Rift Zone in the north and that of the Zoetfontein Fault in the northern part of the Kaapvaal Craton. This zone, with many events around the tectonic boundary between the Kaapvaal Craton and the Limpopo Belt, is mostly related to the 3 April 2017, mainshock and its aftershocks, but it also includes two events prior to the mainshock (in 2014 and 2015). The zone is not recognized in any previous maps of Botswana seismicity.
Double-difference relocation of the 3 April 2017, aftershocks resulted in a narrow distribution with a northwesterly orientation dipping toward the northeast (Figures 3 and 8c). There is no direct evidence for the presence of faults in this region because the Precambrian basement is overlain by a thick cover of Kalahari sands. The association of the mainshock and aftershocks to potential faults therefore relies on the interpretation of aeromagnetic and gravity data. The mainshock is located in the Southern Marginal Zone of the Limpopo Belt as inferred from gravity data (Kolawole et al., 2017;Ranganai et al., 2002). This zone is bounded by the Dinokwe Thrust in the southwest and the Mahalapye Shear Zone in the northeast (Figure 8a) and was formed during the Paleoproterozoic collision of the Kaapvaal and Zimbabwe cratons (Ranganai et al., 2002). Kolawole et al. (2017) analyzed and inverted aeromagnetic data and found that the hypocenter of the mainshock coincided with a faultlike feature in their model, which they named the Moiyabana Fault. Figure 8b presents a cross-section from their paper with our relocated events plotted on top. It shows that the mainshock and most of the aftershocks follow the Moiyabana Fault, confirming the interpretation by Kolawole et al. (2017). Nevertheless, it should be noted that the depth distribution by double-difference relocation is biased by the original event depths, most of which were fixed at 10 km which is too shallow for many of the lower crustal events (Table 1). The actual distribution is therefore likely to be deeper. Moreover, the depth resolution of aeromagnetic inversions is also limited. This means that the lateral locations of the features shown in Figure 8b are better constrained than their vertical positions. The figure further shows that not all the aftershocks follow the interpreted Moiyabana Fault, there is an additional faint branch of low magnitude events with an opposite dip that delineates a potential antithetic fault. This interpretation is strengthened by the aeromagnetic model because this branch approximately outlines the southwestern boundary of a body of high magnetic susceptibility. Thus, the relocated aftershock distribution not only confirms the presence of the Moiyabana Fault but also suggests a more complex fault system than previously inferred from gravity and aeromagnetic data.
We compared the depths of the main-and aftershocks obtained by double difference relocation to those that were obtained by waveform fitting (Table 1) and found them to be consistently 10-12 km shallower. As stated previously, such a shift is expected for events that occurred in the lower crust but were given a fixed depth of 10 km as input for the double difference inversion because the inversion keeps the centroid location of the earthquake cluster fixed. It suggests that the double difference depths of all events are roughly 11 km shallower compared to their more reliably estimated waveform depths.
The focal mechanisms of the waveform-fitted events are also consistent with the Moiyabana Fault. The nodal planes with strikes between 310° and 350° not only match the general strike of this fault, their dips between 32° and 38° toward the northeast also agree with the dip of the Moiyabana Fault in the lower part of the crust (Figure 8b).
The EGF analysis of the mainshock yielded a rupture propagation direction toward the northwest (θ r = 314°). This direction not only matches the strike of the Moiyabana Fault, it is also consistent with the distribution of aftershocks that mostly occurred northwest of the main event. The rupture propagation direction and the aftershock distribution both suggest that faulting progressed toward the northwest, not only for the main event, but also for the aftershocks.
The focal mechanisms of the events obtained by waveform modeling are all of normal faulting type with a NW-SE strike (Figure 6), indicative of NE-SW extension. This extension direction is consistent with the no-net rotation absolute plate motion direction (Argus et al., 2011;Gripp & Gordon, 1990;Wang et al., 2018). In addition, it is consistent with the NE-SW fast polarization directions obtained from shear wave splitting data (Yu et al., 2015) and surface wave data with asthenosperic sensitivity (Adam & Lebedev, 2012). It is therefore tempting to relate the focal mechanisms to a large scale NE-SW extensional crustal stress regime caused by mantle flow and coupling between the mantle and crust. However, such an interpretation is inconsistent with focal mechanism orientations on a larger scale, because they change from normal faulting with a NW-SE strike in central Botswana (this study) to normal faulting with a NE-SW strike in the ORZ (e.g., Bird et al., 2006). This 90° change in orientation over a distance of 300 km is not easily explained by a large-scale stress field in the crust. Instead, we observe that the strikes of the focal mechanisms, in the ORZ as well as those in central Botswana, match the orientations of their fault systems. It implies that local conditions and structural features largely determine the orientations of the focal mechanisms.
The question that then remains is: What caused the central Botswana earthquakes? Global studies on crustal earthquakes in intraplate settings often invoke fluid transport along existing faults Mooney et al., 2012), and several studies have shown that fluids or melt indeed play an important role explaining the seismicity along various segments of the EARS Lee et al., 2017;Lindenfeld et al., 2012;Oliva et al., 2019). Whereas these segments of the EARS are all seismically active and often volcanic in nature, the mainshock in central Botswana occurred in an area without any previously known seismicity. Moreover, it occurred in the lower crust at a depth of 29 km, which is considered to be too deep for brittle failure (W.-P. Chen & Molnar, 1983). Yet, lower crustal earthquakes are not uncommon in cratonic areas of the EARS (Craig et al., 2011;Yang & Chen, 2010). The occurrence of lower crustal earthquakes has been explained by fluid-induced eclogitization, that is transformational faulting caused by the breakdown of a dry metastable granulite facies rock to eclogite (Jackson et al., 2004). In recent years there is growing evidence that supports this mechanism based on observations on exposed rocks (Austrheim, 1987 (Incel et al., 2019), petrological and geophysical modeling (Henry et al., 1997;Hetényi et al., 2021;Malvoisin et al., 2020) and seismic studies (Michailos et al., 2021;Priestley et al., 2008). It can be a top-down process (Jamtveit et al., 2018(Jamtveit et al., , 2016Moecher & Steltenpohl, 2009), but in central Botswana it was a bottom-up process with lower crustal initiation of the mainshock at 29 km depth, and upward migration of the aftershocks in a northwesterly direction.
Eclogitization of the lower crust, as suggested here, requires fluid infiltration to initiate the mainshock by a chemical reaction, where a tiny amount of fluid is sufficient to start the process (Incel et al., 2019;Jackson et al., 2004). The transformation to eclogite will lead to a loss in strength and is accompanied by an increase in porosity and permeability caused by the volume reduction. This allows the fluids to infiltrate further, and the process is repeated, creating sequential aftershocks. Moorkamp et al. (2019) were unable to detect a mantle source of fluids in their seismic and resistivity models, but the more detailed tomographic study by Fadel et al. (2020) using NARS-Botswana data revealed a thin low shear-velocity anomaly in the uppermost mantle that links the main event to a large strong low velocity anomaly in the mantle beneath the eastern ORZ (see their Figure 4). This low velocity anomaly in the uppermost mantle may be interpreted as the mantle source of fluids. The low velocity anomaly beneath the ORZ was also found in several other (lower resolution) tomographic P and S wave velocity models (Ortiz et al., 2019;White-Gaynor et al., 2020), and may be related to decompression melting caused by lithospheric thinning in response to tensional stresses (Yu et al., 2017). The normal-faulting focal mechanisms are consistent with the tensional stress regime that is required for such a scenario.
The earthquakes in central Botswana suggest that the Kalahari Craton, formed by the amalgamation of the Kaapvaal Craton, Limpopo Belt and Zimbabwe Craton during Paleoproterozoic collision (2.7-2.6 Ga), has started to break up along one of its ancient suture zones. How this will develop in the future is unclear. It is generally believed that the African superplume plays an important role in rifting the African continent (Forte et al., 2010;Koptev et al., 2016). Nevertheless, geodynamic studies agree that lithospheric stresses arising from mantle coupling and lithospheric buoyancy forces are insufficient to initiate rifting (Ghosh & Holt, 2012;Stamps et al., 2018) and that pre-existing weak zones are needed for rifting to occur (Celli et al., 2020;Kendall & Lithgow-Bertelloni, 2016;Stamps et al., 2018). The earthquakes in central Botswana investigated in this study suggest that the southern margin of the Limpopo Belt is such a zone that is currently weakened by upward fluid transport from the mantle. The earthquake distribution of our catalog further suggests that the southwestern branch of the EARS may continue from the ORZ southwards along the southern boundary of the Limpopo Belt with the Kaapvaal Craton, in agreement with the hypothesis by Malservisi et al. (2013) based on geodetic data.

Conclusion
The deployment of 21 seismic stations of the NARS-Botswana network (2013-2018) allowed several seismic studies that provide new insights into the seismotectonics of the region. We created a new 4-year catalog of seismic events with an estimated magnitude of completeness of 3.1. The catalog confirms the seismic activity of the Okavango Rift Zone, which is currently interpreted as the terminus of the southwestern branch of the East African Rift System. In addition, the catalog lists the events in central Botswana that are related to the unexpected Mw 6.5 earthquake of 3 April 2017. Double-difference relocation of the aftershocks of this main event shows that a fault zone along the southern margin of the Limpopo Belt is reactivated. This fault system, formed as a thrust system during Paleoproterozoic collision of the Kaapvaal and Zimbabwe cratons (2.7-2.6 Ga), is now reactivated within a tensional stress regime as is evident from the focal mechanisms. The initiation of the mainshock at a lower crustal depth of 29 km with rupture propagation to the northwest along the Moiyabana Fault, together with the upward migration of the aftershocks along its fault system, points to eclogitization of a dry metastable granulite facies rock triggered by fluid infiltration. The presence of a low velocity anomaly in the uppermost mantle  provides additional support for this scenario. Whereas fluid transport has previously been suggested for earthquake occurrence in stable continental regions (Craig et al., 2011), this is, to our knowledge, the first study that interprets eclogitization in the lower crust from a main-and aftershock distribution. Furthermore, we speculate that rifting in the Okavango Rift Zone will continue southwards along the southern boundary of the Limpopo Belt, thus extending the southwestern branch of the East African Rift Zone into central Botswana. New data from the Botswana Seismic Network may help to resolve this hypothesis.

Data Availability Statement
The data used in this study are available from the IRIS Data Management Center: http://service.iris.edu/fdsnws/ dataselect/1/.