Numerous investigations of the mature segments of the East African rift system (EARS) have significantly improved our understanding of the structure and processes associated with well-developed continental rifts. In contrast, knowledge of rifting processes at their early stage is still significantly limited. Here we present results from a teleseismic P-wave tomography investigation of the incipient Okavango rift zone (ORZ), which is located at the southwestern terminus of the EARS. P-wave relative travel-time residuals recorded by 17 recently deployed portable seismic stations were manually picked and inverted for three-dimensional upper-mantle and mantle transition-zone tomographic images beneath the ORZ and its adjacent areas. High-velocity anomalies probably representing cratonic lithosphere are visible under the Congo and Kalahari cratons, extending to depths of ∼250–350 km. The tectonic boundary of the Congo craton is observed along the western edge of the ORZ. A localized low-velocity anomaly of about –1% in magnitude is revealed in the upper asthenosphere beneath the ORZ, which is interpreted to represent decompression melting induced by lithospheric thinning. The results support the notion that the initiation and early-stage development of the ORZ are mostly due to lithospheric stretching resulted from the relative motion between the Archean Congo and Kalahari cratons along preexisting ancient orogenic zones.

Continental rifting plays a significant role in the evolution of the Earth (e.g., Ziegler and Cloetingh, 2004). Our current comprehension of rifting processes is primarily obtained from studying relatively mature rifts such as the Main Ethiopian and Kenyan rifts of the East African rift system (EARS). In contrast, the mechanisms and processes associated with the initiation and early-stage development of rifting are inadequately investigated, preventing a thorough understanding of the whole rifting cycle.

Most mature rifts are characterized by a low-velocity (low-V), potentially partially molten upper mantle, which has been interpreted as the result of active upwelling of mantle plumes or decompression melting due to passive lithospheric stretching (White and McKenzie, 1989). In addition, the time of the emplacement of the sub-rift low-V zone as well as the amplitude of the anomaly beneath emerging young rifts are still poorly constrained. Here we seek to fill the knowledge gap of the evolution of rifting at its earliest stage by conducting a seismic body-wave tomographic investigation of the incipient Okavango rift zone (ORZ) in southern Africa.

The ORZ is believed to be one of the youngest continental rifts in the world (Reeves, 1972; Reeves and Hutchins, 1975; Scholz et al., 1976). The initiation of the ORZ is estimated to have been in the Holocene, approximately between 120 ka and 40 ka (Modisi et al., 2000; Moore and Larkin, 2001). The ORZ is considered to be the terminal section of the southwestern branch of the EARS (Reeves, 1972; Modisi et al., 2000), having primarily developed within the Paleoproterozoic Magondi belt, the Paleo-Mesoproterozoic Rehoboth province, and the Neoproterozoic Damara belt between the Archean Congo craton to the northwest and the Kalahari craton (which includes the Kaapvaal and Zimbabwe cratons and the Limpopo belt in the study area) to the southeast (Fig. 1; Begg et al., 2009).

Based on analyses of Shuttle Radar Topography Mission and aeromagnetic data, Kinabo et al. (2008) proposed that the localization and development of rift-associated faults beneath the ORZ were mainly controlled by preexisting basement fabrics. Similarly, Khoza et al. (2013) suggested that the ORZ was generated in preexisting orogenic zones, on the basis of the lack of elevated mantle conductivity from two-dimensional (2-D) and three-dimensional (3-D) inversions of magnetotelluric data. Gravity and magnetic investigations concluded that hydrothermal metasomatism and strain-weakening of the lithospheric mantle promoted the initiation of the ORZ (Leseane et al., 2015). Recently, using broadband seismic data from Seismic Arrays for African Rift Initiation (SAFARI; Gao et al., 2013), Yu et al. (2015a) showed that the Moho shallows by 4–5 km and that there are no thermal anomalies in the mantle transition zone beneath the ORZ, ruling out the existence of one or more mantle plumes originating from the lower mantle directly beneath the study area (Yu et al., 2015b). The observations, however, cannot readily rule out the influence on early-stage rifting by one or more mantle plumes entering the upper mantle from outside the study area, e.g., farther north along the EARS (Ebinger and Sleep, 1998). In addition, mantle anisotropy beneath the ORZ and its adjacent areas is characterized by dominantly northeast-southwest fast orientations, which are consistent with those observed in the surrounding areas in southern Africa, suggesting a lack of small-scale mantle flow circulation associated with the ORZ (Yu et al., 2015c). The evolution of the ORZ has thus been proposed to follow a passive rifting mode and to be associated with relative movements of the Archean cratons (Yu et al., 2015c).

