Lateral Velocity Gradients in the African Lower Mantle Inferred From Slowness Space Observations of Multipathing

Large low‐velocity provinces (LLVPs) are hypothesized to be purely thermal features or possess some chemical heterogeneity but which exactly remains ambiguous. Regional seismology studies typically use travel time residuals and multipathing identification in the waveforms to infer properties of LLVPs. These studies have not fully analyzed all available information such as measuring the direction and inclination of the arrivals. These measurements would provide more constraints of LLVP properties such as the boundary velocity gradient and help determine their nature. Here, we use array seismology to measure backazimuth (direction) and horizontal slowness (inclination) of arriving waves to identify structures causing multipathing and wavefield perturbation. Following this, we use full‐wavefield forward modeling to estimate the gradients required to produce the observed multipathing. We use SKS and SKKS data from 83 events sampling the African LLVP, which has been extensively studied providing a good comparison to our observations. We find evidence for structures at heights of up to 600 km above the core‐mantle boundary causing multipathing and wavefield perturbation. Forward modeling shows gradients of up to 0.7% δVs per 100 km (0.0005 km s−1 km−1) can produce multipathing with similar backazimuth and horizontal slowness to our observations. This is an order of magnitude lower than the previous strongest estimates of −3% δVs per 50 km (0.0044 km s−1 km−1). As this is lower than that predicted for both thermal and thermochemical structures, lateral velocity gradients capable of producing multipathing are not necessarily evidence for a thermochemical nature.


