Crustal and Upper Mantle Shear Wave Velocity Structure of Botswana: The 3 April 2017 Central Botswana Earthquake Linked to the East African Rift System

Rayleigh wave group and phase velocity measurements obtained from ambient noise and earthquake data at 51 broadband stations were used to construct the first 3‐D crustal and upper mantle shear wave velocity model of Botswana. The model shows low crustal velocities associated with the Passarge and Nosop sedimentary basins, whereas the Kaapvaal, Zimbabwe, Maltahohe, and Congo Cratons are recognized by high mantle velocities. The lowest upper mantle shear wave velocity, beneath northeastern Botswana, is associated with the southwestern branch of the East African Rift System. This low‐velocity mantle anomaly appears to be linked to the crust of the Okavango Rift Zone and the location of the 3 April 2017 Mw 6.5 earthquake in central Botswana. We suggest that fluids or melt at the base of the crust from the southward continuation of the East African Rift Zone triggered the intraplate earthquake in an extensional tectonic setting.

. (a) Topography of southern Africa with the branches of the EARS (Chorowicz, 2005), Lake Kariba, the focal mechanism of the 3 April 2017 earthquake and outlines of Botswana and Figure 1b. (b) Tectonic map of the region with the distribution of seismic stations. (c) Tectonic units of the study area (Okavango Rift Zone with gray outline). (d) Sedimentary thickness map from aeromagnetic data after Pretorius (1984) and earthquake distribution. The 3 April 2017 earthquake and its aftershocks are shown in red.
The tectonics of Botswana today is influenced by the southward propagation of the East African Rift System (EARS). It has long been argued that the EARS terminates at the incipient Okavango Rift Zone in Botswana ( Figure 1) (Leseane et al., 2015;Reeves, 1972;Yu, Liu, et al., 2015 and references therein). It is therefore remarkable that the second largest earthquake in Botswana's history, the magnitude 6.5 earthquake of 3 April 2017 (Figure 1), was located close to the border of the Kaapvaal Craton (Kolawole et al., 2017), approximately 400 km away from the Okavango Rift Zone. The event occurred at approximately 29 km depth and had an extensional mechanism (Gardonio et al., 2018;Materna et al., 2019;Midzi et al., 2018). The mechanism of the earthquake suggests a relation with a rifting system. This raised the question whether the earthquake might be linked to the EARS. To address this, detailed imaging of the structure beneath Botswana is required. In this study, we used data from 51 stations in Botswana and Namibia, to obtain the first 3-D shear wave velocity model for the entire country of Botswana. We used interstation Rayleigh wave dispersion measurements from ambient noise and teleseismic surface wave data to construct our model.

