Asthenospheric Flow of Plume Material Beneath Arabia Inferred From S Wave Traveltime Tomography

Widespread Cenozoic volcanism in the Arabian Peninsula has been attributed to mantle plume activity and/or lithospheric thinning due to rift‐related extension. However, there is discrepancy between geochemical and geophysical studies about which mechanism is dominant over the other for post‐12 Ma volcanism. Plume signals in some volcanic fields in the Arabian shield are not evident in isotope analyses, but low‐velocity anomalies connected to Afar are found beneath the Arabian shield in tomographic studies and interpreted as asthenospheric flow from the Afar plume. To resolve this contradiction, we investigate the upper mantle beneath the Arabian Peninsula and northeastern Africa by inverting relative S and SKS wave arrival times recorded at dense seismic networks to derive a high‐resolution S wave velocity model. Our results clearly show a narrow, elongated low‐velocity anomaly along the Makkah‐Madinah‐Nafud (MMN) volcanic line beneath the Arabian shield at 100–300 km depth which extends northward to Harrat Ithnayn and Harrat Lunayyir, but most likely not further north. The limited extent of the low‐velocity anomaly and variations in lithospheric thickness of the Arabian shield may explain why geochemical studies did not find plume signals in some harrats. Therefore, the timing and plume signals of volcanism in western Arabia may not be age progressive from Afar. We also find a possible connection between the low‐velocity anomalies beneath Harrat Lunayyir and the MMN line, suggesting that the 2007–2009 seismic swarm may be associated with northward asthenospheric flow of plume material from Afar.