Due to the existence of a thick Kalahari sedimentary cover (Cooke and Verstappen, 1984; Ringrose et al., 2005), results from previous geological studies are controversial regarding the location of the southern boundary of the Congo craton (Singletary et al., 2003; Begg et al., 2009; McCourt et al., 2013). Investigations by Begg et al. (2009) show that the boundary of the Congo craton is situated outside of Botswana, while those by Singletary et al. (2003) and McCourt et al. (2013) indicate that the Congo craton occupies part of northwestern Botswana (Fig. 1). In contrast, magnetotelluric studies show that the lithospheric thickness beneath the northern and southern ORZ is ∼160 km (Muller et al., 2009) and 180 km (Miensopust et al., 2011), respectively, and the boundary of the Congo craton approaches the northwestern border of the ORZ (Khoza et al., 2013).

One of the main seismological observables that is still lacking is the upper-mantle velocity structure beneath the ORZ. In this study, we conduct a teleseismic P-wave tomography investigation of the upper-mantle and mantle transition zone velocity structure in the vicinity of the ORZ. The tomographic results shed new light on rifting initiation and early-stage development, as well as on the location of cratonic boundaries in the lithosphere.

Data used in the study were recorded by 17 broadband seismic stations belonging to the SAFARI experiment (Gao et al., 2013). The stations were deployed along a ∼760-km-long northwest-southeast profile traversing the ORZ (Fig. 1) with a sampling frequency of 50 Hz and a continuous recording time of 2 yr starting from May 2012. The epicentral distance of the events used in the study ranges from 25° to 95°, and the minimum body-wave magnitude is 4.3, as determined by the U.S. Geological Survey ( The vertical-component seismograms with a high signal-to-noise ratio are windowed based on the theoretical arrival times computed using the IASP91 Earth model (Kennett and Engdahl, 1991), and the P-wave arrival time is manually picked. The picking accuracy is estimated to be 0.1–0.2 s, as demonstrated in the example shown in Supplemental Figure S11. The picking procedure includes several interactive steps. For a given event, we start with a careful visual check of all the vertical-component seismograms to reject all lacking a clear P-wave onset. Those remaining are then zoomed (Supplemental Fig. S1C [see footnote 1]) and the onset is picked. The seismograms are then aligned according to the picked arrival time to identify mis-picked traces, and necessary corrections are made until the alignment is satisfactory. Each of the events used in the study was recorded by at least five stations. A total of 1736 direct P-wave arrivals from 148 teleseismic events (Fig. 2) are used for the tomographic inversion. The ray paths crisscross reasonably well in both the horizontal and vertical planes down to a depth of 700 km (Fig. 3).

In order to minimize the effects of the uncertainties in hypocenter locations and origin times, as well as velocity heterogeneities outside the study volume, we compute relative travel-time residuals for each event in the Supplemental Data2 by subtracting the mean residual from the raw residuals prior to the tomographic inversion. The tomographic method of Zhao et al. (1994) is then applied to image the 3-D upper mantle velocity structure under the ORZ and its surrounding areas. The IASP91 Earth model is adopted to be the starting one-dimensional velocity model for the tomographic inversion. The effect of surface topography and contributions from crustal thickness variations on the observed residuals are removed by interpolating a crustal model determined by a receiver-function study using data from the same stations (Yu et al., 2015a). An average crustal P-wave velocity of 6.5 km/s is adopted based on a previous seismic refraction study in the Kaapvaal craton (Durrheim and Green, 1992), and the length of the ray path in the crust is calculated by considering the angle of incidence calculated using the ray parameter. The small crustal thickness variation of ∼5 km across the profile leads to a maximum difference of ∼0.1 s in crustal contribution to the observed relative residuals. Most of the corrected relative residuals fall in the range of –0.4 to 0.4 s (Supplemental Fig. S2 [see footnote 1]), and relative to the Kaapvaal craton to the southeast of the rift, the ORZ exhibits a travel-time delay of ∼0.4 s (see Yu et al., 2015b, their figure 7), which is significantly smaller than the amount of delay associated with mature rifts such as Baikal (southeastern Russia, ∼1.3 s; Gao et al., 2003; Zhao et al., 2006), Rio Grande (southwestern USA and northern Mexico, ∼1.5 s; Gao et al., 2004), and the Tanzanian and Ethiopian sections of the EARS (>2.0 s; Ritsema et al., 1998; Bastow et al., 2005).

The mean travel-time residual for each of the events, which is removed from the raw travel-time residuals to obtain relative residuals, ranges from –0.8 to 4.5 s, with an average of 1.5 ± 0.1 s. Given the decent azimuthal coverage of the events (Fig. 2) and the associated small likelihood that the delay is dominantly caused by low-V anomalies near the hypocentral areas, this observation suggests an overall slower-than-normal mantle beneath southern Africa. Interestingly, previous surface-wave studies (e.g., Fishwick, 2010; Priestley et al., 2008) reveal normal to higher-than-normal velocities in the upper mantle beneath the study area. These observations are consistent with the existence of a large-scale low-V anomaly (the so-called African superplume) in the lower mantle beneath southern Africa (Ritsema et al., 1998).

A 3-D grid is set up in the study volume with a lateral grid interval of 1° and a vertical grid interval of 50–100 km. P-wave velocity (Vp) perturbations at the grid nodes are taken as unknown parameters, and the Vp perturbation at any point in the model is obtained by linearly interpolating the Vp perturbations at the eight grid nodes surrounding that point (Zhao et al., 1992). Theoretical travel times and ray paths are calculated using a 3-D ray-tracing technique (Zhao et al., 1992; Zhao, 2015), which is a combination of the pseudo-bending algorithm and Snell’s law. A conjugate-gradient algorithm (Paige and Saunders, 1982) with damping and smoothing regularization is employed to solve the large and sparse system of observation equations (Zhao et al., 1994; Zhao, 2015). The optimal values of the damping and smoothing parameters are determined based on the tradeoff curve between data variance reduction and model norm (Fig. 4; Hansen, 1992). The root-mean-square relative travel-time residuals before and after the inversion with 77 iterations are 0.32 s and 0.25 s, respectively (Supplemental Fig. S2B [see footnote 1]).

In order to evaluate the resolution scale and reliability of our tomographic model, we conducted checkerboard resolution tests (CRTs), which are a convenient way to examine the spatial resolution of tomographic images. The 3-D grid nodes are assigned with alternatively negative and positive Vp anomalies of 3%. The numbers of stations, events, and ray paths in the synthetic data set are the same as those in the real data set. Synthetic travel-time residuals are calculated for the input checkerboard model. Random noise with a standard deviation of 0.1 s is added to the synthetic data to simulate the picking errors in the observed data. The same tomographic algorithm is then applied for inverting the synthetic data to recover the input checkerboard model. Areas where the checkerboard model is well reconstructed are considered to have a high resolution.

Several CRTs with different grid intervals are conducted (Fig. 5; Supplemental Figs. S3–S8 [see footnote 1]). Considering the balance between the resolution scale and model reliability, a grid interval of 1° is chosen to conduct the optimal tomographic inversion of the real data. The CRT results show that the checkerboard pattern is well recovered for the depth range of 50–700 km in and around the ORZ (Fig. 5; Supplemental Figs. S3 and S4 [see footnote 1]). Due to the nearly vertically incident rays near the surface, the input anomalies in the shallower lithosphere (≤100 km depth) are generally resolved beneath the areas in and around the seismic array. With increasing depth, the crisscrossing of seismic rays gradually improves (Fig. 3) and the area with an acceptable resolution becomes larger (Fig. 5). However, the spatial resolution becomes significantly lower with depth in and around the ORZ area, where a limited ray-path coverage exists at greater depths (Fig. 3). This is the consequence of the fact that the linear array, which was designed to be perpendicular to the ORZ, is not optimally oriented toward the world’s main seismic zones. The resulting limited azimuthal coverage of the incoming rays could lead to biased velocity anomalies for nodes with strong seismic anisotropy. However, the study area is dominated by northeast-southwest–oriented azimuthal anisotropy (Yu et al., 2015c) which is parallel to the back-azimuth of the majority of the events used in this tomography study (Fig. 2). Consequently, the influence of azimuthal anisotropy may not be significant, although a quantitative evaluation of the degree of influence is impossible because the depth distribution of the azimuthal anisotropy is not known. As a whole, the CRT results (Fig. 5) suggest that the obtained tomographic images (Figs. 6 and 7) in the top 700 km depth are reliable.

Both low- and high-velocity bodies are clearly revealed in the upper mantle, and there are few or no significant low-V anomalies within the mantle transition zone, an observation that is consistent with the receiver-function results (Yu et al., 2015b). The velocity anomalies are dominantly in the range of –2% to 2% (Fig. 6). The distribution of the observed velocity anomalies is closely associated with the tectonic terranes. For the Congo craton which is in the northwestern part of the study area, low-V anomalies are constrained in the shallow lithosphere (<100 km depth), while the underlying depth range of 100–300 km is dominated by a high-V body (Fig. 7B). The most prominent feature observed on the cross-section (Fig. 7B) is a low-V zone in the depth range of 150–350 km beneath the ORZ. The anomaly covers the entire ORZ on the depth slices at 200 km and 300 km (Fig. 6) and possibly extends to the north of the ORZ. High-V anomalies, which are characteristic of the lithosphere, exist from the surface to the depth of ∼150 km beneath the ORZ area and to ∼250 km beneath the Kalahari craton (Fig. 7).

To confirm the main features of the tomographic images, we also conduct a restoring resolution test (RRT; Zhao and Hasegawa, 1993; Zhao et al., 1994). The RRT is similar to the CRT in principle, but instead of using the alternating velocity anomalies used in the CRT, the obtained tomographic results from the real data set are taken as the RRT input model. The resulting RRT image shows that the main features such as the low-V anomalies beneath the ORZ and the high-V bodies beneath the Congo and Kalahari cratons are well recovered (Figs. 7C; Supplemental Fig. S9 [see footnote 1]). More synthetic tests with different input models beneath the ORZ are conducted (Supplemental Figs. S10–S12 [see footnote 1]), and the test results suggest that the input models are well recovered in the upper mantle, whereas the resolution becomes lower in and below the mantle transition zone.

Lithospheric Thickness

The prominent high-V anomalies clearly show the depth extent of the cratonic roots and the lithosphere beneath the rift zone (Figs. 6 and 7). Under the commonly used assumption that high-V anomalies represent the cratonic lithosphere (James et al., 2001), our results suggest that the lithosphere beneath the Archean Congo craton may extend to an estimated maximum depth of 300–350 km, which is deeper than what was obtained from the 2-D and 3-D inversions of magnetotelluric data (Fig. 7B; Khoza et al., 2013). The inferred lithospheric thickness of the Kalahari craton is ∼250 km (Fig. 6), which is similar to the tomographic results in southern Africa, where the high-V mantle roots extend to depths of at least 250 km and locally 300 km (e.g., James et al., 2001). It is also consistent with results of a magnetotelluric experiment on the Zimbabwe craton, which shows a lithospheric thickness of 220 km (Miensopust et al., 2011). Our tomographic results suggest that the Congo craton terminates near the northern edge of the ORZ (Figs. 6C and 6D). The boundary of the Congo craton described by McCourt et al. (2013) and Singletary et al. (2003) is >100 km from what we determined from our tomographic images (Figs. 6C and 6D).

Our results show a significant change in the lithospheric thickness between the Congo craton and the ORZ where the high-V anomalies terminate at a depth of ∼150 km (Fig. 7), an observation that is consistent with the magnetotelluric results (Khoza et al., 2013). Magnetotelluric studies have also been conducted in the northern and southern parts of the ORZ within the Damara belt, resulting in the observation of a resistive 180-km-thick lithospheric mantle (Miensopust et al., 2011) and a lithospheric thickness of 160 km (Muller et al., 2009), respectively (Fig. 7B).

The depth of the upper-mantle high-V anomalies increases toward the northeast away from the ORZ (Figs. 6C–6E) and is consistent with the spatial variation of the apparent depths of the mantle transition zone discontinuities (Yu et al., 2015b). A recent shear wave splitting investigation indicates that the optimal depth of the center of the anisotropic layer beneath the central part of the seismic profile is approximately between 240 and 280 km (Yu et al., 2015c). This depth range provides the bulk of the observed azimuthal anisotropy through simple shear between the rheologically contrasting lithosphere and asthenosphere, and thus approximately represents the bottom of the lithosphere (Yu et al., 2015c). Such a depth of the anisotropic layer is consistent with our tomographic results (Figs. 6 and 7).

Thermal Status

High seismic velocities are usually found in Archean lithosphere because of its iron-depleted composition (Jordan, 1978). Thus the observed low-V zone beneath the ORZ, which occupies part of the Proterozoic orogenic belts, could reflect lateral contrast in chemical composition. However, this model requires low velocities in the entire orogenic belts rather than only beneath the ORZ. The spatial correspondence between the observed low-V zone and the ORZ suggests that the low-V zone has a thermal origin. The lithosphere of old cratons is usually characterized by low temperature and high-V anomalies in the tomographic images, whereas mature rift zones such as the Afar and Kenya rifts usually exhibit high temperature and significant low-V anomalies in the upper mantle (Zhao, 2015). Forward modeling shows that an increase of temperature by 100 °C can result in a Vp anomaly of –0.75% ± 0.15% at 200 km depth and –0.23% ± 0.05% at 800 km depth (Cammarano et al., 2003), although different values have also been proposed (e.g., Goes et al., 2000). Assuming that the resulting velocity anomalies are mainly caused by temperature changes, the prominent localized low-V anomaly with an amplitude of about –1% beneath the ORZ corresponds to a temperature anomaly of ∼100 °C. The low-V anomalies (<–0.5%) beneath the Limpopo belt may result from chemical modification of the mantle during magmatic emplacement, which could be related to the Bushveld Complex (James et al., 2001).

The observed low-V anomalies within the shallow lithosphere of the Congo craton may be related to the granitic-gneissic Quangwadum Complex, which is a Paleoproterozoic basement of Congo craton affinity (Hanson, 2003; Singletary et al., 2003). Mesoproterozoic granite complexes have been proposed to explain the observed low-density anomalies within the crust of the Congo craton based on a recent joint inversion of receiver functions and gravity data (Yu et al., 2015a). Consistent with tomographic studies (James et al., 2001) of the adjacent Kalahari craton, there is no evidence to support a significant low-V asthenospheric layer beneath the Archean keel.

There is a significant low-V zone in the depth range of 300–500 km situated at the northern edge of our study area (Fig. 6). In the vicinity of this area, xenolith-defined temperature profiles suggest that the Congo craton lithosphere experienced melt-induced metasomatic episodes ca. 32 Ma (Batumike et al., 2009). Numerical modeling shows that for a cratonic lithosphere, long-term processes related to metasomatism and fertilization associated with magma percolation may modify the composition of the lithosphere through time (Manatschal et al., 2015).

Decompression Melting and Rifting Implications

The localized low-V anomaly beneath the ORZ could be a result of partial melting in the lithospheric mantle, which can be attributed to increased temperature, reduced pressure, or presence of volatiles (Shankland et al., 1981; Aivazpourporgou et al., 2015). The most commonly cited source leading to an increased mantle temperature is the presence of a mantle plume (Zhao, 2015). However, the expected significant low-V anomalies associated with lower-mantle plumes beneath the ORZ are not observed (Figs. 6 and 7). Instead, our results support the existence of an estimated –1% Vp anomaly in the upper asthenosphere beneath the ORZ. The absence of plumes in the upper mantle and mantle transition zone beneath the ORZ is also supported by the observed normal thickness of the mantle transition zone (Yu et al., 2015b), the lack of elevated mantle conductivity (Khoza et al., 2013), the lack of rift-related anomalies in the dominantly northeast-southwest mantle anisotropy (Yu et al., 2015c), and the localized uplift of the Moho (Yu et al., 2015a).

A possible cause of the observed low-V anomalies beneath the ORZ is decompression melting. A recent geodetic study (Malservisi et al., 2013) shows that the South African block is probably rotating clockwise relative to the African continent along the Damara belt with an estimated spreading rate of ∼0.5–1.5 mm/yr in the vicinity of the ORZ. Such a relative motion can exert a trans-extensional force to the ORZ, leading to mechanical stretching of the lithosphere. Mechanical stretching can in turn induce lithospheric thinning which will reduce the overlying pressure, enhancing decompression melting (Corti et al., 2003; Acocella, 2014). The existence of lithospheric decompression melting beneath the ORZ has also been suggested based on the observed relatively lower subcrustal mantle density from inversion of gravity data (Yu et al., 2015a). In addition, the high heat flow observed within the ORZ may have a mantle origin, as suggested by a 3-D inversion of aeromagnetic data (Leseane et al., 2015).

Our tomographic images show that the lithosphere beneath the ORZ, which is embedded within Proterozoic orogenic belts between the Archean Congo and Kalahari cratons, is thinner than that of the adjacent cratons. The cratonic roots are outlined in the tomographic images by high-V anomalies, which extend to 250–350 km depths beneath the Congo and Kalahari cratons. Our observations provide the most convincing evidence thus far that the southern boundary of the Congo craton extends to near the ORZ. The localized low-V anomaly observed in the upper asthenosphere is most likely caused by decompression melting from lithospheric stretching, and corresponds to an estimated temperature anomaly of 100 °C. Our results suggest that decompression melting is induced during the earliest stage of continental rifting. Results from this teleseismic P-wave tomography study of the ORZ, when combined with those from recent crustal and mantle structure and anisotropy investigations, indicate that the ORZ is characterized as a passive rift initiated in a preexisting zone of lithospheric weakness, probably due to relative movements of rigid cratonic blocks.

Equipment and logistical support from the Incorporated Research Institutions for Seismology (IRIS) Program for Array Seismic Studies of the Continental Lithosphere (PASSCAL) Instrument Center is greatly appreciated. We thank the Data Management Center of IRIS for archiving the SAFARI seismic data sets used in this study. We are grateful to Angela M. Reusch of PASSCAL and Keletso Kaisara of the University of Botswana for field assistance. Constructive comments from two anonymous reviewers and the Associate Editor significantly improved the manuscript. The field and data analysis work was supported by the U.S. National Science Foundation under award EAR-1009946 to S.S.G. and K.H.L., and was partly funded by the National Natural Science Foundation of China through grant 41606043 to Y.Y., and the National Program on Global Changing and Air-Sea Interaction through grant GASI-GEOGE-05.

1 Supplemental Figures S1–S12. Please visit or the full-text article on to view the Supplemental Figures.
2Supplemental Data. Travel-time residuals. Please visit or the full-text article on to view the Supplemental Data.