Seismological Data
We used data recorded by 51 broadband stations from four seismic networks ( Figure 1b; Table S1 of the supporting information). The largest part of the data is obtained from the NARS-Botswana network (NARS, 2019 https://doi.org/10.7914/SN/NR), deployed from November 2013 to March 2018. It consists of 21 seismic stations that are distributed over the different tectonic units in Botswana. In addition, we used 17 stations from the Seismic Arrays for African Rift Initiation (Gao et al., 2013, https://doi.org/10.7914/SN/XK_2012), a network that was deployed between May 2012 and June 2014 to image the crustal and mantle structure underneath the Okavango Rift Zone. Furthermore, to enhance ray coverage to the west, we used 10 stations from the AfricaArray Namibia network that has been recording between May 2015 and January 2019 Finally, we used data from two stations from the Global Seismological Network (Butler et al., 2004, https

Methodology
We used two complementary methods to derive the 3-D shear wave velocity structure underneath Botswana: (1) ambient noise tomography (Bensen et al., 2007) where interstation Rayleigh wave group-and phase velocity measurements are obtained from ambient noise data at relatively short periods and (2) Helmholtz tomography (Jin & Gaherty, 2015) that yields Rayleigh wave phase velocity maps at longer periods from teleseismic earthquakes. The combined group and phase velocity maps provide a data set that is used to invert for the crustal and upper mantle shear wave velocity structure. The two methods to construct the group and phase velocity maps and the inversion scheme are briefly described in the following subsections.

Ambient Noise Tomography
Ambient noise tomography has become a well-established imaging technique since the pioneering work of Shapiro and Campillo (2004) and Shapiro et al. (2005). Here, we adopted the methodology of Poli et al. (2013) to extract the phase and group velocity measurements for all possible 1,225 interstation paths. The data were corrected for instrument response, filtered from 2.5 to 50 s period, decimated to 1 s, and subjected to spectral whitening. For each station pair, the cross-correlation function was calculated for 4 hr segments, excluding segments with amplitudes larger than 10 times the standard deviation of the windowed data and stacked over 3 month periods to minimize the effects of seasonal variations of the ambient noise sources (Bensen et al., 2007). Then, an automated frequency-time analysis was applied to cross-correlation functions with a signal-to-noise ratio of 10 or larger (Bensen et al., 2007) to extract the phase (3-35 s) and group (3-30 s) velocity dispersion curves. The dispersion curves for each interstation path were then averaged to obtain the final curves with their standard deviation. Examples of measured dispersion curves are shown in Figure S1 of the supporting information.
We inverted all interstation dispersion data for group and phase velocity maps on a 1 • by 1 • grid using the method of Barmin et al. (2001). A penalty function is minimized depending on data misfit, model smoothness, and perturbation to the reference model. The spatial smoothing is governed by weight and Gaussian smoothness width . The perturbation to the reference model, or norm damping, is controlled by weight and the path density contribution which is governed by parameter . We varied the inversion parameters systematically considering the data misfit and model resolution to find optimum values, seeking a balance between the smoothness of the maps and data misfit and finally selected = 100, = 100 km which is approximately equal to the cell size, = 1.0 and = 0.5. Figure S1 of the supporting information shows the phase velocity inversion results for a period of 22 s as an example. Checkerboard tests were carried out that confirmed that the maximum resolution of the ambient noise group and phase velocity maps is 1 • ( Figure S2 of the supporting information). Figures 2a-2h show the group and phase velocity maps for periods of 5, 10, 20, and 30 s.

Helmholtz Tomography
We used the Automated Surface Wave Phase Velocity Measuring System that was developed by Jin and Gaherty (2015) to estimate phase velocities for periods between 30 and 120 s.
Automated Surface Wave Phase Velocity Measuring System is applied to the vertical components of teleseismic earthquake records at stations with interstation distances between 50 and 600 km to calculate the multichannel cross correlations of the fundamental mode Rayleigh wave. We used records between June 2012 to May 2017 for earthquakes with M w > 6, at 30-160 • distance from the center of the study area with depths shallower than 50 km. A window function is applied to the cross-correlation waveforms to determine the coherency of the fundamental mode Rayleigh wave signal. The signals with a coherence larger than 0.5 are subjected to narrow band filters around target frequencies using low and high cutoffs of ±10% of the central 19 target periods between 30 and 120 s for which phase delays and amplitude variations are measured.
Initial phase delays for each period are estimated by minimizing the misfit between the coherent band-filtered waveforms and a Gaussian-shaped wavelet that includes phase and group delays. Phase delays with a difference larger than 10 s than predicted for the average phase velocity at each period are removed from the data set to avoid cycle skipping and eliminate poor quality data. Autocorrelations are then used to correct for biases in the initial phase delays that are caused by the window function. Data from 323 earthquakes (with nonuniform azimuthal coverage) were finally used in the tomographic inversion ( Figure S3 of the supporting information). The phase delays are inverted to obtain initial phase velocity maps using the Eikonal equation with a smoothing parameter that gives a smoothing length equal to 25% of the modeled wavelength at each period. Amplitude information is used in the inversion using the Helmholtz equation. The (frequency-dependent) amplitudes of the (frequency-dependent) wavelets provide a good approximation of the amplitude of the power spectrum. These amplitude measurements are used to correct the phase velocity maps for the effects of focusing and defocusing using the Helmholtz equation. For more details the reader is referred to Jin and Gaherty (2015). The final phase velocity maps for periods of 40 to 120 s are shown in Figures 2i-2p.

Combining the Ambient Noise and Helmholtz Phase Velocity Maps
At this stage, we had 3-35 s phase velocity maps from ambient noise and 30-120 s from Helmholtz tomography. The phase velocity dispersion maps in the overlapping range of 30-35 s were merged to provide a smooth transition from the ambient noise to the Helmholtz tomography maps. For the overlapping period range, relative weights were used (0.7 to 0.3 for ambient noise and 0.3 to 0.7 for Helmholtz tomography for periods of 30 to 35 s) to facilitate a smooth transition between the ambient noise and the Helmholtz tomography phase velocity data. This scheme was preferred over other merging schemes because it prevented creating abrupt changes in the dispersion curves that might lead to artifacts at the later inversion stage. Figure S4 of the supporting information shows an example of the merging scheme for the 35 s phase velocity maps.

Joint Inversion of Group and Phase Velocity Dispersion Curves
The group and phase velocity maps provide the local variations of the group and phase velocity at each grid point. Group and phase velocity curves at each grid point were jointly inverted to produce the local 1-D shear wave velocity structure at that location. The inversion of all grid points resulted in a 3-D model of the crust and upper mantle. The inversion was carried out using iterative damped least squares optimization based on the Levenberg-Marquardt algorithm (Moré, 1978) similar to Greve et al. (2014). The 1-D inversion approach aimed to minimize an objective function depending on measured data d and 1-D subsurface model m that consists of three misfit terms: the data misfit, the vertical smoothness, and the norm of the perturbation from the reference model, m ref : (1) where the first term represents data misfit. G(m) is the calculated dispersion data for model m, and d is the measured dispersion data. The model smoothness term is calculated using the first difference of the 1-D model following Julià et al. (2000), and is a parameter controlling the contribution of the model smoothness term to the objective function. The model perturbation term is calculated using the difference between the estimated model and the 1-D reference model, and is a parameter controlling the contribution to the norm of the perturbation from the reference model to the objective function.
We inverted the group (3-30 s) and phase (3-120 s) dispersion curves to estimate the 3-D subsurface structure of Botswana in two steps. We first inverted the average group and phase velocity curves to estimate the average 1-D shear wave velocity structure in Botswana. We used a linear interpolated version of AK135 (Kennett et al., 1995) as a starting model with a transition between the crust and mantle at 40 km depth based on the average crustal thickness estimates over Botswana (Fadel et al., 2018; Figure S5 of the supporting information), similar to Kgaswane et al. (2009). We did not impose a hard discontinuity for the Moho as suggested by Ma and Clayton (2014) since surface waves lack strong sensitivity to discontinuities. The inversion was carried out using = 0.5 based on several systematic tests to keep the balance between the model smoothness and dispersion curves fitting. We used = 0 to keep the inversion of the average 1-D model as independent as possible from the starting model (Julià et al., 2000;Ma & Clayton, 2014). The inversion was carried out for shear wave velocity only. For the crustal part with V s < 4.1 km/s, the P wave velocity and density were coupled to shear wave velocity using the empirical relations found by Brocher (2005). For the deeper part of the model with V s > 4.1 km/s, the P wave velocity and density were kept in fixed ratios to shear wave velocity as the AK135 starting model. In the second step, the average 1-D subsurface model was used as a starting and reference model to invert the group and phase velocity curves at all grid points in Botswana. We used = 0.5 and = 0.5. The choice for was made to prevent large deviations from the average 1-D model of Botswana.
The 1-D models were discretized into 22 layers from the surface until the mantle discontinuity at 410 km depth: 2.5 km steps for the top 15 km of the crust, 5 km from 15 to 50 km depth, 10 km from 50 to 100 km, and 20 km layers for the rest of the model. Whereas the 1-D inversions were carried out to 410 km depth to prevent leakage of deeper structure artifacts to shallower depths, the final 3-D shear wave velocity model in Figure 3 is presented to 200 km depth, the region with the largest data sensitivity.
The resolution of the model was checked by a synthetic test. Using our 3-D shear wave velocity model as input, we calculated interstation group and phase velocity curves for the original data distribution. These synthetic data were used to produce group and phase velocity maps, which were subsequently inverted to obtain the output model. The results of this test for the four cross sections of Figure 4 are presented in Figure S6 of the supporting information. They demonstrate the robustness of the retrieved structures and show that the inversion mainly reduces the amplitude of the anomalies. This implies that the retrieved anomalies of our 3-D model are likely underestimated.

Results and Discussion
The final 3-D shear wave velocity model of Botswana (Figure 3) has interesting features some of which are better recognized in the cross sections of Figure 4.

The Crust
The first distinctive features in the crust are the low-velocity anomalies of the deep sedimentary Nosop (NB) and Passarge (PB) Basins. In the southwest, the Nosop Basin has a low shear wave velocity to ∼15 km depth (Figure 4a and A-A ′ ). The Passarge Basin in central Botswana has a low velocity that extends to a depth of ∼10 km (Figure 4a, B-B ′ and C-C ′ ). To our knowledge, the latest attempt to map these sedimentary basins was done by Pretorius (1984) using aeromagnetic data by Reeves and Hutchins (1982). The location and depth extent of our model agree with those of Pretorius (1984) (Figure 1d) and the global sediment map of Laske (1997).
Another distinctive feature in the upper and middle crust (0-30 km) is the crustal low-velocity anomaly of the Okavango Rift Zone (Figures 4a and 4b, C-C ′ and D-D ′ ). The low shear wave velocity is compatible with the high Vp/Vs ratio found by recent receiver function studies Fadel et al., 2018) and the shallow Curie depth of 8-15 km from 3-D inversion of aeromagnetic data (Leseane et al., 2015). Our results and the previous observations are consistent with the hypothesis of ascending hot fluids in the crust from the mantle through weak zones in the lower crust (black lines in Figure 4, C-C ′ and D-D ′ ) (see Leseane et al., 2015;Fadel et al., 2018;. It has been suggested that these ascending fluids trigger the earthquakes in the Okavango Rift Zone (Leseane et al., 2015, and references therein).

The Upper Mantle
The upper mantle structure of Botswana is dominated by high shear wave velocities that are related to the cratons and low velocities that extend from northern Botswana toward the south-southwest, until the border of the Kalahari Craton (Figures 3 and 4). The Kalahari Craton is characterized by high shear wave velocities in the east and southeast. It consists of two cratonic blocks, the Kaapvaal Craton (KC) and the Zimbabwe Craton (ZC), that are sutured by the Limpopo Belt (LB) in between (Figures 4a-4d). These cratons have deep roots that extend to (at least) 200 km depth (Figures 4c and 4d, A-A ′ , B-B ′ , and C-C ′ ). For Kaapvaal, we can compare our results with surface wave studies using data from the Southern Africa Seismic Experiment with a limited number of stations in southeastern Botswana (Carlson et al., 1996). The depth-dependent structure of our model is in general agreement with previous studies (Adams & Nyblade, 2011;Li & Burke, 2006;Ravenna et al., 2018) with mantle shear wave velocities that increase to 4.7-4.8 km/s at ∼150 km before decreasing again. The gradual velocity increase to 150 km may be explained by an increase in garnet content as argued by Ravenna et al. (2018). Furthermore, our results suggest that the deep root of the Zimbabwe Craton (at ∼200 km depth) has a higher shear wave velocity than the Kaapvaal Craton (Figure 4d). This is in agreement with the study of Li and Burke (2006) and may reflect a difference in composition or temperature between the two cratons.
The northwestern part of Botswana shows a high-velocity anomaly at 60-120 km depth and another, potentially related, high-velocity anomaly at 170-200 km depth (CC in Figure 4, C-C ′ and D-D ′ ). These anomalies may be associated to the Congo Craton as suggested from the high P wave anomaly in the study by Yu et al. (2017). The anomalies are separated by the low-velocity anomaly from the EARS, entering Botswana from the northeast and extending to western and central Botswana at depths of 140-180 km (Figures 4i-4k). This low-velocity anomaly may have affected the edge of the craton.
The western part of section A-A ′ of Figure 4 shows a positive shear wave velocity anomaly, approximately between 50 and 200 km depth, that seems to extend toward the center of the section. Some studies have suggested the presence of the Maltahohe microcraton below the Rehoboth Province (Begg et al., 2009;Wright & Hall, 1990). This seems confirmed by the high P and S wave velocity anomalies of the study by Ortiz et al. (2019) that was attributed to a thick cratonic lithosphere. A major part of the Rehoboth Province was formed during the Proterozoic (2.2-1.9 Ga) around an Archean nucleus (van Schijndel et al., 2011). Furthermore, the Rehoboth Province has a similar lithospheric conductivity as the western block of the Kaapvaal Craton, as shown by a magnetotelluric study of Muller et al. (2009). The high-velocity anomaly of this study, together with these pieces of evidence, supports the interpretation of the presence of the Maltahohe microcraton underneath the Rehoboth province.
Between the Kaapvaal Craton in the east and the deep cratonic root of the Rehoboth Province in the west, there is a zone of relatively low shear wave velocities (KB in Figure 4 A-A ′ ). It was also imaged as a small but negative velocity anomaly in the P and S wave tomography study by Youssof et al. (2015) and is probably associated to the Kheis Belt and the Okwa Block. Petrological analysis of rock samples from the Okwa Basement Complex by Mapeo et al. (2006) has shown that the Okwa block was emplaced during a narrow time span at 2, 056 ± 3 Ma. This timing corresponds with the intrusion of the Bushveld Complex in the Kaapvaal Craton in South Africa at ∼2.05 Ga, suggesting a link between the two events.
The strongest negative mantle anomaly of our model is found in the northeast of Botswana, beneath the northeastern tip of the Okavango Rift Zone (ER in Figures 4c and 4d, B-B ′ and D-D ′ ). The Okavango Rift Zone is an incipient rift and is often interpreted as the southern terminus of the EARS (e.g., Bufford et al., 2012;Ortiz et al., 2019;Kinabo et al., 2008;Leseane et al., 2015;Modisi et al., 2000;Reeves, 1972;Scholz et al., 1976;. Our model is in agreement with this interpretation because the location of the low-velocity anomaly matches the continuation of the southwestern branch of the EARS with Lake Kariba as its last surface expression (Figure 1a). The low-velocity anomaly in the mantle seems to be connected to more shallow low-velocity anomalies (∼30-50 km) in the lower crust and uppermost mantle of the Okavango Rift Zone (black solid lines in Figure 4, C-C ′ and D-D ′ ). This can be interpreted as support for ascending fluids that have been proposed for the Okavango Rift (Fadel et al., 2018;Leseane et al., 2015;. Toward the south and southwest, at 120-200 km depth, the low-velocity anomaly extends to the border of the Kaapvaal and Zimbabwe Cratons (Figures 3h-3l). This suggests that the temperature or compositional changes related to the EARS anomaly end at these cratons at these depths.

The 3 April 2017 M w 6.5 Earthquake
Whereas the EARS low-velocity anomaly does not extend beyond the Kaapvaal and Zimbabwe Cratons at depths larger than ∼120 km, the pattern is less clear at shallower depths (Figures 3d-3g). In fact, Figure 4 B-B ′ shows a low-velocity zone extending from the EARS anomaly to the location of the 3 April 2017 M w 6.5 earthquake (dashed line) which suggests a link between the two. Our resolution test shows that this anomaly is resolved (Figure S6 B-B ′ of the supporting information).
The earthquake was a normal faulting event in the lower crust, close and parallel to the border of the Kaapvaal craton (Materna et al., 2019;Midzi et al., 2018;Moorkamp et al., 2019;Kolawole et al., 2017). Here, we suggest that it is associated to the EARS, where the low-velocity zone acts as a conduit of fluids or melt from the low-velocity mantle, activating brittle failure in the crust. A low-velocity zone of 3% at 50 km depth would require a 150 • C temperature increase, a 10% increase in Fe content, or a 0.3-4% increase in melt (depending on the interconnectivity and geometry of the melt inclusions) (Goes et al., 2000). Such temperature or compositional variations are difficult to explain within a narrow zone of 20-30 km, making partial melt, or another fluid phase, most likely. The scenario of fluid activation was recently proposed by Gardonio et al. (2018), but our model shows the link to the EARS.

Conclusion
We presented the first 3-D shear wave velocity model of Botswana using Rayleigh wave dispersion data obtained from ambient noise and teleseismic earthquakes largely recorded by stations of temporary networks in the region. The group (3-30 s) and phase velocity (3-120 s) data were jointly inverted to derive a 3-D shear wave velocity model to 200 km depth. The model shows upper crustal low-velocity anomalies that are associated to the Passarge Basin in central Botswana and the Nosop Basin in the southwest. Low velocities in the upper crust beneath the Okavango Rift Zone may be caused by fluids in the crust as suggested by Leseane et al. (2015). The Kaapvaal and Zimbabwe Cratons are characterized by high velocities extending to at least 200 km depth. In addition, we find a positive mantle velocity anomaly beneath the Rehoboth Province, which is interpreted as the enigmatic Maltahohe Craton, and there are high mantle velocities related to the Congo Craton. Lastly, our model shows a low-velocity anomaly in the mantle that, through its location in northeast Botswana, is related to the EARS. This anomaly seems to be connected to the crust of the Okavango Rift Zone through zones of low velocity. These zones may represent the conduits of fluids which reduce the seismic velocity in the crust and trigger the local seismicity. Another low-velocity zone links the EARS mantle anomaly to the location of the 3 April 2017 M w 6.5 earthquake close to the border of the Kaapvaal Craton in central Botswana. We suggest that fluids associated with this southward continuation of the EARS anomaly triggered this intraplate lower crustal event in an extensional tectonic setting.