Introduction
Large low velocity provinces (LLVPs) are roughly antipodal, low-velocity features of the lower mantle located beneath Africa and the Pacific and are surrounded by high-velocity material hypothesized to be slab remnants (Bijwaard et al., 1998;Grand, 2002;Grand et al., 1997), shown in Figure 1. Since first observed, LLVPs have remained enigmatic features of the lower mantle with their origin, composition, and therefore their influence remaining uncertain.
suggests that LLVPs are influential on whole Earth dynamics. Despite being very significant for our understanding of global dynamics, many properties of the LLVPs are still unknown and there are several hypotheses of their origin. These hypotheses can be approximately split into those where LLVPs are purely thermal features and those in which they are chemically distinct relative to the surrounding mantle (Garnero et al., 2016). For a purely thermal feature, a common hypothesis is that LLVPs are a cluster of plumes (Schubert et al., 2004) which appear as one large slow feature because of the inherent resolution limitations from seismic tomography (Bull et al., 2009;Davies et al., 2012;Davies, Goes, & Lau, 2015;Ritsema et al., 2007). The thermochemical origin hypothesis requires a source of material chemically unique to the current lower mantle either from the primordial Earth or material that has accumulated over geological time. Material from the primordial Earth is hypothesized to start as a basal layer of material that is swept into piles forming the LLVPs. Mechanisms for the origin of this base layer include a basal magma ocean (Labrosse et al., 2007), accumulation of dense melts (Lee et al., 2010), or an ancient, iron-enriched crust which was then subducted and is stable at CMB conditions (Tolstikhin & Hofmann, 2005). This basal layer could then have been swept into piles observed as LLVPs, which has been shown numerically (Tackley, 1998) and experimentally (Davaille, 1999). Alternatively, they could have accumulated over geological time as subducted lithosphere in the lower mantle (Christensen & Hofmann, 1994;Hirose et al., 1999Hirose et al., , 2005 which is swept into piles, forming the LLVPs (Mulyukova et al., 2015;Tackley, 2011). However, there is some question of the feasibility of producing negative velocity perturbations (Deschamps et al., 2012) and for the slab material to accumulate at the same rate as it is stirred into the mantle (Li & McNamara, 2013).
Depending on the origin of the LLVPs, our understanding of how the Earth evolved from its primordial state changes. If LLVPs are a short-lived cluster of mantle plumes, they do not need to exist in early Earth history. If they are long-lived piles of primordial Earth remnants, their formation and survival would need to be accounted for. Constraining the origin of LLVPs therefore has implications for our understanding of the Earth's history as well as whole Earth dynamics.
To reduce the number of hypotheses, there has been a focus on determining whether LLVPs are purely thermal or thermochemical features. Their relative density could provide constraints but conflicting observations have suggested both higher and lower relative density (Ishii & Tromp, 1999;Koelemeijer et al., 2017;Lau et al., 2017). Anticorrelation of S wave velocity and bulk sound speed (Masters et al., 2000;Su & Dziewonski, 1997) is commonly used as evidence for compositional heterogeneity for LLVPs, but this has also been interpreted as the presence of postperovskite (Davies et al., 2012;Koelemeijer et al., 2015). The presence of strong lateral velocity gradients has been attributed to a thermochemical origin (Ni et al., 2002;, but these gradients can also be replicated with purely thermal structures (Davies et al., 2012;Schuberth et al., 2009).
Most of these studies use observations or constraints from seismological studies. Seismic tomography provides global, broad observations of LLVP location, morphology and relative velocity (e.g., French & Romanowicz, 2014;Grand, 2002;Grand et al., 1997;Koelemeijer et al., 2015;Ritsema et al., 2011;Simmons et al., 2010). The agreement of the long wavelength structure of LLVPs in tomography models shows they are a result of lower-mantle structure and not the different datasets or methodologies used (Lekic et al., 2012). In addition to these global observations, regional seismology studies combine travel time residuals, multipathing observations in the waveform and forward modeling to recover the location,  (French & Romanowicz, 2014) with an isosurface of −1% δV s shown in red and an isosurface of +1% δV s in blue. The isosurface is plotted below 80% of the Earth's radius (5,097 km, 2,205 km above the CMB). (b) Multipathing at LLVP boundaries. As the wavefront moves over a strong lateral velocity gradient, different parts travel at different speeds and arrive at the stations at different times as two distinct arrivals (1). The gradients can cause the wave to diffract and the structure can cause the wave to refract as it passes through it. As a result, multipathed arrivals can arrive from different directions and inclinations (2). gradient and inclination of LLVP boundaries (e.g., Frost & Rost, 2014;He & Wen, 2009;He et al., 2006;Ni et al., 2002;Ritsema et al., 1998;Roy et al., 2019;Sun & Miller, 2013;. Multipathing occurs when a wavefront is incident on a strong lateral velocity gradient that causes the wavefront to move at different speeds and arrive at a recording station with different travel times as two arrivals. In addition to this, the boundary structure causes the wave to diffract and the structure causes the waves to refract as they pass through it, so the multipathed arrivals arrive from different directions and inclinations as well as arrival times. Figure 1 illustrates the multipathing phenomena at LLVP boundaries and how they can be observed at the surface. LLVP boundary studies using travel time residuals and waveforms are common and, from their observations, have estimated the gradients at the boundaries of LLVPs to range from 3% δV s per 50 km (0.0044 km s −1 km −1 ) (Ni et al., 2002) to 2% δV s per 300 km (0.00048 km s −1 km −1 ) (Ritsema et al., 1998) Table 1 for published estimates of African LLVP S wave velocity gradients). Combining travel time residuals, multipathing identification, and forward modeling to observe and infer the properties of structures is well established and has been applied to a variety of structures (Silver & Chan, 1986;Sun et al., 2019Sun et al., , 2010Sun et al., , 2017 and algorithms developed to identify multipathing automatically in the waveforms (Sun et al., 2009;Zhao et al., 2015). Although regional seismology studies only use the waveform to infer the effects of deep Earth structure on the wavefield, they do not analyze all information available such as the direction and inclination of the arrival. Current observations have not been sufficient to constrain LLVP properties and therefore their composition, origin, and influence remain ambiguous. Both purely thermal and thermochemical structures can replicate properties such as the strong gradients, velocity reduction, morphology, and anticorrelation between S wave velocity and bulk sound speed (Davies et al., 2012;McNamara et al., 2010;McNamara & Zhong, 2004Schuberth et al., 2009;Tackley, 1998). Because current seismic observations are not enough to constrain LLVP properties, we explore what can be inferred from analyzing the direction and inclination of the arrivals.
This study uses array seismology to measure the backazimuth (direction) and horizontal slowness (a proxy for inclination) to identify multipathing and regions of diffraction and refraction. This is applied to data sampling the lower mantle beneath Africa, where several studies have identified multipathing and sharp travel time residuals (e.g., Ni et al., 2002;Sun et al., 2009;Wen et al., 2001). Different frequency bands are used to hypothesize differences in the African LLVP boundary structure such as gradient, depth, and inclination. Using these observations, we estimate the gradients capable of producing multipathing with similar backazimuth and horizontal slowness deviations as our observations.

Slowness Vector Grid Search and Beamforming
To measure the backazimuth and horizontal slowness, we search over a range of slowness vectors each with its own backazimuth and horizontal slowness and use beamforming (Rost & Thomas, 2002) to measure the power of the coherent signal. If there are multiple arrivals, we detect high-power values at different backazimuths and horizontal slownesses. The results are referred to as θ-p plots as they describe how the power of coherent signal varies with backazimuth (θ) and horizontal slowness (p). Figure 2 shows examples of clear, possible, and null multipathing observations in the θ-p plots and clear example of multipathing in the  Ni and Helmberger (2003c) −3% per 100-150 km 0.0022-0.0015 Ni and Helmberger (2003a) −3% per 50 km 0.0044 Sun and Miller (2013) −3.5% per 200 km 0.0013 Ritsema et al. (1998) −2% per 300 km 0.00048 This study −0.7% per 100 km 0.00050 Note. The gradients for km s −1 km −1 were calculated using the V s value for PREM (Dziewonski & Anderson, 1981) at the CMB.

10.1029/2020GC009025
Geochemistry, Geophysics, Geosystems waveforms. The analysis is conducted within a time window selected from visual inspection of record section, typically on the order of tens of seconds. Information such as the time windows, stations, measurements, and multipathing identification are in the supporting information.
Most array techniques assume that energy propagates as a plane wavefront (Rost & Thomas, 2002). If the array aperture is small, this assumption holds and the effect of a curved wavefront is negligible. We use data from the Kaapvaal array , which has a large aperture (spread over ≈20°in northwest-southeast orientation) so the plane wave assumption breaks down and can contribute to some deviation from the predicted backazimuth and horizontal slowness.
We alter the travel time calculation of beamforming to account for a circular wavefront given a backazimuth and horizontal slowness ( Figure 3). To calculate the travel times of a circular wavefront moving over a spherical Earth from event to station locations, the radial distances are calculated using the Haversine formula. This distance is then multiplied by an angular slowness value in s/°. From these estimates, the traces are shifted, stacked, and the power of the coherent signal estimated. To search over backazimuth, the event is relocated keeping the epicentral distance constant. From this new location, the radial distance to each station is calculated relative to the mean distance and the travel times calculated (see supporting information).
We test our correction on synthetic data arriving from a known backazimuth and horizontal slowness (see supporting information). We find that our correction reduces the backazimuth deviation from 2.37°to 0.40°and the horizontal slowness deviation from 0.20 s/°to 0.03 s/°.

Multipathing Identification and Slowness Vector Measurements
Multipathed arrivals are identified as power maxima separated in backazimuth and horizontal slowness and with a power value above the background noise and at least 10% of the maximum power value ( Figure 2). We calculate the orientation of the locus between the multipathed arrivals clockwise from north where multipathing is identified. To calculate this, the locations of the multipathed arrivals in the θ-p observation are recorded and the angle of the vector connecting the two points from north is calculated. This angle is then rotated by 90°as the locus is orthogonal to the vector connecting the two points.
Our data set includes SKKS phases (section 2.4) at distances where other phases such as S3KS could arrive at similar times and horizontal slownesses, which make it challenging to identify multipathing. For SKKS observations where multipathing could be present, we analyze synthetics generated using SYNGINE and the 1-D model prem_i_2s Krischer et al., 2017) as an estimate of the relative power of SKKS and S3KS. If there is any power for an S3KS arrival in the synthetic θ-p plots and there are multiple arrivals in the recorded data, the observation is labeled as "possible" multipathing. See supporting information for more details.
In addition to identifying multipathing, several measurements can be made from each observation. The backazimuth residual (ΔΘ) between the observed (Θ observed ) and the backazimuth predicted by the great circle path between the event and mean station location (Θ predicted ) is given by ΔΘ = Θ observed − Θ predicted . The horizontal slowness residual (Δp) between the observed (p observed ) and the PREM (Dziewonski & Anderson, 1981) predicted horizontal slowness (p predicted ) is given by Δp = p observed − p predicted . The vector from the predicted location to the observation location in the θ-p plot is recorded as a measure of the direction and strength of the perturbation the wave has experienced. Figure 4 illustrates the meaning of this vector residual, locus between the arrivals and visualizes backazimuth and horizontal slowness deviations.

Frequency Analysis
To analyze the frequency dependence of multipathing and its wavefield effects, the data are filtered in five frequency bands and analyzed separately (frequency bands: 0.07-0.28 Hz, 0.10-0.40 Hz, 0.13-0.52 Hz, 0.15-0.60 Hz, 0.18-0.72 Hz, 0.20-0.80 Hz) each with a width of two octaves. The frequencies will affect the size of the Fresnel zone, which gives an approximation of the area contributing to the observation. For both the main and multipathed arrival to have enough power to be observed, there needs to be a significant enough velocity change over the Fresnel zone. The frequency variation of multipathing could be indicative of differences in sharpness, depth, or inclination between boundaries. Fresnel zones for each frequency band were calculated at the CMB using velocity values from PREM (Dziewonski & Anderson, 1981), details in the supporting information.

Data and Preprocessing
SKS and SKKS data ( Figure 5) from events located 70°to 140°away from the center of the array and with magnitudes between 5.5 and 7.5 recorded at the Kaapvaal array are used to analyze Africa LLVP boundary structure ( Figure 5). We deconvolve the instrument response, remove the mean amplitude, taper, and apply a band-pass filter between 0.05 and 1.0 Hz (period of 1-20 s) for visual inspection. The horizontal components are rotated to radial and tangential components for clear SKS and SKKS identification. Following this, the signal-to-noise ratio (SNR) is estimated in a 70 s time window around the predicted arrival time and used to roughly sort the data into traces that should be kept (SNR > 3), removed (SNR < 2.5) and could be used (2.5 < SNR < 3). Events with more than 10 traces sorted into "keep" or more than half between the "keep" and the potentially usable bins were sorted by hand after visual inspection of the record section aligned on
The frequency bands we use are limited by the station spacing of the array. If the interstation spacing is too large, spatial aliasing could occur in the θ-p plot and be misidentified as multipathing. The Nyquist criterion for the station spacing of each frequency band is used to limit the frequencies used. The lower frequencies will likely have higher amplitudes and influence the stacking significantly more than the higher frequencies, so we only limit the lower-frequency cutoffs for the frequency analysis using this criterion. The lower-frequency cutoff is limited to 0.20 Hz.

Noise Reduction Techniques
Multipathed arrivals could arrive with a lower SNR and stack to a similar power as incoherent signal at other backazimuths and horizontal slownesses. To aid multipathing identification, several techniques to improve the SNR of coherent arrivals are implemented. We use phase weighted stacking (Schimmel & Paulssen, 1997), F statistic (Blandford, 1974) and deconvolve the array response function (ARF) using the Richardson-Lucy deconvolution method (Lucy, 1974;Richardson, 1972) as done in previous studies (Maupin, 2011;Picozzi et al., 2010). These are detailed further in the supporting information with examples of their effectiveness. We use the outputs of all these methods to identify multipathing in the data with criteria for clear, potential, and no multipathing explained in section 2.2. Measurements of horizontal slowness and backazimuth deviations are taken using the phase-weighted (Schimmel & Paulssen, 1997) stack points as they most consistently have lower noise than the other methods.

Subarrays
To better constrain the location of multipathing and its wavefield effects, the available stations in the Kaapvaal array are grouped into subarrays. Data from all available stations are also analyzed. We group the traces using their waveform properties, backazimuths, and epicentral distances. We accept that we are adding our own bias to the observations by grouping the subarrays this way. Whole array observations are used to identify multipathing, but because the large area of the combined Fresnel zones of the Kaapvaal array, not used to analyze backazimuth and horizontal slowness deviations. The 317 different subarray geometries were used; stations for each subarray are given in the supporting information.

Method Strengths and Limitations
Other studies have developed a method to automatically detect multipathing in the waveform (Sun et al., 2009). In comparison to this method, there are several limitations and advantages. The multipathed arrivals need to be present in enough traces to stack coherently and produce clear arrivals on the θ-p plot. Arrivals of similar slowness may not be resolved as separate arrivals. On the other hand, noisier traces can be used because the stacking methods improve the SNR. The observations themselves also allow measurements of backazimuth and horizontal slowness deviations, which can be used to analyze structures affecting the wavefield.

Multipathing
This section describes our multipathing observations and discusses the frequency dependence (section 3.1) and spatial variation (section 3.2) with interpretations of possible boundary locations. Clear multipathing is observed in 16% of our whole array observations and 6.6% of our subarray observations. Figure 6 shows the variability of multipathing with different frequency bands. Some observations show clear multipathing within specific frequency bands while others in all frequency bands, which could be due to

10.1029/2020GC009025
Geochemistry, Geophysics, Geosystems differences in boundary nature such as the velocity gradient, inclination, or depth. As explained in section 2.3, to observe multipathed arrivals, enough of the Fresnel zone needs to sample different velocities. This requires that the lateral velocity gradient needs to be sufficiently strong and sufficiently sampled by the wavefield. Analyzing the power spectra of data with and without multipathing, using Welch's method (Welch, 1967), shows evidence for an increase in the power of higher frequencies (Figures S10-S14). We hypothesize that this is caused by the focusing of higher frequencies due to diffraction from lateral velocity gradients. Further work is needed to constrain the exact relation between velocity gradient and frequency content and other possible causes, which is not the focus of this study.
Observations of multipathing at different frequencies could be due to differences in wavelength. Multipathing observations at higher frequencies could be indicative of strong velocity gradients, while multipathing at low frequencies is indicative of a significant velocity change over a wider boundary. If the boundary is at an angle to the incidence of the wave, the boundary will not be sampled for as long and appear smoother.
Sampling boundaries at different depths could cause frequency variation in our observations due to changes of wavelength with velocity. At the same depth, and therefore the same 1-D velocity, the boundaries need to have different gradients or inclinations for multipathing to occur at different frequencies. At different depths, the boundaries could be the same sharpness and inclination but observed at different frequencies due to different Fresnel zone sizes.
The size and station density of the array could contribute to the frequency variation. Larger, denser arrays will be sensitive to a larger area and will record multipathed arrivals in more waveforms. Lower frequencies with larger Fresnel volumes are more sensitive to weaker velocity gradients, but the weaker gradients may mean that the multipathed arrivals will have a smaller amplitude also. Whole array observations should have more multipathing observations at lower frequencies (Figure 7) because weaker multipathed arrivals will be recorded in more waveforms and stack to an observable power.

Spatial Analysis
Our observations show that multipathing is not limited to one region and occurs in different frequency bands depending on the region. In Figure 8, the loci and the tomography velocity contours for both whole and subarray observations align well to the east of Africa (25°S, 32°W) with a boundary trending northwest-southeast which then curves to trend approximately west-east as the boundary moves southward. Contours of S40RTS (Ritsema et al., 2011) and pierce points are mapped at 2,400 km depth because for several paths this depth has an increase in lateral velocity gradient (e.g., Figure 12) and features at this depth provide possible explanations for most observations. It is possible, due to the 3-D nature of the structure, structures at other depths could explain our observations. However, we have tested several tomographic, and therefore data-based, models at different depths and found this to be the most satisfactory. In these regions, multipathing is observed in all frequency bands over both whole and subarray observations. The range of frequencies could be interpreted as an LLVP boundary being sampled at several depths or a boundary with both a strong lateral velocity gradients and a significant velocity change.

Geochemistry, Geophysics, Geosystems
The circular low-velocity feature to the southeast of Africa (35°S, 30°W) marked by −1.5%δV s velocity contour aligns well with the loci in the area. Multipathing is observed at a range of different frequencies here with arguably more observations at frequencies above the 0.15-0.60 Hz band. Observing multipathing at higher-frequencies imply sampling of a relatively sharp gradient and observations in a broad frequency range suggest large and sharp velocity changes or sampling boundaries at several depths.
To the west of Africa (25°S, 15°W), particularly in the subarray observations, there is a lot of scatter in loci orientations and multipathing is mainly observed in the higher-frequency bands. The scattered loci are possibly due to the waves traveling through the body of the LLVP boundaries and sampling boundaries at several depths. Depending on the depth, the boundaries could have different orientations, therefore leading to scattered loci. Observing multipathing in higher-frequency bands could be due to strong lateral velocity gradients or the depths the boundaries have been sampled.
There are significant differences between the whole and subarray observations. Whole array observations will be sensitive to a larger area meaning small structures affecting a small part of the wavefield may not be resolved. This could explain why more multipathing is observed to the west in subarray observations. Where multipathing is observed in whole array observations only, the velocity gradient may not be sampled for long enough along raypaths from the event to the subarrays. The multipathed arrivals would then not arrive with enough amplitude to stack to an observable power. Studies using travel time and waveform observations have reported a boundary to the southwest of Africa with an approximate northwest-southeast strike (Ni et al., 2002). The orientation of the locus of our multipathed arrival in this region approximately agrees (Figure 8) supporting these previous results.  (b) show 1-D ray path pierce points at 2,400 km depth (≈500 km above the CMB) for events showing clear multipathing for whole array and subarray observations, respectively. The size and color of the circles correspond to the frequencies at which multipathing is observed. The locus between the arrivals is marked for each frequency to represent the approximate orientation of the boundary causing the multipathing. Subfigures (c) and (d) show the pierce points at 2,400 km depth for clear (red) possible (orange) and no (blue) multipathing for whole and subarray observations, respectively. Velocity contours are shown at 2,400 km depth from tomography model S40RTS (Ritsema et al., 2011). 10.1029/2020GC009025 Geochemistry, Geophysics, Geosystems Sun et al. (2010) find evidence for a mantle plume in the midmantle of this region too. We do not find evidence for this, most likely because of resolution and sampling limitations.
To further explore the spatial distribution of multipathing, we compare the locations of clear, possible, and no multipathing observed at any frequency (Figure 8). Multipathing is not limited to one region and the pierce points of clear multipathing are very close to pierce points that show no or possible multipathing. Our interpretation is the boundary structure needs to be sampled in a specific way for the multipathed arrivals to arrive with observable amplitudes.

Seismic Anisotropy
There have been several studies analyzing seismic anisotropy in the region of this study (e.g., Cottaar & Romanowicz, 2013;Ford et al., 2015;Lynner & Long, 2014;Reiss et al., 2019;Wang & Wen, 2007a). Shear wave splitting could complicate the waveforms and be misinterpreted as multipathing. Therefore, we measure SKS splitting in splitting time and direction of the fast axis, remove the measured effect and repeat the analysis for a selection of events. After the anisotropy correction, we still observe multipathing. For low SNR events, correcting for anisotropy reduced the quality of the observation. Since anisotropy alone is not the cause of observed multipathing and can reduce the quality of some observations, we do not correct for shear wave splitting.

Slowness Vector Residuals
This section analyses the slowness vector deviations and how they vary spatially. We focus our spatial analysis on the full slowness vector deviation from the predicted to the observed arrival in slowness space (Figure 4). Descriptions of backazimuth and horizontal slowness deviations are given in the supporting information. When spatially analyzing these deviations, the pierce point location is moved to match observed backazimuth and horizontal slowness.
Backazimuth residuals shows little variation between frequency bands (supporting information). The majority of the observations lie between 8°and −14°and maximum values are 10°to −22°for positive and negative deviations, respectively. There are more negative residual observations with on average ≈64% negative residual observations compared to 36% positive. This is possibly due to heterogeneous sampling from limited event station configurations.
The horizontal slowness deviations vary little with frequency, with most observations lying between 1.2 s/°a nd −1.0 s/°(supporting information). Outliers are present in these observations, but show no clear pattern and range from a maximum of 2.1 s/°and a minimum of −1.6 s/°. Like the backazimuth residuals, the observations are not evenly distributed about 0 s/°with 60% positive residuals and 40% negative. This variation could be due to the dominantly slow mantle structure beneath Africa causing them to refract and arrive at a shallower angle.
The magnitude of the slowness vector deviations does not vary greatly with frequency with slightly more high-magnitude deviations at higher frequencies and with minimum and maximum observed values from less than 0.1 to 2.1 s/°(supporting information).

Spatial Analysis Slowness Vector Deviation
The full slowness vector deviation is a vector from the predicted arrival in the θ-p plot to the observed arrival. The azimuth of the vector indicates the direction of perturbation, and the length is indicative magnitude. This vector combines the backazimuth and horizontal slowness perturbations giving a clear picture of how the wavefield is being affected. Figure 9 shows how these vectors vary spatially.  (Figure 4). The slowness vector describes the full perturbation of the wavefield essentially combining the information from backazimuth and horizontal slowness. The contours from S40RTS (Ritsema et al., 2011) and the pierce points are marked at a depth of 2,400 km to outline potential structures contributing to the observations. The frequency band used is from 0.13 to 0.52 Hz. The pierce points have been relocated according to the observed backazimuth and horizontal slowness.

Geochemistry, Geophysics, Geosystems
The radial pattern and magnitude of the vectors around the circular feature southeast of Africa (35°S, 30°W) support our interpretation that this structure is the cause of our observations. Northeast of this region, the azimuths change their orientation to be approximately orthogonal to the velocity contours striking northwest-southeast. The vector residuals west of Africa (25°S, 15°W) are more scattered than in other regions but generally have an azimuth pointing away from the array and arrive at a shallower inclination. The scattered vector residuals, the scattered loci, and the presence of multipathing in this region suggests that the wavefield is being affected by several boundaries at different depths. The trend of slowness vectors pointing away from the array suggest that they are being refracted by the body of the LLVP.
East of Africa (25°S, 40°W) the vectors have opposite azimuths shown by the color change from red to green in the vector heads. The paths of these waves suggest they may not sample the LLVP boundary. The arrivals sampling this region arrive with negligible backazimuth residuals, but have opposite horizontal slowness residuals ( Figure 10). We hypothesize the cause of these observations are adjacent fast and slow structures causing the wavefield to vertically refract. The location of fast structures relative to the LLVP boundary at the core-mantle boundary in tomography model SEMUCB-WM1 (French & Romanowicz, 2014) aligns well with the transition (Figure 10), implying that these structures could be the cause of our observations. Given the size of the subarrays and the size of the Fresnel zone at these frequencies, it is possible this fast structure is causing the waves to refract and arrive at a steeper inclination with negligible backazimuth deviation.
Previous studies have analyzed similar regions and show some evidence for similar structures. Sun et al. (2009) analyze regions of the lowermost mantle similar to areas where we find boundaries between slow and fast structures and a quasi-circular structure. Using their multipath detector method with S diff data, they identify a region with strong gradients southeast of Africa similar to our hypothesized boundary in Figure 10. Their travel time residuals transition from negative to positive over this region supports our interpretation of a transition from slow to fast. Another Sun et al. (2009) event shows evidence for smaller-scale structure southeast of the Kaapvaal array with a similar structure and approximate location as our observed circular structure. Mega ULVZs have been shown to cause backazimuth deviations  and are a possible explanation of our observations. However, a 20 km tall ULVZ with a 20% δV s reduction lead to travel time residuals of 0.7 s, which is well below the observed difference between multipathed arrivals. We tested travel time contributions of stronger ULVZ models (Yu & Garnero, 2018) and found they also cannot explain our results.

Forward Modeling and Comparison to Tomography Models
In this section, we explore possible conditions for multipathing to be observed in slowness space and compare these to previous studies (Table 1).
Our observations are sensitive to 3-D structure, show frequency dependence, and of relatively high frequencies (up to ∼0.20 Hz). Therefore, forward modeling strategies based on approximations, such as ray theory, are not appropriate to reproduce them. Instead, we must use a full-physics numerical method to capture the finite-frequency effects of wave propagation in a 3-D medium, which is necessarily much more computationally costly. We use SPECFEM3D (Komatitsch & Tromp, 2002a, 2002b to create synthetic data for three Figure 10. Pierce points at the CMB colored with (a) horizontal slowness deviations and (b) backazimuth deviations. Negative contours −1.0%, −1.5%, −2.0% δV s and positive contours 0.5%, 1.0%, 1.5% δV s of tomography model SEMUCB-WM1 (French & Romanowicz, 2014) are shown to highlight the transition from fast to slow structures east of Africa. The events have been relocated so the 1-D paths arrive from the observed backazimuth and horizontal slowness.

10.1029/2020GC009025
Geochemistry, Geophysics, Geosystems earthquakes (section 5.1). These events were chosen as they show the clearest multipathing, and therefore likely sample the strongest gradients along their raypaths. As such, the analysis of these events will provide us with an upper bound of the required seismic velocity gradients to explain the data in this locale. As the modeling is computationally expensive, we limit ourselves to these events and model frequencies up to ≈0.18 Hz. We test the effects of ellipticity and topography and find they have a negligible effect.
The loss of small-scale heterogeneity and reduction of velocity amplitude and gradients in seismic tomography from regularization, smoothing, and limited sampling coverage is well documented (Bull et al., 2009;Foulger et al., 2013;Ritsema et al., 2007;Schuberth et al., 2009). Given the large parameter space of a 3-D structure that could cause multipathing, we take the structure of tomography as an approximation of long-wavelength Earth structure and accept the mentioned limitations. From this starting point, we increase the velocity perturbations and gradients linearly to approximately account for the reduction through tomographic filtering and recreate conditions for multipathing to be observed in our method. Figure 11. Analysis of multipathing for three events in the observed data (top row) with synthetics from models (second row) M1 to (fifth row) M4 in the rows beneath (labeled on the right). For each event, the same frequency bands are used for the observed and synthetic data.

10.1029/2020GC009025
Geochemistry, Geophysics, Geosystems S40RTS (Ritsema et al., 2011) is used as a starting point as the velocity contours shown in the figures in sections 3 and 4 provide possible explanations for our observations. In each model, the velocity perturbations have been amplified at depths greater than 1,000 km with depths shallower than 660 km are unchanged. The transition from the amplified lower mantle to the upper mantle is tapered to avoid artifacts. No crustal model is used in our modeling as tests show no identifiable effect of crustal structure on our observations. Three models are used where perturbations at depths greater than 1,000 km have been doubled (labeled as M2), trebled (M3), quadrupled (M4), and we use S40RTS (Ritsema et al., 2011) with no amplification (M1). Of course, any single mantle velocity model is not a unique fit to the data and many other possible models exist. However, using tomography-based, and therefore data-based, models means we incorporate structures which other observations have already identified.

Gradients of Boundaries
We compare observations of SKS data from events on the 25 May 1997, 28 March 1998, and 28 May 1997 to runs using all models described earlier. Figure 11 shows the θ-p plots of the synthetic data with the observations. For these events, the S40RTS (Ritsema et al., 2011) velocity perturbations are not sufficient to cause detectable multipathing, indicating that stronger gradients are required. In models with stronger gradients, whether multipathing is observed and how similar it is to the observation varies with the event likely due to the different sampling geometry. Synthetic data for the 25 May 1997 in model M3 shows clear multipathing where the relative power and location of the two arrivals are similar to the observation. In model M2, there is no clear multipathing and the location of the arrival is approximately the average of the locations of the observed multipathed arrivals. As the only difference between M2 and M3 is the strength of amplitudes in the lower mantle, we argue that it is the lower mantle structure causing the observed multipathed arrivals in this event.
The 29 May 1997 event shows some weak multipathing in all amplified models in similar locations to the observation, but the arrivals do not have the same relative power in the θ-p plot. This suggests that there is a boundary being sampled, but the gradient in the model is weaker or the path length along the boundary is shorter than that in the data. The 29 March 1998 event shows no multipathing in most of the models except for M4, but this has much weaker multipathing and both arrivals are different to their location in the observation. The strength of the velocity gradient of the boundary or its location in the models is not enough to reproduce the observation.
These varying results are to be expected with the inherent limitations of tomography and show that the models are not representative of the LLVP everywhere. Due to the good agreement between synthetic data from model M3 and real observation for the 25 May 1997 event, we analyze the gradients sampled along this path. Figure 12a shows a cross section through model M3 between 70°and 140°distance along the great-circle path for the 25 May 1997 event. The 1-D raypath to the mean station location from the exit point is shown on the cross section with the lateral velocity gradient, and velocity perturbations sampled are shown in  Figure 12. The largest gradients sampled are not at the CMB but ≈600 km above it, a similar depth to the maximum misfit found by Zhao et al. (2015) in their analysis of waveform broadening and the Pacific LLVP. The path in the cross section samples the edge of a particularly low-velocity region at 4,100 km radius (shown by velocity contour of −4% δV s ), which we hypothesize is the same feature causing the circular contours shown previously. The maximum gradient sampled is 0.7% δV s per 100 km (0.0005 km s −1 km −1 ) about 600 km above the CMB. This is an order of magnitude lower than that found in some previous studies, which we discuss further below.
Although the modeled θ-p observation is similar, the modeled SKS data arrives much earlier than in the observations as shown in Figure 13. The difference in travel times is a reflection of the velocity perturbations sampled, whereas the observation of multipathing is indicative of the gradients sampled. For this example, the gradient sampled over the raypath is sufficient to create comparable multipathing to the observations, but the velocity perturbations are not sufficient to replicate the observed travel time residuals.
In the synthetic waveforms, there is evidence for diffracted phases such as SPdKS. We do not think this is the cause of the second arrival in the θ-p plots. SPdKS is expected to arrive within a narrow slowness window, which the majority of our arrivals are not observed in. Furthermore, multipathed arrivals only appear in the θ-p plots once the velocity gradients have been increased. If SPdKS arrivals are present, they should be observed in the θ-p plots of all models and not just those with amplified gradients. We use a relative stack to isolate SPdKS using data in Figure 13 at distances larger than 119°where SPdKS is most clear and separated from SKS. We find that SPdKS is not visible when SKS is included in the analysis time window (supporting information) and suggest it is because it arrives with lower amplitude. When not including SKS, SPdKS is only visible in the synthetics using PREM (Dziewonski & Anderson, 1981) and not in the real data or synthetics from M3. Because of this, we are confident that the multiple arrivals observed are multipathing and not SPdKS.  (Dziewonski & Anderson, 1981) predicted arrival times for SKS (red line) and SPdKS (blue line) are added. The multipathed arrivals in the observed data and M3 synthetics have a visibly different moveout to SPdKS. We discuss the possible presence of SPdKS further in the supporting information. The modeled waveforms arrive significantly earlier than the observations. We suggest this is a reflection of the velocity perturbations in the model rather than the lateral velocity gradients.
10.1029/2020GC009025 Miller, 2013;Wang & Wen, 2007b). We assume the gradient of the boundary is the main cause of the observed multipathing. As only one of our models matches well with the observation, we only compare the gradient we found to produce multipathing for the 25 May 1997 event with other studies (see Table 1). Gradients up to 0.7% δV s per 100 km (0.00050 km s −1 km −1 ) can produce multipathing for the 25 May 1997 event which is an order of magnitude lower than the strongest estimated gradients of −3% δV s per 50 km (0.044 km s −1 km −1 ) (Ni et al., 2002), though similar to that found by Ritsema et al. (1998) −2% δV s per 300 km (0.00048 km s −1 km −1 ).
The apparent disparity between studies may be explained by the differing sensitivity of the methods used. Our observations analyze coherent signals across the array by stacking many waveforms together and not analyzing them individually. Each measurement is sensitive to a larger region and could lead to boundary structures being sampled for longer, therefore weaker gradients are sufficient to produce multipathing. Most estimates of the stronger gradients used 2-D forward modeling to replicate their observations (Ni & Helmberger, 2003b;Ni et al., 2002), therefore estimating the gradients along the great-circle path. Any travel time delay or multipathing would have to be from in-plane structures and contributions from out of plane structure would not be accounted for. We use 3-D full wavefield modeling, thus accounting for contributions from out of plane structures, which more fully represents the ability of weaker gradients to lead to the same effect.
The presence of strong velocity gradients at LLVP boundaries causing multipathing and sharp changes in travel time residuals is commonly used as evidence for a thermochemical origin of LLVPs (Ni et al., 2002;Ritsema et al., 1998;. We require gradients an order of magnitude lower than previous estimates to produce multipathing similar to our observations. The gradients of 0.7% δV s per 100 km (0.00050 km s −1 km −1 ) are well below those evident in purely thermal models (2.25% δV s over 50 km (0.0032 km s −1 km −1 ) (Schuberth et al., 2009) and 3.5-4.5% δV s per 100 km (0.0025-0.0032 km s −1 km −1 ) (Davies et al., 2012)). This modeling implies that the presence of lateral velocity gradients capable of producing observable multipathing cannot distinguish between thermal and thermochemical LLVPs.

Conclusions
Through measuring the backazimuth and horizontal slowness of SKS and SKKS data sampling the lower mantle beneath Africa, we identify clear multipathing in ≈16% of our whole array observations and 8.0% of our subarray observations. We find evidence for wavefield perturbation from backazimuth deviations of up to 22°and horizontal slowness deviations of up to 1.2 s/°. Spatial analysis of these measurements relative to structure resolved by seismic tomography gives evidence for a circular feature to the southeast of Africa, adjacent fast and slow structures, and an LLVP boundary. This suggests that tomography models are able to recover the shape but not strength of the features explaining our observations. We conduct full wavefield forward modeling to constrain what lateral velocity gradients are needed to observe multipathing in slowness space. We find gradients of up to 0.7% δVs per 100 km (0.00050 km s −1 km −1 ) sampled ≈600 km above the CMB can reproduce our multipathing observations. This is an order of magnitude lower than previous estimates of −3% δV s per 50 km (0.0044 km s −1 km −1 ) (Ni et al., 2002), commonly used to argue for a thermochemical origin of LLVPs. As the gradients we predict are well below the largest estimates for both thermal and thermochemical structures (Davies et al., 2012), we argue that observing multipathing caused by lateral velocity gradients of LLVP boundaries is not necessarily evidence for a thermochemical composition.