Introduction
The separation of the Arabian plate from the African plate initiated at~30 Ma as rifting systems along the Red Sea and the Gulf of Aden were activated (e.g., Bohannon et al., 1989;Bosworth et al., 2005;Le Pichon & Gaulier, 1988). Around the same time, massive flood basalts were produced in northern Ethiopia and Yemen to form the Ethiopia-Yemen traps (e.g., Hofmann et al., 1997). The precise origin of the Red Sea is still debated, whether the dominant mechanism is passive rifting due to slab pull from the Zagros subduction zone (e.g., McGuire & Bohannon, 1989;McQuarrie et al., 2003), active rifting by upwelling of a hot mantle plume (e.g., Bellahsen et al., 2003;Hansen et al., 2007;Ritsema et al., 1999), or a mix of active and passive rifting, where plume activity initially drives continental breakup, but preexisting lithospheric heterogeneities and plate-boundary forces provide a passive influence (e.g., Courtillot et al., 1999).
A similar debate is ongoing over the origins of Cenozoic volcanism in western Arabia. The Arabian plate consists of Proterozoic basement rocks, Phanerozoic sedimentary rocks, and Cenozoic volcanic fields. The Cenozoic volcanic fields, also known as "harrats," are widely spread over the Arabian shield (Figure 1a), while the Arabian platform is covered by sedimentary rocks up to 10 km in thickness (Chen et al., 2015;Stoeser & Camp, 1985). The Cenozoic volcanic fields may have originated from plume activity and/or lithospheric thinning due to rift-related extension and can be broadly divided into older and younger groups (Bertrand et al., 2003;Camp & Roobol, 1992;Daradich et al., 2003). The older basalts (30-20 Ma) are tholeiitic to transitional in composition and oriented in a northwestern direction parallel to the Red Sea, whereas the younger basalts (12 Ma to present) are transitional to strongly alkalic trending northward (Camp & Roobol, 1992). This had led researchers to attribute the origin of the older basalts to the rifting of the Red Sea, while the origin of the younger basalts remains enigmatic.
Geophysical studies have focused on imaging mantle velocity structure under the Ethiopia-Afar region and the Arabian Peninsula to find the source of the post-12 Ma basalts. Debayle et al. (2001) performed a waveform inversion and found that a low-velocity anomaly beneath Afar is deeply rooted down to the mantle transition zone, inferring the presence of mantle plume upwelling. Sicilia et al. (2008) reached the same conclusion of a deeply rooted low-velocity anomaly beneath Afar using surface wave tomography. Benoit, Nyblade, and VanDecar (2006), Benoit, Nyblade, Owens, and Stuart (2006), and Bastow et al. (2008) found strong low-velocity anomalies down to~400 km depth beneath Ethiopia using teleseismic traveltime tomography, and they asserted that the anomalies may be caused by upwelling of hot material from the lower mantle across the mantle transition zone. Recently, Civiero et al. (2015Civiero et al. ( , 2016 found two low-velocity anomalies down to~800 km beneath Afar and west of the Main Ethiopian Rift using teleseismic P and S wave traveltime tomography. On the other hand, Ebinger and Sleep (1998) proposed that a single large plume could feed volcanism in Africa via channeled flow along thin lithosphere using numerical simulations, which has been also supported by some global tomographic models (e.g., Hansen & Nyblade, 2013;Ritsema et al., 1999). Park et al. (2008) showed narrow low-velocity anomalies at~200 km depth that extend northward along the southern Red Sea coastline and the so-called Makkah-Madinah-Nafud (MMN) volcanic line of post-12 Ma basalts in the Arabian shield using surface wave tomography. They interpreted the anomalies as northward propagation of warm mantle material from Ethiopia. Other regional tomographic models with surface wave measurements also showed a northward trend of low-velocity anomalies at asthenospheric depths, extending from Afar to the Arabian shield along the MMN line Tang et al., 2018Tang et al., , 2019Yao et al., 2017). The Ha'il Arch along the MMN line, formed in pre-Permian times (Bosworth et al., 2005) and reactivated~15 Ma (Coleman, 1993), may play a role as a natural lithospheric channel for the mantle flow. Benoit et al. (2003) and Park et al. (2007) also showed low-velocity anomalies beneath the Arabian shield using teleseismic traveltime tomography. Koulakov et al. (2016) obtained P and S velocity tomographic models utilizing traveltime data and found similar patterns at asthenospheric depths. High attenuation in the upper mantle is found beneath the Arabian shield in attenuation tomography studies (e.g., Al-Damegh et al., 2004;Pasyanos et al., 2009), which may be related to the low-velocity anomalies through compositional variation such as fluid/melt or hot temperature. Some researchers found northward alignments of fast axes of shear-wave splitting at stations installed in Arabia (Elsheikh et al., 2014;Hansen et al., 2006;Wolfe et al., 1999), providing evidence that the observed low-velocity anomaly beneath Arabia in tomographic studies may be caused by northward asthenospheric flow from the Afar plume. Some studies have attributed the formation of harrats in the Arabian shield to local plumes beneath the Arabian shield (Camp & Roobol, 1992), the Arabian platform (Koulakov et al., 2016), or Jordan (Chang & van der Lee, 2011). Unlike aforementioned geophysical studies emphasizing roles of mantle plumes for the post-12 Ma volcanism in Arabia, some geochemists have proposed lithospheric thinning due to extension related to the Red Sea rift as the source of post-12 Ma thermal activity in the Arabian shield (Bertrand et al., 2003;Moufti et al., 2012Moufti et al., , 2013. Bertrand et al. (2003) opposed the plume origin for the harrats because a high 238 U/ 204 Pb ratio (HIMU) pervades the entire Arabian plate from Yemen to Syria and Jordan, which may be more compatible with the assumption of presence of HIMU in the mantle domain of the continental lithosphere. In addition, northward attenuation of the HIMU signature was not found, which should have been observed if there was asthenospheric flow from Afar. Based on isotope data in Harrat Al-Madinah (the northern part of Harrat Rahat), Moufti et al. (2012) argued that the HIMU source existed within the lithosphere and was produced by interaction of the lithosphere with metasomatic fluids. Thus, they suggested that the harrats were formed by decompression melting triggered by lithospheric thinning related to the Red Sea rift. By presenting new 40 Ar/ 39 Ar ages in Harrat Al-Madinah, Moufti et al. (2013) challenged the idea of northward migration of volcanic activity by channeled flow of plume material, because they found that the whole duration of volcanic activity in Harrat Al-Madinah is very similar to that of the entire Harrat Rahat. Altherr et al. (2019) and Konrad et al. (2016) also did not find any plume signals in Harrat Uwayrid and Harrat Hutaymah. However, a study by Duncan et al. (2016) argued for plume sources for the origin of the harrats because a thin lithosphere beneath the MMN line may be attributed to active mantle flow that provides thermal buoyancy for uplift and thermal erosion of lithosphere. By measuring new 40 Ar/ 39 Ar ages, they suggested that Harrat Hutaymah formed rapidly and recently due to a plume-influenced source, and they proposed that the coupled effects of northward asthenospheric flow from the Afar plume and rift-related extension may be responsible for the variability of depth and degree of melting and differences in timing among harrats. Thus, whether asthenospheric flow from the Afar plume is responsible for the formation of harrats in western Arabia is hotly contested, and it is not clear how far plume flow has propagated laterally, or to what extent the harrats have been affected by the asthenospheric flow.
Understanding upper mantle structure and mantle flow below the Arabian Peninsula plays a key role in finding sources for the post-12 Ma volcanism on the surface. To better understand the link between the post-12 Ma volcanism and mantle flow, we derive an S velocity model of the mantle beneath Arabia and northeastern Africa by inverting 17,155 teleseismic S and SKS wave relative arrival times measured from seismic data collected by 407 seismometers at regional seismic networks ( Figure 1b). We show that our model has better resolution in the upper mantle, especially beneath the Arabian shield, compared to previous studies (e.g., Benoit et al., 2003;Park et al., 2007) owing to the dense seismic network operated by the Saudi Geological Survey (SGS). We image linear low-velocity anomalies which we interpret as asthenospheric flow beneath the Arabian shield and discuss its extent of propagation and its relationship with the post-12 Ma volcanism. We also discuss the role of the mantle flow in relation to the seismic swarm in Harrat Lunayyir in 2007-2009.

Data Collection
We collected data for events occurring between 2007 and 2014 from the SGS, and we used data collected by IRIS and GEOFON for the years of 1995, 2000-2002, 2007-2014, and 2017-2018 to include permanent as well as temporary seismic networks (Figure 1b and Supporting Information Table S1), thereby increasing data coverage over Arabia and northeastern Africa. The azimuthal distribution for S phases is predominantly biased to the east of the study region because of earthquakes occurring in Pamir, Tibet, and the western Pacific (Figure 2). To alleviate this azimuthal bias, we included SKS phases, which increased the number of data from the west by including events occurring in subduction zones beneath North and South America. Furthermore, the addition of SKS phases can increase model resolution at greater depths, since crisscrossing S and SKS rays can be found at deeper depths due to their different ray parameters. We collected S phases with epicentral distances of 30°to 90°and SKS phases between 100°and 140°. To ensure high-quality signals, we used a minimum magnitude of M W 5.5. SGS data were obtained from 143 broadband seismic stations located in Saudi Arabia, and IRIS and GEOFON data were from 192 and 72 stations, respectively, both distributed across the Arabian Peninsula and northeastern Africa (Figure 1b). After preprocessing, we band-pass filtered the event data between 0.04 and 0.1 Hz using a zero-phase two-pole Butterworth filter to extract S and SKS phases properly. Relative S and SKS wave arrival times were estimated using 297 and 134 events, respectively ( Figure 2).

Relative Traveltime Residuals
We derive an S velocity mantle model for the study region by inverting S and SKS wave relative traveltime residuals that are calculated using the multi-channel cross-correlation (MCCC) method (VanDecar & Crosson, 1990). We utilize the IASP91 model (Kennett & Engdahl, 1991) as a 1-D reference model to compute predicted traveltimes. The MCCC method is a semi-automated technique to determine relative arrival times of S or SKS phases recorded at regional seismic networks from teleseismic events (VanDecar & Crosson, 1990). To avoid cycle skipping and time-consuming manual phase picking, we first apply the iterative cross-correlation stack algorithm to automatically align traces with narrow time windows (Lou et al., 2013). In this algorithm, all waveforms are stacked according to predicted arrival times using a 1-D reference model. Then each trace is cross-correlated with the array stack to obtain a lag time at peak correlation, and the arrival time is updated. By repeating this process until a convergence criterion is satisfied, this algorithm makes preliminary picks of arrival times before using MCCC. The MCCC method is then applied to refine the relative arrival times.
The MCCC method proceeds as follows. First, by taking advantage of waveform coherency between waveforms from the same teleseismic event, we find a relative delay time (Δt ij ) that maximizes the cross-correlation function between the ith and jth waveforms (t i and t j ), represented as Δt ij ¼ t i − t j . Since waveforms are usually contaminated with noise, relative delay times are not perfectly consistent (i.e., Δt 13 ≠ Δt 12 +Δt 23 ). Because of this inconsistency, a least squares method is used to optimize the Figure 2. Distribution of earthquakes used for measuring S and SKS traveltimes. White and black circles are the event locations for S and SKS phases, respectively. Gray lines represent 30°interval distances from the center of study region, and the red lines indicate plate boundaries (Bird, 2003). cross-correlation-derived delay times subject to the constraint that the mean of the optimized relative delay times is zero (∑t i ¼ 0). This procedure thus removes the background mean of the traveltimes from the relative traveltime residuals. We used a cross-correlation coefficient threshold of 0.8 to obtain high-quality relative delay times, and we removed residuals larger than ±10 s as outliers. Based on the described methods, we estimated 13,441 S phase residuals from 297 events and 3,714 SKS phase residuals from 134 events.

Crustal Correction
Body waves from teleseismic events arrive at near-vertical incidence angles, causing shallow structure to have limited resolution in body-wave tomography. To avoid contamination from shallow structure, a "station correction" term can be added to the tomographic inversion as an unknown for each station, which corrects traveltime residuals by removing a common delay time. However, this approach ignores prior knowledge about crustal structure and may cause large velocity anomalies under a station, such as those due to a craton or mantle plume, to be mapped into large station corrections making them disappear from the resulting model (Nolet, 2008). Instead, a priori 3-D crustal model can be used to calculate the traveltime differences between the 3-D crustal model and 1-D reference model and then correct traveltime residuals to remove the crustal effect. Although there are high-resolution regional crustal models for parts of our study region (e.g., Hammond et al., 2011;Kaviani et al., 2020;Tang et al., 2019), there are few regional models that cover the entire area from northeast Africa to Arabia. Therefore, for consistency we use the CRUST1.0 model (Laske et al., 2013), a global crustal model made from a compilation of various sources of information on a 1°× 1°grid, for the crustal correction. We extract crustal properties from CRUST1.0 beneath each station to create a corresponding 1-D velocity model.
The traveltime through 1-D crustal structure can be represented as where p is the ray parameter, c is shear-wave velocity, and Δz is layer thickness. However, this sums up traveltimes for segments along a perturbed ray path in the 3-D crustal model (1-D model corresponding to each station) with the same slowness p as in the unperturbed ray path in the 1-D reference model. The traveltime for the same epicentral distance as for the unperturbed ray path (denoted as τ) is obtained by correcting for the extra epicentral distance δΔ: Therefore, the crustal correction between CRUST1.0 and IASP91 can be defined as follows:

Model Parameterization and Inversion Procedure
We use spherical shells of grid points to parameterize the S velocity distribution in our model. We used 16 shells from the surface to 920 km depth. The depth intervals between shells are 30 km between 0 and 120 km, 50 km between 120 and 470 km, and 100 km between 470 and 770 km, with the deepest shell at 920 km depth. The grid point locations for each shell are derived through triangular tessellation of a sphere (Wang & Dahlen, 1995). The center of our grid is located at 21°N, 44°E and extends 25°in all directions from the center point. Each shell contains 3,406 grid points, which are 85 km apart on average at the surface. The total number of grid points is 54,496.
We use the IASP91 model (Kennett & Engdahl, 1991) as a 1-D reference model for the inversion. The ray paths along event-station pairs are calculated using 1-D ray tracing through IASP91 based on ray theory. Differences between observed and predicted relative delay times are arranged in a data vector, d. S velocity perturbations are assigned to model parameters, m. Partial derivatives of delay times with respect to model parameters constitute the sensitivity kernel matrix, G. These three components form a set of linearized equations:

10.1029/2020JB019668
Journal of Geophysical Research: Solid Earth Gm ¼ d: (4) We then define a misfit function as follows: where F is a flattening matrix, I is the identity matrix, and λ F and λ D are flattening and damping weights, respectively. The flattening matrix is a first-order differencing operator between adjacent grid points in the lateral and vertical directions. We minimize Equation 5 by applying the iterative conjugate-gradient algorithm LSQR (Paige & Saunders, 1982) with up to 100 iterations. Optimal regularization values were chosen by analyzing the trade-off between model roughness and data fit (Figure 3). The damping and flattening values of the final S velocity model were 1.5 and 1.0, respectively, resulting in~82% variance reduction. We present various S velocity models with different damping parameters in Figure S1. Although anomalies become weaker or stronger, depending on the damping parameters, the patterns of anomalies remain essentially unchanged, therefore confirming the robustness of anomalies resolved in our preferred model. We also considered event mislocations and structural effects outside of the model domain by including source terms (latitudes, longitudes, depths, and origin times of events) as unknowns for every event in the inversion. Moreover, our study region covers a wide area that includes high-velocity anomalies beneath the Arabian platform as well as low-velocity anomalies beneath the Arabian shield, Afar, and the Main Ethiopian Rift, thereby alleviating the possibility of removing common structure such as cratons or hotspots (Bastow, 2012).

Resolution Tests
We performed several resolution tests to assess the resolving power of our S velocity mantle model. Synthetic data vectors were obtained by multiplying hypothetical models with the sensitivity kernel matrix; resulting structural models are then returned by inverting the synthetic data vectors using the same regularization and parameterization as in the inversion with real data. In all resolution tests, we added Gaussian noise to the synthetic data using maximum values twice as large as the corresponding data uncertainties.
First, we tested a checkerboard model with box sizes and amplitudes of 200 km × 200 km and ±200 m/s, respectively. Depth slices from 100 to 600 km of the recovered model are presented in Figure 4. We find that checkerboard anomalies from 100 to 400 km depth are recovered well under most of the Arabian shield, Israel, Jordan, and Ethiopia, corresponding to regions with densely distributed seismic stations. Smearing effects begin to appear under the eastern Arabian Peninsula at 500-600 km depth. Checkerboard anomalies with larger sizes of 300 km × 300 km are well recovered throughout the study region down to 600 km depth ( Figure S2).
Second, we carried out vertical resolution tests using one horizontal layer of box-shaped anomalies with sizes and amplitudes of 200 km × 200 km and ±200 m/s at 200 and 400 km depth aligned along four transects of interest, which are presented in Figures 5 and 6, respectively. Our results show that the input anomalies are generally well recovered, although vertical smearing effects occur due to quasi-vertically incident S and SKS phases, especially when using box-shaped anomalies at 400 km depth. In Section A-a across Harrat Lunayyir, Harrat Ithnayn, and Harrat Hutaymah (  depth ( Figure 6) show lower resolution with stronger smearing. However, it seems that the depth extent of the resolved anomalies is confined mostly to the deep upper mantle (300-500 km depth) except for the mantle beneath the Arabian platform, thereby indicating that anomalies at 400 km depth are not strongly mapped into artifacts at 200 km depth.
Third, we conducted several structural resolution tests along Cross Sections A-a, B-b, and C-c, which are shown in Figures S3-S5, respectively. For Transect A-a, we investigated the robustness of two low-velocity anomalies beneath Harrat Lunayyir and the MMN line (Harrat Ithnayn) down to 300 km depth separated by a small distance ( Figure S3), which are resolved well with weak vertical smearing. In Section B-b, we performed a structural test with an elongated low-velocity anomaly within 100-300 km depth along the MMN volcanic line connected to a low-velocity anomaly extending to the mantle transition zone beneath Yemen ( Figure S4). The recovered model shows that the input anomalies are relatively well resolved with some smearing given the limitation of teleseismic traveltime tomography. In Section C-c, we tested the resolving power of our model to retrieve a continuous low-velocity anomaly within 100-300 km depth along the Main Ethiopian Rift, linked to a large low-velocity anomaly down to the mantle transition zone beneath Afar ( Figure S5). In the retrieved model, the anomalies are well resolved, except for the southern end of the Main Ethiopian Rift due to lack of data there. To test whether propagating asthenospheric flow can be imaged beneath the Arabian Peninsula, we additionally performed a structural resolution test with three linear low-velocity anomalies as an input model (Figure 7). The resolution test results show that our model is able to resolve narrow (~150 km), linear low-velocity anomalies in the Arabian Peninsula south of Harrat Ithnayn.

Results
Depth slices through our S velocity mantle model are presented from 100 to 600 km depth in Figure 8. Remarkably, we found a rather sharp velocity contrast between the Arabian shield (low velocity) and the Arabian platform (high velocity) at all depths. High-velocity anomalies are also weakly observed beneath part of the Sinai Peninsula and east of the Main Ethiopian Rift at 100-300 km depth. Beneath Afar, Yemen, and the Main Ethiopian Rift, a large, strong low-velocity anomaly is observed from 100 to 600 km depth with a peak of around −270 m/s, which may indicate the presence of the upwelling Afar plume. The low-velocity anomaly beneath Afar seems to linearly extend toward the MMN line in the Arabian shield at 100-300 km depth, as previously imaged Park et al., 2007Park et al., , 2008Tang et al., 2018Tang et al., , 2019Yao et al., 2017), but with higher resolution (~150 km width). The amplitude of this linear low-velocity anomaly becomes weaker with depth, implying that it exists only at asthenospheric depths of 100-300 km. Strong, linear low-velocity anomalies are found along the Main Ethiopian Rift at 100-300 km depth with a peak of around −200 m/s, while patches of low-velocity anomalies exist below 300 km depth.
Vertical cross sections through our S velocity mantle model are shown in Figure 9. In Section A-a across Harrat Lunayyir, Harrat Ithnayn, and Harrat Hutaymah, a sharp velocity contrast between the Arabian shield and the Arabian platform is observed. This velocity difference may reflect a contrast between hot mantle upwelling from the asthenosphere beneath the Arabian shield responsible for harrat formation and cold, thick lithosphere beneath the Arabian platform (Al-Amri, 2015; Hansen et al., 2007; Tkalčić  , 2006). The depth extent of the lithosphere beneath the Arabian platform is estimated to be shallower than 200 km by previous receiver function studies (Al-Amri, 2015;Hansen et al., 2007;Tkalčić et al., 2006) and 200-250 km by a body-wave tomography study (Koulakov et al., 2016). Therefore, high-velocity anomalies deeper than 300 km beneath the Arabian platform are likely caused by strong vertical smearing due to poor resolution for the region (Figures 5 and 6).
In Section A-a, two low-velocity anomalies with peaks of around −180 m/s are located beneath the MMN line and Harrat Lunayyir with a depth extent of~300 km, likely connected to each other at 200-300 km depth. To refute the possibility of this connection being an artifact of smearing effects, we performed a structural resolution test with two low-velocity anomalies beneath Harrat Lunayyir and the MMN line (Harrat Ithnayn) separated by a small distance ( Figure S3). The result shows that the two low-velocity anomalies are well resolved and separated, with weak vertical smearing. Therefore, the observed connection between the two low-velocity anomalies beneath Section A-a seems to be robustly resolved and not caused by smearing effects. In Section B-b along the MMN volcanic line, a wide and strong low-velocity anomaly is found extending down to the mantle transition zone beneath the region between Afar and Yemen. We also observe an elongated low-velocity anomaly with a peak of~−180 m/s confined to depths of 100-300 km below the MMN line (hereinafter referred to as "the MMN line anomaly"), which is connected to the large low-velocity anomaly beneath the region between Afar and Yemen. We performed an additional resolution test along this transect ( Figure S4) and found that the anomalies are fairly well resolved, albeit with some smearing.
In Section C-c along the Main Ethiopian Rift, strong low-velocity anomalies with a peak of~−270 m/s are observed beneath Afar down to the mantle transition zone, which are linked to low-velocity anomalies with a peak of~−200 m/s at 100-300 km depth along the Main Ethiopian Rift, as observed beneath the Arabian shield in Section B-b. We investigated the robustness of the low-velocity anomaly in this transect with a resolution test ( Figure S5) and found that it is better resolved with higher amplitudes than the case for Section B-b ( Figure S4), which may result in stronger amplitudes of the low-velocity anomaly in Section C-c in our S velocity model (Figure 9).
In Section D-d across Harrat Lunayyir and Harrat Uwayrid, two low-velocity anomalies are located beneath Harrat Lunayyir and Harrat Uwayrid, but the connection between the two anomalies is weak even though this transect shows the best resolution among the four transects in the resolution tests ( Figure 5). Therefore, the low-velocity anomalies beneath Harrat Lunayyir and Harrat Uwayrid may be either weakly connected or not connected at all. The low-velocity anomaly beneath Harrat Lunayyir (~−180 m/s) is much stronger than the anomaly beneath Harrat Uwayrid (~−100 m/s).

Spatial Extent of Afar Plume Material Propagation Beneath Saudi Arabia
Various studies have been conducted to explain the source of the post-12 Ma volcanism in the Arabian shield. Previously, some geochemical studies have argued that the post-12 Ma volcanism originated from decompressional partial melting triggered by lithospheric thinning due to extensional stress related to the Red Sea rift (e.g., Bertrand et al., 2003;Moufti et al., 2012Moufti et al., , 2013. However, geophysical studies have found a low-velocity anomaly along the MMN line (e.g., Park et al., 2007Park et al., , 2008Tang et al., 2018Tang et al., , 2019Yao et al., 2017), and shear-wave splitting results show fast-axis directions aligned
Our results (Figures 8 and 9), coupled with shear-wave splitting studies, are consistent with the idea that plume material may be flowing in a northward direction from Afar beneath the Arabian shield. We image the MMN line anomaly much more clearly than in previous studies Park et al., 2007Park et al., , 2008Tang et al., 2018Tang et al., , 2019Yao et al., 2017), revealing a narrower lateral width (~150 km). In Cross Section B-b (Figure 9), the MMN line anomaly, confined within 100-300 km depth, is well observed from Afar via the MMN line to Harrat Ithnayn. The volume of the asthenospheric flow seems to be larger, moving upward to a shallower depth of~50 km just beneath Harrat Khaybar and Harrat Rahat. In Figure 10, we show variations of lithospheric depths in the Arabian Peninsula using LITHO1.0, which is a global model of the crust and upper mantle constructed using surface wave dispersion data with a period range of 25 to 200 s (Pasyanos et al., 2014). Lithospheric depths in LITHO1.0 are consistent with punctuated estimates of lithospheric depths using receiver function analyses (Hansen et al., 2007). Lithospheric thickness is observed to be thinner along the MMN line than in directions along the Red Sea or toward the Arabian platform from Afar, which may have guided northward flow. Along the MMN line, abrupt lithospheric thinning from 75 to 50 km thickness is found beneath Harrat Khaybar and Harrat Rahat (Figure 10), which may have caused decompression melting by upwelling of mantle material and consequently the observed large low-velocity anomaly and surface volcanism. A correlation between thin

10.1029/2020JB019668
Journal of Geophysical Research: Solid Earth lithosphere and the large low-velocity anomaly beneath Harrat Khaybar and Harrat Rahat would therefore seem natural. However, this small variation of lithospheric thickness is unable to be conclusively resolved in our model due to vertical smearing, given that a resolution test with a uniform-thickness low-velocity input anomaly is recovered with some thickness variation ( Figure S4).
Our model shows that the MMN line anomaly beneath the Arabia shield (Section B-b in Figure 9) does not extend below 300 km depth and that there is no low-velocity anomaly in the upper mantle beneath the Arabian platform (Figure 8). This finding refutes the hypotheses of local plumes beneath Arabia (Camp & Roobol, 1992;Koulakov et al., 2016) but supports a recent study that found a thin mantle transition zone only beneath the Ethiopia-Yemen traps (Kaviani et al., 2018), which means that hot upwelling from the lower mantle occurs only beneath the Afar region. In our model, it seems apparent that the low-velocity anomaly beneath Afar is rising at least from the mantle transition zone (Sections B-b and C-c in Figure 9), and the asthenospheric flow within 100-300 km depth is likely to start from this upwelling. We cannot resolve the existence, or absence, of the Jordan plume  due to limited resolution in northern Arabia and Syria.
Given the evidence for northward flowing plume material from Afar, we now try to answer the question of the extent of this flow. At 200 km depth (Figure 8b), the MMN line anomaly seems to be continuous from Afar to Harrat Ithnayn with a connection to Harrat Lunayyir in the west. Resolution test results with linear low-velocity input anomalies (Figure 7) confirm that our model can resolve asthenospheric flow from Afar to Harrat Ithnayn, but whether there is any further northern asthenospheric flow remains unresolved. Instead, plume material may have migrated westward from Harrat Ithnayn along a thin lithospheric channel ( Figure 10) to Harrat Lunayyir. The low-velocity anomaly beneath Harrat Lunayyir may not or may be weakly connected to the lowvelocity anomaly beneath Harrat Uwayrid (Section D-d in Figure 9), but definitely not further north along the Red Sea based on our resolution test (Figure 7). Since a geochemical study reports that no plume signal is found in Harrat Uwayrid (Altherr et al., 2019), the low-velocity anomaly beneath Harrat Uwayrid may have been caused by lithospheric thinning due to the Red Sea rift rather than asthenospheric flow of plume material. We also found a weak low-velocity anomaly (~−70 m/s) beneath Harrat Hutaymah, which seems to be weakly linked to the low-velocity anomaly beneath the MMN line (Figure 8 and Section A-a in Figure 9). We have shown that our model could recover the link between the MMN line and Harrat Hutaymah if it was real ( Figure 7). However, there is a weak low-velocity link between the MMN line and Harrat Hutaymah unlike the strong low-velocity anomaly between the MMN line and Harrat Lunayyir, which may indicate no plume flow to the east from the MMN line. This interpretation is consistent with a geochemical study that did not find any plume signal in Harrat Hutaymah (Konrad et al., 2016).
Geochemical studies that did not find a plume signal interpreted those results as absence of mantle flow across the entire Arabian Peninsula. Our model suggests that a lack of plume signal in certain harrats is more likely due to the limited extent of mantle flow. The MMN line anomaly in our model has a narrow width (~150 km), and the actual width is likely even narrower if we consider smearing effects and limited tomographic resolution. The size and spatial extent of the MMN line anomaly, as well as the timing of volcanism, appears to be controlled by lithospheric thickness (Figure 10). Hot asthenospheric material can easily rise beneath thin lithosphere, causing decompression melting and subsequent volcanism on the surface. In contrast, surface volcanism would be delayed under thick lithosphere, as it would take a longer

10.1029/2020JB019668
Journal of Geophysical Research: Solid Earth period of time for hot mantle upwelling to erode and then penetrate the lithosphere. Variations in lithospheric thickness may act as spatial controls on the northward propagation of plume material, which would have precluded a simple relationship between mantle flow and the timing of volcanism in western Arabia (Bertrand et al., 2003;Hopp et al., 2004).

Implications for the Harrat Lunayyir Seismic Swarm
Harrat Lunayyir (Figure 1a) is one of the post-12 Ma volcanic fields in the Arabian shield, and its last volcanic eruption occurred about 1,000 years ago (Camp et al., 1987). From April to June in 2009, an intense seismic swarm of over 30,000 earthquakes occurred in Harrat Lunayyir following episodic swarms since 2007, with the largest event reaching M W 5.4. However, this area was only recently considered a seismic zone in hazard maps for the region . After the seismic swarm, a surface rupture propagated 10 km across the northern area of Harrat Lunayyir (Pallister et al., 2010). Because of this unusual activity, many studies have been carried out to determine the origin of the swarm and the possibility of future seismic activity (Baer & Hamiel, 2010;Duncan & Al-Amri, 2013;Pallister et al., 2010).
The presence of several open fissures in Harrat Lunayyir and the absence of a major shock at the beginning of the seismic swarm indicate that it may have been caused by a magmatic source (Baer & Hamiel, 2010;. Strong thermal activity remains, as evidenced from high temperature groundwater and the presence of fumaroles in Harrat Lunayyir (Duncan & Al-Amri, 2013). The combination of low and high frequencies in the seismic waves has been interpreted as intrusion of magma and fracturing of the crystalline basement rocks (Pallister et al., 2010). A 3-D attenuation model using coda-wave attenuation data distinguishes a low-attenuation zone corresponding to the rigid basaltic field beneath Harrat Lunayyir but shows high attenuation below 6 km depth, which is interpreted as the top portion of a conduit that brings fluids and/or melts from deeper magmatic sources (Koulakov et al., 2014). Koulakov et al. (2015) estimated Vp/Vs ratios by creating tomographic models of P and S velocities and found two separate low-velocity anomalies with a high Vp/Vs ratio: one at 7-15 km depth likely corresponding to a steady-state magma reservoir and the other below 15 km depth interpreted as a conduit through which overheated material has been injected from deeper sources to activate the magma reservoir. Saibi et al. (2019) found a low-density zone in the middle-to-upper crust by 3-D inversion of the residual gravity field, which may be a magmatic source area that intruded into the upper crust beneath Harrat Lunayyir.
The earthquake swarm in Harrat Lunayyir has been associated with the Red Sea rift-related extension (Pallister et al., 2010;Sanfilippo et al., 2018). Pallister et al. (2010) proposed that crustal stress is primarily controlled by asthenospheric flow away from the Red Sea rift axis, comparing InSAR modeled dyke directions with the rift axis of the Red Sea. Sanfilippo et al. (2018) analyzed incompatible trace elements and proposed that the post-12 Ma volcanism in Harrat Lunayyir was caused primarily by decompression melting from lithospheric thinning due to the Red Sea rifting. However, it is unusual to observe active extension by means of dyke intrusion along a passive margin nearly 200 km away from the oceanic spreading center. Alternatively, the dyke intrusion could be tied to the presence of hot mantle rocks that are upwelling further east  as Ebinger and Belachew (2010) suggested.
Our tomographic results show that the low-velocity anomaly beneath Harrat Lunayyir is as strong as the anomaly beneath the MMN line (Figures 8 and 9). The anomaly beneath Harrat Lunayyir reaches down to~300 km depth and is likely connected with the low-velocity anomaly beneath the MMN line at 200-300 km depth (Section A-a in Figure 9), which is robustly resolved (Figure S3), implying an influx of plume material to basalts in Harrat Lunayyir. This interpretation is supported by previous geochemical studies, which found similarities between basalts from the Harrat Lunayyir and the MMN line. For example, Murcia et al. (2013) reported that both Harrat Lunayyir and northern Harrat Rahat in the MMN line have the same olivine composition. Based on trace element analyses, Duncan et al. (2016) found that Harrat Lunayyir and Harrat Rahat were formed from a plume-influenced source. On the other hand,  found magnetic lineations in Harrat Lunayyir aligned NE-SW using aeromagnetic data, which may partly correspond to dykes caused by asthenospheric flow from Harrat Ithnayn. Yao et al. (2017) found high-velocity anomalies at a depth range of deeper than 60 km beneath the Red Sea, which implies that the Red Sea is not an active spreading center but is instead passively rifting at shallow depth due to slab pull from the Zagros subduction zone. Therefore, it is unlikely that crustal stress beneath Harrat Lunayyir is controlled by mantle flow away from the Red Sea rift as Pallister et al. (2010) proposed.

Journal of Geophysical Research: Solid Earth
Based on these observations, we present a new geodynamic model beneath western Arabia. We posit that northward asthenospheric flow of plume material from Afar has reached Harrat Ithnayn and then extended to Harrat Lunayyir west of the MMN line guided along thin lithosphere as a natural channel (Figures 10 and  11). However, this asthenospheric flow likely does not extend to Harrat Uwayrid and Harrat Hutaymah, which is consistent with geochemical studies that found no plume signal there (Altherr et al., 2019;Konrad et al., 2016). Due to thin lithosphere, plume material may have risen and caused decompression melting, leading to the formation of Harrat Lunayyir. Furthermore, plume material may have contributed to triggering the seismic swarm in Harrat Lunayyir by sudden injection of unstable liquid volatiles through a conduit to a magma reservoir in crust (Sychev et al., 2017). Lithospheric thinning due to the Red Sea rifting-related extension may also be partly responsible for the seismic swarm.

Conclusions
We estimate upper mantle S velocity structure beneath the Arabian Peninsula and northeastern Africa between 100 and 600 km depth by inverting S and SKS relative traveltime residuals using the MCCC method. By incorporating SGS data, we have obtained denser data coverage beneath the Arabian shield than previous studies, which has allowed us to better infer reasons for the formation of harrats in western Arabia. Our model clearly shows a narrow (~150 km), elongated low-velocity anomaly beneath the Arabian shield from Afar at 100-300 km depth, which can be interpreted as northward asthenospheric flow from Afar when coupled with shear-wave splitting studies. This northward flow reaches up to Harrat Ithnayn and then turns to the west to Harrat Lunayyir following thin lithosphere as a natural channel. Low-velocity anomalies beneath Harrat Uwayrid and Harrat Hutaymah do not seem to be convincingly connected with the asthenospheric flow along the MMN volcanic line, which may support geochemical studies that show no plume signal to be found in these harrats (Altherr et al., 2019;Konrad et al., 2016). Therefore, the timing of volcanism and plume signals in western Arabia may not be age progressive from Afar due to the narrow width and limited extent of the asthenospheric flow from Afar as well as variations in lithospheric thickness.
Harrat Lunayyir, one of the post-12 Ma volcanic fields in the Arabian shield, is located away from the MMN volcanic fields on the surface, but the low-velocity anomaly beneath Harrat Lunayyir is likely connected with the low-velocity anomaly beneath the MMN line. Therefore, northward flow from the Afar plume may have contributed to the seismic swarm in Harrat Lunayyir in 2007-2009 together with lithospheric thinning due to Red Sea rift-related extension.