Home Research, Reports & Papers

Crustal Structure And Seismic Anisotropy Of Rift Basins In Somaliland


This study examines the crustal structure and seismic anisotropy of rift basins in Somaliland using teleseismic receiver function analysis, Bayesian inversion, common conversion point imaging, and 2D forward gravity modeling.

The results show 36.8-38.2 km of crust in southern Somaliland, while the central and northern regions have thinned crust with 5-6 km thick sediments. The study also reveals fast axis directions trending to 50-56° in the upper mantle, indicating azimuthal anisotropy is oriented in the regional Africa-Arabia plate motion.

The fast polarization of the crust in central Somaliland is oriented at -15°, indicating fossil deformation related to the NW-SE trending Late Jurassic rift event. Fast polarization for stations near the Gulf of Aden is oriented at 75-80°, suggesting crustal deformation associated with the Oligocene rift event. The crustal anisotropy at southern Somaliland is oriented at -85°, indicating a preserved far-field response of the WNW-ESE trending Late Cretaceous rift event.

By Mohammed Y. Ali, Mohammad Ismaiel & Ayoub Kaviani

Crustal Structure And Seismic Anisotropy Of Rift Basins In Somaliland

Mohammed Y. Ali,  Department of Earth Sciences, Khalifa University of Science and Technology, Abu Dhabi, UAE

Mohammad Ismaiel, Department of Earth Sciences, Khalifa University of Science and Technology, Abu Dhabi, UAE

Ibrahim M. Yusuf, School of Geography, Earth and Environmental Sciences, University of Birmingham, Birmingham, UK &

Ayoub Kaviani, Institute of Geosciences, Goethe University Frankfurt, Frankfurt, Germany




Rift margins provide insights into the processes governing the rupture of the continental lithosphere and the subsequence formation of sedimentary basins. The Proterozoic basement underlying Somaliland has been affected by multiple rifting; however, the crustal structure of these rifted basins remains unknown. This study utilized teleseismic receiver function analysis, Bayesian inversion, common conversion point imaging, and 2D forward gravity modeling to examine the crust and upper mantle of Somaliland. The results indicate 36.8–38.2 km of crust in southern Somaliland, while the central and northern regions feature thinned crust (~ 21 km) with 5–6 km thick sediments. The joint analysis of radial and transverse components of receiver functions and shear wave splitting revealed fast axis directions trending to 50–56° in the upper mantle, indicating that azimuthal anisotropy is oriented in the regional Africa-Arabia plate motion. Such orientation may have resulted from lattice preferred orientation of olivine from the asthenospheric flow. Additionally, the fast polarization of the crust in central Somaliland is oriented at − 15°, indicating fossil deformation in the thinned crust related to the NW–SE trending Late Jurassic rift event. Further, the fast polarization for stations near the Gulf of Aden is oriented at 75–80°, suggesting crustal deformation associated with the Oligocene rift event. The crustal anisotropy at southern Somaliland revealed fast polarization oriented at − 85°, indicating a preserved far-field response of the WNW-ESE trending Late Cretaceous rift event. Overall, the study provides for the first-time insight into the rift-related extensional strain fabric in the crust and upper mantle anisotropy induced by asthenospheric flow in Somaliland.



Rift initiation and formation of sedimentary basins are the key steps responsible for the deformation of the Earth’s surface. Seismic anisotropy is commonly used to examine the ongoing dynamic processes and paleo-deformation history of the continental rifts e.g.,1,2,3. Many studies on continental rift systems, such as the East African Rift, show a fast azimuthal anisotropy axis parallel to the rift axis with a splitting time as large as 2.5 s2,4,5,6. The anisotropy within the lithosphere in rift zones is produced by intruded melt and spreading-parallel fabric2,6,7,8. Asthenospheric flow beneath the rift shoulders produces lattice-preferred orientation (LPO) of olivine crystals where the fast axes of the mineral aggregates line up in the divergence direction7,9.

Somaliland is part of the Arabian-Nubian Shield, which is comprised of Proterozoic terranes that amalgamated during the Neoproterozoic10. Somaliland has been impacted by multiple rift events, which have led to the formation of sedimentary basins, thus providing a unique opportunity to study failed rift systems. The first event occurred in the Late Jurassic period, resulting in NW–SE trending grabens, which are attributed to be a result of the late-stage breakup of Gondwana and the detachment of Madagascar from Africa11,12,13,14. This rift led to the development of the NW–SE trending basins in Somaliland, such as the Guban, Dood Arale, and Daroor basins, as well as the Balhaf and Masilah basins in Yemen12,13,15,16. The second rift event occurred in the Late Cretaceous, and it was likely caused by a transcurrent motion that resulted in the separation of India and Madagascar12,13,17. This event led to the formation of WNW-ESE trending basins, such as the Nogal basin, which are primarily filled with 6–7 km of Upper Cretaceous sediments15. The final event took place in the Oligocene and is associated with the separation of Arabia from Africa and the subsequent opening of the Gulf of Aden12,13,18. This rift occurred along the Gulf of Aden and reactivated the Guban, Nogal, and Daroor basins. These rift events played a significant role in shaping the crustal architecture of the basins and caused thinning of the crust and upper mantle by up to a factor of two in some basins12,13,15.

The tectonic evolution of the sedimentary basins in Somaliland is well explained by a multi-rift extensional model12,13,15. As a result, it is anticipated that the crustal structure will show the signature of aborted rift systems and the crust will exhibit a ‘frozen-in’ deformation, which will provide useful insights into the rupture of the continental lithosphere. In this study, we utilized teleseismic records from five seismic stations, along with gravity data acquired across Somaliland, to determine the crustal and upper mantle structure and deformation inherited from the rifts in the region. Receiver function analysis provides an estimate of the Moho depth and VP/VS ratio. Bayesian inversion of receiver functions offers additional information on the VS structure beneath the stations. The common conversion point (CCP) imaging and 2D forward gravity modeling along a north–south oriented transect are used to better constrain the crustal architecture in Somaliland. Furthermore, splitting analysis of receiver functions and core-refracted shear waves is employed to explore azimuthal anisotropy within the crust and upper mantle. The findings of this study demonstrate that anisotropy within the crust is controlled by structural inheritance developed during the rift phases, while anisotropy in the upper mantle is generated by large-scale mantle flow. This is the first time, we report the detailed crustal structure, fossil crustal fabrics, and asthenospheric flow beneath the rift basins of Somaliland.


Geophysical data

From December 2019 to August 2022, we collected passive seismic data using broadband seismometers at five stations along a north–south transect in Somaliland (Fig. 1). At each station, a Nanometrics Trillium 120 three-component seismometer with a flat response between 120 s and 150 Hz was deployed at a depth of 1 m. The gravity data were acquired in two separate surveys. The northern part of the gravity data (from Eil Daraad to DS-2 well) was carried out in 201011. The remaining part of the data was acquired in July–August 2021 (Fig. 1). The data were collected using a Scintrex CG5 gravimeter with an instrument reading resolution of ± 0.001 mGal. The distance between the gravity stations in the northern and southern parts was approximately 1 km and 5 km, respectively, but some modification in the sampling distance was needed due to access issues in some places.

Figure 1

figure 1

(a) Geology map of northern Somalia (Somaliland), modified from13. The map shows the location of the five broadband seismometers (green triangles) deployed in Somaliland. SRTM15 + data is used to display the elevation and bathymetry19. The location of two drilled wells BD-1 (Biyo Dader-1) and DS-2 (Dagah Shabel-2) are shown. Inset is a location map of Somalia (the green box illustrates the study area) and the surrounding region within the tectonic framework of surrounding plates. The red arrow shows the absolute plate motion of major tectonic plates20. (b) Azimuthal equidistant projection map centered at the seismic station in Burao showing geographic distribution of teleseismic events (yellow circle). (c) and (d) Rose diagrams for distance windows of 30–98° and 90–130° showing the number of teleseismic events as a function of backazimuth. Between September 2020 and August 2020, a total of 354 events occurred in the azimuthal distance window of 30–98°, while 198 events occurred in the 90–130° distance window. The figure was generated using Generic Mapping Tools (GMT)21.

H-κ stacking

Receiver function analysis of passive seismic data reveals crustal thickness variation along the north–south transect in Somaliland (Fig. 2 and Table 1). The total number of earthquakes with a moment magnitude (Mw) ≥ 5.5 and epicentral distance of 30–98° that occurred between September 2020 and August 2022 is 354 (Fig. 1b), although the seismic stations only recorded 150 to 260 events. After adopting a signal-to-noise ratio (SNR) 3 and a visual assessment of receiver functions, the backazimuth coverage for these teleseismic events is concentrated between 60° and 90°, and the number of selected radial (R) component receiver functions for H-κ stacking per station is 20–30. In the R-component, the coherent PS phase appears between 2 and 4 s after the P-wave arrival (Fig. 2). The multiple reflected waves, PPPS and PSPS + PPSS arrive between 8–13 s and 10–18 s, respectively, and are weaker than PS phase. The strong positive phases that arrive at approximately 2 s in Eil Daraad and Burao, 3 s in Hagal, and 4 s in Haydh Duato and Dharyeley are the Moho PS, inferring crustal–mantle discontinues (Moho) at depths of 21–38 km and VP/VS ratios of 1.6–1.76. It is interesting to note that the northern station of Eil Daraad, which is at the shoreline of the Gulf of Aden, reveals a Moho depth of 21.1 km, while Burao station, which is located 150 km inland from the coast, shows an exceptionally shallow Moho depth of 22.3 km with a relatively high VP/VS of 1.76. This contrasts with the southern stations of Haydh Duato and Dharyeley, which reveal Moho depths of 36.8 km and 38.2 km, respectively. The Hagal station shows a Moho depth of 32.0 km (Fig. 2d and Table 1).

Figure 2

figure 2

(a)–(e) Receiver functions for all five stations in Somaliland. The top panels display the receiver functions plotted with increasing backazimuth angle. The stacked receiver functions are shown in grey panels and calculated arrival times for PS conversion phases (PS, PPPS, PSPS + PPSS) for the optimal solution in H-κ stacking are plotted with green vertical lines. The bottom panels show the Moho depths (H) versus VP/VS ratios (κ) from the H-κ stacking technique. The optimum solution for H and κ for each station is shown in the bottom panels. Refer to Fig. 1 for the location of the stations.

Table 1 Moho depth and VP/VS (κ) for the five stations located in Somaliland.

Station Latitude (°) Longitude (°) Elevation (m) Moho depth (km) VP/VS
Eil Daraad 10.75869 45.60419 12 21.1 ± 1.1 1.61 ± 0.04
Hagal 10.25256 45.74044 382 32 ± 0.9 1.70 ± 0.03
Burao 9.51717 45.56362 1030 22.3 ± 0.5 1.76 ± 0.02
Haydh 8.71734 45.73534 799 36.8 ± 2.0 1.61 ± 0.04
Dharyeley 8.39776 45.83372 762 38.2 ± 1.7 1.60 ± 0.03

V S profiles

The VS profile at station Eil Daraad increases from 2.2 to 2.6 km/s in the top layer down to 6 km depth, and then, after a slight decline, it rises sharply from 2.2 to 3.5 km/s in the 6–12 km depth range (Fig. 3a). The VS steadily increases in the depth range of 12 to < 32 km with no obvious change in the Moho interface. However, the VS reaches up to 4.0 km/s at a depth of 32 km, which may indicate that the Moho boundary is relatively deeper compared to the H-κ stacking results (Fig. 2a). The Vs profile for Hagal station, in contrast to that of Eil Daraad station, does not reveal the significant thickness of the top low-velocity layer, but instead shows a gradual increase in velocity from 3.2 to 3.6 km/s in the top 12 km layer and from 3.6 to 4.2 km/s up to 31 km depth (Fig. 3b). The Moho depth, which is slightly shallower than the H-κ stacking result, is noted at the interface with a VS of 4.2 km/s. The presence of 1–2 km of sediments is consistent with drilling data of BD-1 and DS-2 wells, which penetrated the crystalline basement at depths of 1,471 and 1,444 m, respectively. At Burao station, the top 5–6 km layer has a low VS (2.2 km/s), and at 6 km depth, the velocity suddenly increases from 2.2 to 3.4 km/s (Fig. 3c). The VS increases from 3.2 to 4.0 km/s, after a slight decrease, up to a depth of 23 km, at which point it remains constant. The Moho depth is marked at the interface with a VS of 4.0 km/s and matches with the H-κ analysis. The low VS (2.2 km/s) is not seen in the Haydh Duato and Dharyeley stations. The VS profile at Haydh Duato station remains constant in the upper crust, whereas the lower crust, down to a depth of 36 km, shows a stepwise increase from 3.2 to 4.0 km/s (Fig. 3d). The Moho depth determined from the inversion results agrees well with receiver function analysis. The VS profile at Dharyeley station shows high and low VS layers within the upper crust. The VS is 3.6 km/s between 12 and 34 km depth and increases to 4.0 km/s at 38 km depth. Dharyeley, like Haydh Duato, appears to have a limited (1–2 km) sedimentary cover with low VS (2.2–2.5 km/s) (Fig. 3e).

Figure 3

figure 3

(a)–(e) VS structure beneath the five stations obtained from Bayesian inversion of P-wave receiver functions. The top panels show the observed (blue curves) and modeled receiver functions. The modeled curves are the ones (out of 21 chains) that converge successfully and retain acceptance rates in the process of multiple independent parameter search paths. The bottom panels show the VS posterior distributions and corresponding probabilities for interface depth. The red curves are the initial VS model and the white arrows represent the Moho depth from H-κ stacking. A low VS top layer is noticed under Eil Daraad and Burao stations. At Hagal, Haydh Duato, and Dharyeley stations, the VS values in the top layer match with the upper crustal velocity. The Moho depth from Bayesian inversion matches well with H-κ stacking results. Refer to Fig. 1 for the location of the stations.

Crustal structure: CCP and gravity modeling

The CCP imaging along the profile of seismic stations indicates a high-velocity gradient identified as the Moho boundary, which coincides well with the H-κ results (Fig. 4). The crustal thickness increases from 23 km in Eil Daraad to 30 km at Hagal and then reduces to 24 km at Burao (Fig. 4a). The Moho depth is interpreted at 30 km beneath Haydh Duato and Dharyeley stations. This Moho topography is shallower than that determined by the H-κ results (37–38 km). At a depth of 10 km, a strong middle-lower crustal interface is found beneath Burao, Haydh Duato, and Dharyeley stations. However, the wide gaps between the stations precluded the resolution of shallow intra-crustal discontinuities (white gap in Fig. 4a).

Figure 4

figure 4

(a) CCP imaging along Eil Darrad–Dharyeley transect. Pink inverse triangles show the location of the stations and black-filled circle symbols represent the Moho depths estimated from the H-κ stacking. CCP image reveals a major high-velocity gradient between 20 and 30 km depth (yellow dotted line), which is interpreted as the Moho boundary. A good match is noticed for Moho depths observed from H-κ stacking and CCP imaging beneath Eil Daraad, Hagal, and Burao stations, while shallow Moho depths are seen at Haydh Duato and Dharyeley. Due to wide gaps between stations (white gaps), the shallow intra-crustal discontinuities are unresolved. The scale bar represents the amplitude of positive (red) and negative (blue) polarities of arrivals. (b) Two-dimension density model along the transect-A running across Somaliland (see Fig. 1 for location). The Moho depths (black-filled circles) are derived from H-κ stacking and densities of sedimentary column constrained are obtained from drilled wells BD-1 and DS-2, which penetrated the Proterozoic basement at 1,471 m and 1,444 m, respectively. The crustal model suggests an asymmetric sedimentary basin in the Burao region with a maximum depth of 6 km. A sauce-shaped basin with a maximum depth of 2 km is noticed beneath Hagal. Beneath Eil Daraad and Burao, the Moho gets shallower (~ 22 km), and the basement is at a 5–6 km depth. The fit between observed and calculated gravity anomaly is particularly good.

Two-dimensional forward gravity modeling along the north–south transect delineates the geometry of the crust and rifted sedimentary basins (Fig. 4b). The Bouguer anomaly between Burao and Dharyeley remains relatively constant (− 160 mGal) for 160 km, whereas the gravity anomaly increases to − 30 mGal towards Eil Daraad. The crustal model reveals thick sedimentary basins with densities of 2300–2400 kg/m3 beneath Eil Daraad and Burao. The geometry of the basin in Burao seems asymmetrical, with maximum sediment thickness on the northern side. At Hagal station, a saucer-shaped basin with a maximum sedimentary thickness of 2 km is observed. The sediment thickness is thin (< 1 km) at Haydh Duato and Dharyeley stations. The crustal thickness is 36–38 km beneath Haydh Duacato and Dharyeley and it reduces to 24 km at Burao. The crustal thickness increases again to 32 km beneath Hagal and thins to 20 km towards the Gulf of Aden in Eil Daraad.

Seismic anisotropy

A joint analysis technique on the radial (R) and transverse (T) components of the receiver function dataset allowed for determining crustal anisotropy at the five seismic stations. 25–60 P-wave receiver functions were selected from among 150–260 events that were recorded based on SNR 3 to offer an accurate assessment of crustal anisotropy. The binned stacked R and T components arranged in ascending backazimuth for each station are shown in Fig. S1. Figure 5 shows three computed individual objective functions and a joint objective function for each station. The anisotropy calculated at Eil Daraad station using R cosine energy maximization, T energy minimization, and R cross-correlation maximization shows good agreement with a fast direction (ϕ) of 78–80°, (Fig. 5a). The R cosine energy and R cross-correlation techniques estimate a splitting time of (δt) 0.4 s, while the T energy method yields a value of 1.3 s. The joint analysis shows ϕ of 80° with a δt of 0.4 s. At Hagal station, R cosine energy, R cross-correlation, and T energy minimization show good agreement with ϕ of 72–83°. (Fig. 5b). The joint analysis yields ϕ of 75° and δt of 0.4 s. In the case of the Burao station, R cosine energy, R cross-correlation, and T energy methods show a good match with ϕ of − 22 to − 15° and δt of 0.4 s (Fig. 5c). The joint objective function reveals a ϕ of − 15° and a δt of 0.4 s. For the Haydh Duato and Dharyeley stations, the joint method shows ϕ of 95° and a δt of 0.3–0.4 s (Fig. 5d,e). However, at Dharyeley station, the T energy minimization deduced a δt of 1.3 s. The difference in δt deduced from the T energy comparison to R cosine energy and R cross-correlation methods is due to the noise level associated with recorded in these components or the effect of dipping anisotropic layer.

Figure 5

figure 5

(a)–(e) Seismic anisotropy results were calculated from three objective functions (R energy maximization with cosine moveout correction, Radial cross-correlation maximization, and T energy minimization) and a joint objective function on R- and T- components of receiver functions. The objective functions are plotted in a 2-D plane of ϕ and δt in the range of 0–360° and 0–1.5 s, respectively, with an interval of 1° and 0.02 s, respectively. The ϕ and δt are marked by a white cross symbol. Color bars indicate variation in objective functions.

Splitting analysis is also performed on core-refracted shear phases (SKS, SKKS, and PKS, collectively called XKS henceforth) to determine the upper mantle anisotropy from teleseismic data recorded at the five stations. To compute XKS splitting parameters, rotation-correlation, and T-component energy minimization approaches were applied on a selected time window to XKS waveforms from teleseismic events that occurred at epicentral distances of 90–130° and had Mw > 5.7. Around 150–160 seismographs were successfully associated with earthquake events (a total of 198) at individuals linked to seismic stations. The elliptical particle motion is best linearized by the minimum energy method over rotation-correlation (Fig. S2S4). Therefore, splitting measurements were obtained from the minimum energy method (Table 2). To examine the effect of crustal anisotropy on the XKS splitting measurements, good and fair non-null measurements were fitted with synthetic curves generated by a two-layer anisotropic model (Fig. 6). The splitting parameters (ϕ, δt) from crustal anisotropy study were used to constrain the upper layer and the ϕ and δt for the lower layer (i.e., upper mantle) were varied to best fit with measurements. The non-null measurements at Eil Daraad and Hagal stations fit well for an upper layer with a ϕ of 70–80° and a δt of 0.4s, and a lower layer with a ϕ of 50°–55°, and δt of 1.1–1.3 s. For the Burao station, the data fit well to an upper layer with a ϕ of − 28°and a δt of 0.4 s, and a lower layer with a ϕ of 50° and a δt of 1.0 s. At Haydh Duato and Dharyeley stations, the predicted values of ϕ and δt for the upper layer are − 85° and 0.4 s, respectively. The splitting parameters (ϕ and δt) of the lower layer are 54–56° and 0.7–1.1 s. The best-fitting model of SKS splitting data for the five stations has a strong lower anisotropic layer with an average ϕ of 50–56°.

Table 2 Crustal anisotropy from joint analysis of R and T- components of receiver function and SKS splitting results for the five teleseismic stations in Somaliland.

Station Crustal anisotropy SKS splitting
Crust Upper layer Lower layer
ϕ (°)


ϕ (°)


ϕ (°)


Eil Daraad 80 ± 8 0.4 ± 0.10 80 0.4 50 1.3
Hagal 75 ± 10 0.4 ± 0.20 70 0.4 55 1.1
Burao − 15 ± 6 0.4 ± 0.05 − 28 0.4 50 1.0
Haydh − 85 ± 8 0.3 ± 0.20 − 85 0.4 54 1.1
Dharyeley − 85 ± 10 0.4 ± 0.10 − 85 0.4 56 0.7


Figure 6

figure 6

(a)–(e) SKS splitting analysis results for five stations. Two-layer modeling is performed to fit the observed measurements from SKS splitting. The upper layer anisotropy varies only around the direction of ϕ estimated from crustal anisotropy. The top panels show the ϕ versus backazimuth and the bottom panels display δt versus backazimuth plots. Good and fair splitting measurements are shown with red plus and blue cross symbols, null measurements are shown with red circles, while theoretical curves for the best fit are shown in green color curves. The parameters of the upper and lower layers from the minimum energy method are placed below the top panel. The direction of ϕ ranges from 50 to 56° and delay time is from 0.7 to 1.3 s for the lower layer.


Moho topography and V P/V S ratio

This is the first study to use receiver function analysis to show lateral variation in the Moho depth across Somaliland (Fig. 2 and Table 1). A thick unstretched crust of 36.8–38.2 km is observed beneath Dharyeley and Haydh Duato in southern Somaliland. The Moho boundary reaches a shallow depth of 22.3 km under Burao, which indicates a decrease towards the north. The Moho depth again deepens and reaches 32.0 km at Hagal and then reduces to 21.1 km at Eil Daraad. The Moho topography variation suggests that the continental crust is stretched more beneath Burao compared to the northern (Hagal) and southern (Haydh Duato and Dheryeley) stations. Moreover, continental extension took place from Hagal (32.0 km) to Eil Daraad (21.1 km), towards the Gulf of Aden. The rapid change in Moho depth is due to the transition of the continental crust beneath the Somaliland margin to the oceanic crust in the Gulf of Aden12,22. The results of receiver functions suggest that the unaltered continental crust with a thickness of 37–38 km has been strongly thinned to ~ 21 km to form the rift basins with up to 5–6 km of sedimentary cover. These results agree well with process-oriented gravity and flexure modeling across the conjugate margins of Somaliland and Yemen12. Previous receiver function studies from eastern Afar in Ethiopia and Djibouti revealed a crustal thickness of 40 km beneath the African shield e.g.,23,24, and Moho depths of 15–25 km along the rift e.g.,25,26. Moreover, the inferred crustal thickness beneath the Yemen Plateau is ~ 35 km and it reduces to 22 km near the coast of the Gulf of Aden and 5–10 km within the Gulf of Aden27. Additionally, receiver functions and wide-angle seismic data analysis from onshore Yemen, Oman, and central Socotra Island reveal an average crustal thickness of 28 –35 km, and crustal thinning towards the Gulf of Aden28,29,30.

The VP/VS ratios range from 1.6 to 1.76 at the five stations (Fig. 2 and Table 1). The average VP/VS ratios of the crust beneath the western Yemen Plateau, the northern margin of the eastern Gulf of Aden, and Socotra Island are between 1.67 and 1.9127,28,29. The crustal composition in onshore Yemen is felsic with VP/VS ratios of 1.67–1.77, while the VP/VS ratio increases up to 1.92 close to the Red Sea coast27. The VP/VS ratio (1.76) at the Burao location is close to the global average for the continental crust (VP/VS = 1.73;31). The high VP/VS ratio indicates a reduction in the VS and can be indirect evidence for the presence of fluids or melts within the rift system due to asthenospheric upwelling26. The thin crust with a slightly high VP/VS ratio observed at station Burao implies an enrichment of mafic material. A likely explanation is magmatic underplating/mafic lower crust or the presence of melts generated during the process of rifting. The other possible reason for VP/VS = 1.76, is the sediment infill in the rifted graben. Furthermore, the inversion of receiver functions reveals a low VS (2.2–2.6 km/s) at the top layer at Eil Daraad and Burao stations, while the low VS layer is absent beneath Hagah, Haydh Duato, and Dharyeley stations (Fig. 3). Moreover, a sharp change in VS at the crust–mantle boundary is not observed under Eil Daraad. However, other stations show strong agreement between the Moho depth from the H-κ stacking and the inversion of receiver functions. As a result, crustal deformation from the rift phase may have affected the base of the crust and therefore does not result in a notable change in the VS profile at Eil Daraad.

Crustal structure across Somaliland

CCP imaging is performed across all the stations to image the Moho and intra-crustal discontinuities (Fig. 4a). The Moho boundary is observed to match very well with the H-κ stacking results for Eil Daraad, Hagal, and Burao stations. Additionally, a notable intra-crustal discontinuity is observed at a depth of 10 km. The upward trend in the Moho depth is noticed towards the Gulf of Aden and thinned continental crust at Burao. The crustal architecture across Somaliland is also defined using 2D forward gravity modeling (Fig. 4b). The model outlines a 110 km wide sedimentary basin with an average density of 2400 kg/m3 and a maximum depth of 6 km in the Burao region. A smaller basin with a density of 2300 kg/m3, a width of 40 km, and a maximum sediment thickness of around 1 km, occurs north of Hagal. The Eil Daraad region also has a rift basin with 5–6 km thick sediments that have a density of 2400 kg/m3. These findings are consistent with the 3D inversion of gravity and magnetic data in the Dood Arale basin, east of Somaliland, which indicates an NW–SE trending basin with an average density of 2400 kg/m3,17. The Dharyeley and Hagal regions appear to have an unaltered continental crust, which has a thickness of 36–38 km, while the Burao and Eil Daraad regions have a thinned continental crust (20–24 km). The stretched continental crust in the interior of Somaliland is reported in the Guban, Nogal, and Daroor basins11,13,15. In the Nogal basin, backstripping of biostratigraphic data and process-oriented gravity modeling suggest that the crust thinned down from 35 km at the edge of the basin to 20 km at the basin axis due to the Late Jurassic, Late Cretaceous, and Oligocene rift events15. Furthermore, the crustal structure derived from backstripping studies of biostratigraphic data and process-oriented gravity and flexure modeling suggests Moho depths of about 30 km at BD-1 and DS-2 wells and thinning of the crust (~ 20 km) along the Somaliland margin of the Gulf of Aden12. The thinned continental crust confirms the extension of the continental lithosphere and the rift events that occurred in the region.

Seismic anisotropy beneath the rift systems

Evidence is presented here for the first time for crustal anisotropy produced by rift-related extensional strain and upper mantle anisotropy induced by asthenospheric flow in Somaliland. Table 2 summarizes the results of PS splitting analysis for the crustal anisotropy and SKS splitting analysis for the upper mantle anisotropy. Figure 7 shows a schematic illustration summarizing the azimuthal anisotropy direction observed in the crust and upper mantle beneath the seismic stations. The stations close to the Gulf of Aden (Eil Daraad and Hagal) display a crustal anisotropy with a fast polarization direction aligned in the rift axis of the Gulf of Aden (Fig. 5a,b). This can be explained by the Oligocene rift event that caused the opening of the Gulf of Aden. The crustal anisotropy in Burao is parallel to the NW–SE trending Late Jurassic rift basins13, which suggests a thinned and deformed crust because of the Late Jurassic rift event. Under the unthinned continental crust beneath Haydh Duato and Dharyeley stations (thickness 36–38 km), the crustal anisotropy is observed approximately in the WNW-ESE direction, which is parallel to the Late Cretaceous rift event15. This may indicate a preserved far-field response of the Late Cretaceous rift event.

Figure 7

figure 7

Schematic diagram illustrating the anisotropy observed in the crust and upper mantle beneath the seismic stations. Crustal thickness is from H-κ stacking results. The figure was generated using Generic Mapping Tools (GMT)21.

The SKS splitting parameters were inverted for two-layer models, with the upper layer constrained by the parameters of the crustal anisotropy obtained from the PS splitting analysis of receiver function data (Fig. 6). For the lower layer (i.e., upper mantle), a relatively homogenous pattern in the fast anisotropy directions of 50–56° (NE-SW orientation) is observed with delay times of 0.7–1.3 s. The estimated fast direction of the crust anisotropy in the rift basins of Somaliland is consistent with the observations from other rift systems, such as the main Ethiopian rift, where the splitting direction fast polarization direction is parallel to the rift axis4,6,25,32. Laboratory experiments and xenolith studies indicate upper mantle anisotropy is produced by the LPO of olivine crystal resulting from shear strain8,33,34. In general, the LPO of olivine crystal develops in response to upper mantle flow or frozen-in lithospheric fabric from past tectonic events35,36. The upper mantle anisotropy of all five stations aligns with the regionally occurring NE-oriented direction of the absolute plate motion (APM) of the African and Arabian Plates20,37. We believe that the large-scale mantle flow induced by the drag at the base of the lithosphere due to the plate motion plays the main role in producing the upper mantle anisotropy beneath Somaliland. On the other hand, the upper mantle anisotropy directions are also subparallel to the opening direction of the East African Rift, such as the Malawi rift zone38. Therefore, it is plausible that the anisotropy of the upper mantle could be due to the superposition of two sub-parallel layers, where the upper layer originates from finite extensional deformation in the lithosphere, and the lower layer is associated with the large-scale mantle flow induced by the APM. The inferences made here suggest that the continental lithosphere of Somaliland has evolved through substantial deformation in the past. A weak lithosphere with elastic thickness, Te, of 5 km is observed in the Nogal basin and Somaliland–Yemen conjugate margins12,15. High heat flow at the sub-lithospheric level and substantial internal deformation from multiple rift events may be responsible for the heterogeneous nature of the lithosphere and rheological weakening.

Previous anisotropy studies from East Africa also revealed NE-SW fast polarization directions linked with the LPO induced by plate motion e.g.,5,6,39,40. The upper mantle anisotropy of the northern conjugate margin of the Gulf of Aden exhibits NW–SE, NE-SW, and mixed orientations, which were linked to fossil lithospheric mantle anisotropy41. However, the upper mantle anisotropy beneath the Arabian Peninsula provided evidence for plate-driven viscous drag in the asthenosphere42.


This study employed a comprehensive array of geophysical techniques, encompassing teleseismic receiver function analysis, Bayesian inversion, common conversion point imaging, and 2D forward gravity modeling, to thoroughly investigate the crustal structure and azimuthal anisotropy along an N-S transect within Somaliland. Our findings reveal a 36.8–38.2 km crust in southern Somaliland, while the central and northern regions exhibit a thinned crust of approximately 21 km overlain by 5–6 km sediments. Joint analysis of receiver functions and shear wave splitting elucidates fast axis directions (50–56°) in the upper mantle, reflecting azimuthal anisotropy aligned with Africa-Arabia plate motion, possibly due to olivine lattice orientation influenced by asthenospheric flow. Furthermore, insights into crustal polarization directions have revealed compelling information. The observed -15° orientation in central Somaliland indicates fossilized deformation linked to the NW–SE trending Late Jurassic rift event. Similarly, seismic stations proximate to the Gulf of Aden exhibit fast polarization directions aligned at 75–80°, indicative of crustal deformation associated with the Oligocene rift event, which played a pivotal role in the inception of the Gulf of Aden’s formation. Moreover, the − 85° orientation of crustal anisotropy in southern Somaliland unveils a well-preserved far-field response to the WNW-ESE trending Late Cretaceous rift event. In essence, this study presents insights into rift-induced extensional strain within the crust and asthenospheric flow-induced upper mantle anisotropy. These findings enhance our understanding of rift-related geodynamics and the geological history of Somaliland.


Receiver function analysis

For receiver function analysis, the origin time and location of teleseismic events were downloaded from the United States Geological Survey global earthquake catalogue (Fig. 1b,c). The projected arrival time for the primary P-phase is calculated using the TauP module in the ObsPy package43. After matching the proper events with seismograms, the raw waveform is cut to 200 s around the projected arrival time. Then, the pre-processing of waveforms is performed by removing the linear trend and average from the raw seismograms before downsampling the data to 25 Hz. Subsequently, the waveform was filtered by applying a Butterworth bandpass filter with 0.05 Hz and 2.0 Hz to remove the high-frequency artefacts. The receiver function analysis is performed in the time-domain deconvolution44 and the grid-search algorithm45 is used to evaluate the Moho depth (H) and VP/VS (κ) ratio from P-wave receiver functions beneath the station. For deconvolution, a 30 s time window is selected and radial components with SNR > 3.0 are retained for H-κ stacking. The assumed crustal velocity is 6.5 km/s46 and weights to converted PS phases, PS, PPPS, and PSPS/PPSS are w1 = 0.6, w2 = 0.3, and w3 = 0.1, respectively. The Moho depth is varied from 20 to 50 km and κ value from 1.5 to 2.0 in the grid-search algorithm. Sensitivity analysis with a distinct set of weights (w1 = 0.5, w2 = 0.3 and w3 = 0.2; and w1 = 0.7, w2 = 0.2 and w3 = 0.1) was performed, and good convergence is observed for each set. The error estimation for the Moho depth and κ is determined from a bootstrap resampling method47. The principle of the bootstrap method to quantify the uncertainty is by generating a resampled group randomly, estimating the statistics for resample groups, and 1σ standard deviation is calculated from distribution analysis The receiver function analysis is performed using SplitRFLab toolbox48.

A Markov chain Monte Carlo transdimensional Bayesian inversion of P-wave receiver functions is adopted to construct a VS depth profile beneath each station49. The initial velocity model consists of the top sedimentary layer, upper, middle, and lower crust and upper mantle with VS range from 2 to 5 km/s for a vertical depth from the surface of 60 km (Fig. 3). The inversion was performed for a maximum number of layers 20, VP/VS ratios of 1.6–1.9 and 21 chains, where each chain performing 100,000 burn-in iterations and 50,000 main-iterations with 40–45% acceptance rate. 5% of the complete dataset was rejected as outlier chains, and the final posterior distribution of 100,000 models was assembled based on likelihood. For the inversion of P-wave receiver functions, the radial component for the time window − 5 to 35 s is obtained from the stacked receiver function (see Fig. 2) analysis. The present model does not deal with complex crustal structures such as dipping or anisotropic layers at the crustal level.

Crustal anisotropy and SKS splitting

Splitting analysis of the Moho converted P-to-S (PS) phase on receiver functions is used to investigate seismic anisotropy in the curst21,50,51. On the other hand, splitting analysis of XKS phases is used to probe mantle anisotropy52. The splitting analysis (of both Ps and XKS phases) assumes that a seismic shear wave splits into fast and slow components with orthogonal polarization directions by traveling through an anisotropic medium. Then, a delay time is accumulated between the two components that are a function of the strength of the anisotropy of the medium and the length of the ray path traveled in the anisotropic layer. In a shear-wave splitting analysis, the polarization direction (ϕ) of the fast component and the delay time (δt) between the two components are measured as proxies for the direction and strength of anisotropy in the underlying medium.

An integrated technique is used to determine the crustal anisotropy using radial and transverse components of receiver functions. We used the Seispy package to estimate the crustal anisotropy from three individual methods based on radial energy maximization, radial correlation coefficient maximization, and transverse energy minimization21, and a joint analysis performed on the average of three individual objective functions53. In the case of radial energy maximization, the estimation of ϕ and δt is based on the ratio of peak energy derived by stacking all the radial components of receiver functions after and before the time correction of the cosine moveout in the Ps arrival time. The mathematical formula for the ratio of peak energy is presented in53. The peak energy ratio is calculated by varying ϕ in the range of 0–360° with an interval of 1° and δt from 0.0 to 1.5 s with an interval of 0.02 s. The second technique is based on maximizing the radial correlation coefficient for an anisotropy-corrected R-component receiver function. To achieve the anisotropy-corrected components, radial and transverse receiver functions are projected in the assumed fast and slow polarization directions and time shifting (half of the splitting time, δt/2) is applied in the forward direction for the fast component and in the backward direction for the slow component. The analysis involves a grid search for ϕ and δt to evaluate the best fit for the crustal anisotropy. The third method, transverse energy minimization, exploits the basic concept of P– to S-wave conversion at the Moho boundary that reverberation is fully polarized in the radial direction, therefore, the presence of energy in transverse component points to anisotropic nature of the crust53,54. For assumed values of ϕ and δt, the anisotropy-removed transverse receiver functions are stacked and used to evaluate the ratio with the initial T energy. The above-mentioned methods are based on the different properties of the PS converted phase in response to inhomogeneous structure in the crust and sensitivity to noise present in receiver function data. Thus, a joint analysis is performed to obtain more consistent results of crustal anisotropy from receiver functions53. A bootstrap method is implemented to obtain the uncertainty in our estimates.

In XKS splitting analysis, the underlying anisotropic medium is probed using multiple methods such as rotation-correlation55, minimum energy, and Eigenvalue technique52. The basic principle of all these techniques is to remove the effect of splitting between the two polarized waves in the fast and slow directions by performing a grid search for the parameters ϕ and δt. We used the Matlab code SplitLab56 for data pre-processing and splitting measurements. This code also allows the examination and modeling of the azimuthal variation of the individual splitting parameters for depth-dependent anisotropy beneath seismic stations.

CCP imaging

CCP stacking technique was carried out along a transect between Eil Daraad and Dharyeley. Receiver functions with SNR > 3.0 were utilized to convert the reverberations phase (PS, PPS, PSS) to depth section following routine procedures: (a) using a 1-D IASP91 velocity model57, the receiver functions converted from time (s) to depth (km) and (b) estimation of piercing points of the rays and then projection and stacking amplitudes occurring within a rectangular bin58. The Ps arrival time and corresponding geographic location are calculated for each receiver function and a linear moveout correction is applied by multiplying the standard slowness factor (± 0.15 s per degree with an increment of 0.01 s per degree) and relative epicentral distance with reference to median epicentral distance. Subsequently, the receiver functions that share common conversion points were stacked based on their epicentral distances. The maximum depth for time-to-depth conversion is considered 800 km with a vertical spacing of 1 km, and a bin spacing interval of 5 km with a width of 100 km.

Gravity data processing and modeling

The gravity data were corrected for drift, tide, elevation, and latitude. The data were reduced to Bouguer anomaly using an average density of 2670 kg/m3. The data was tied to the base station (24.420139° N, 54.502822° E and g of 978,883.210 mGal) in Abu Dhabi, United Arab Emirates. The Geosoft® Oasis Montaj software was used to construct a density cross-section along the north–south transect. The model is based on a 2D flat Earth model for gravity computations of structural block and it is extended along the profile to avoid edge effects. We assumed an initial crustal model using results from the H-κ stacking, inversion of receiver functions, and CCP imaging. The depth to the basement, thickness, and densities of the sedimentary layers are constrained by BD-1 and DS-2 wells. The outcrop geology is utilized to constrain the density cross-section along the transect. The gravity response of the initial model is compared with the observed gravity anomaly and the crustal model was modified until the best fit was achieved between the observed and calculated gravity anomalies.

Data availability

The passive seismic and gravity data were deployed and acquired by the authors. The data are available from the corresponding author.


  1. Gao, S. et al. SKS splitting beneath continental rift zones. J. Geophys. Res. Solid Earth 102, 22781–22797. (1997).

  2. Kendall, J. M., Stuart, G. W., Ebinger, C. J., Bastow, I. D. & Keir, D. Magma-assisted rifting in Ethiopia. Nature 433, 146–148. (2005).

  3. Eilon, Z., Abers, G. A., Jin, G. & Gaherty, J. B. Anisotropy beneath a highly extended continental rift. Geochem. Geophys. Geosyst. 15, 545–564. (2014).

  4. Bagley, B. & Nyblade, A. A. Seismic anisotropy in eastern Africa, mantle flow, and the African superplume. Geophys. Res. Lett. 40, 1500–1505. (2013).

  5. Reed, C. A., Liu, K. H., Yu, Y. & Gao, S. S. Seismic anisotropy and mantle dynamics beneath the Malawi Rift Zone, East Africa. Tectonics 36, 1338–1351. (2017).

  6. Tepp, G. et al. Seismic anisotropy of the upper mantle below the Western Rift, East Africa. J. Geophys. Res. Solid Earth 123, 5644–5660. (2018).

  7. Silver, P. G. Seismic anisotropy beneath the continents: Probing the depths of geology. Annu. Rev. Earth Planet. Sci. 24, 385–432. (1996).

  8. Tommasi, A., Mainprice, D., Canova, G. & Chastel, Y. Viscoplastic self-consistent and equilibrium-based modeling of olivine lattice preferred orientations: Implications for the upper mantle seismic anisotropy. J. Geophys. Res. Solid Earth 105, 7893–7908. (2000).

  9. Blackman, D. K. & Kendall, J.-M. Sensitivity of teleseismic body waves to mineral texture and melt in the mantle beneath a mid–ocean ridge. Philos. Trans. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 355, 217–231. (1997).

  10. Whitehouse, M. J. et al. Precambrian basement character of Yemen and correlations with Saudi Arabia and Somalia. Precambrian Res. 105, 357–369. (2001).

  11. Ali, M. Y. Petroleum geology and hydrocarbon potential of the Guban basin, Northern Somaliland. J. Pet. Geol. 38, 433–457. (2015).

  12. Ali, M. Y. & Watts, A. B. Subsidence history, crustal structure, and evolution of the Somaliland-Yemen conjugate margin. J. Geophys. Res. Solid Earth 118, 1638–1649. (2013).

  13. Ali, M. Y. & Watts, A. B. Tectonic evolution of sedimentary basins of northern Somalia. Basin Res. 28, 340–364. (2016).

  14. Coffin, M. F. & Rabinowitz, P. D. Reconstruction of Madagascar and Africa: Evidence from the Davie fracture zone and western Somali basin. J. Geophys. Res. Solid Earth 92, 9385–9406. (1987).

  15. Ali, M. Y. & Watts, A. B. Crustal structure of the Nogal basin, northern Somalia. J. Afr. Earth Sci. 184, 104385. (2021).

  16. Bott, W. F., Smith, B. A., Oakes, G., Sikander, A. H. & Ibraham, A. I. The tectonic framework and regional hydrocarbon prospectivity of the Gulf of Aden. J. Pet. Geol. 15, 211–243. (1992).

  17. Ali, M. Y., Geng, M., Fairhead, J. D. & Adan, A. 3D constrained inversion of gravity and magnetic data to image basement and plutonic bodies: A case study from Dood Arale Basin, eastern Somaliland. Geophysics 86, B321–B334. (2021).

  18. Leroy, S. et al. From rifting to oceanic spreading in the Gulf of Aden: A synthesis. Arab. J. Geosci. 5, 859–901. (2012).

  19. Tozer, B. et al. Global bathymetry and topography at 15 arc sec: SRTM15+. Earth Sp. Sci. 6, 1847–1864. (2019).

  20. ArRajehi, A. et al. Geodetic constraints on present-day motion of the Arabian Plate: Implications for Red Sea and Gulf of Aden rifting. Tectonics (2010).

  21. Xu, M. & He, J. Seispy: Python module for batch calculation and postprocessing of receiver functions. Seismol. Res. Lett. (2022).

  22. Nonn, C., Leroy, S., Lescanne, M. & Castilla, R. Central Gulf of Aden conjugate margins (Yemen-Somalia): Tectono-sedimentary and magmatism evolution in hybrid-type margins. Mar. Pet. Geol. 105, 100–123. (2019).

  23. Dugda, M. T. & Nyblade, A. A. New constraints on crustal structure in eastern Afar from the analysis of receiver functions and surface wave dispersion in Djibouti. Geol. Soc. Lond. Spec. Publ. 259, 239–251. (2006).

  24. Hammond, J. O. S. et al. The nature of the crust beneath the Afar triple junction: Evidence from receiver functions. Geochem. Geophys. Geosyst. (2011).

  25. Kumar, U., Legendre, C. P. & Huang, B. S. Crustal structure and upper mantle anisotropy of the Afar triple junction. Earth Planets Sp. 73, 166. (2021).

  26. Reed, C. A. et al. Receiver function constraints on crustal seismic velocities and partial melting beneath the Red Sea rift and adjacent regions, Afar depression. J. Geophys. Res. Solid Earth 119, 2138–2152. (2014).

  27. Ahmed, A. et al. Crustal structure of the rifted volcanic margins and uplifted plateau of Western Yemen from receiver function analysis. Geophys. J. Int. 193, 1673–1690. (2013).

  28. Ahmed, A. et al. Crustal structure of the Gulf of Aden southern margin: Evidence from receiver functions on Socotra Island (Yemen). Tectonophysics 637, 251–267. (2014).

  29. Tiberi, C. et al. Crustal geometry of the northeastern Gulf of Aden passive margin: Localization of the deformation inferred from receiver function analysis. Geophys. J. Int. 168, 1247–1260. (2007).

  30. Watremez, L. et al. The crustal structure of the north-eastern Gulf of Aden continental margin: Insights from wide-angle seismic data. Geophys. J. Int. 184, 575–594. (2011).

  31. Christensen, N. I. Poisson’s ratio and crustal seismology. J. Geophys. Res. Solid Earth 101, 3139–3156. (1996).

  32. Gao, S. S., Liu, K. H. & Abdelsalam, M. G. Seismic anisotropy beneath the Afar Depression and adjacent areas: Implications for mantle flow. J. Geophys. Res. Solid Earth (2010).

  33. Zhang, S. & Karato, S.-I. Lattice preferred orientation of olivine aggregates deformed in simple shear. Nature 375, 774–777. (1995).

  34. Mainprice, D. & Silver, P. G. Interpretation of SKS-waves using samples from the subcontinental lithosphere. Phys. Earth Planet. Inter. 78, 257–280. (1993).

  35. Silver, P. G. & Savage, M. K. The interpretation of shear-wave splitting parameters in the presence of two anisotropic layers. Geophys. J. Int. 119, 949–963. (1994).

  36. Wüstefeld, A., Bokelmann, G., Barruol, G. & Montagner, J.-P. Identifying global seismic anisotropy patterns by correlating shear-wave splitting and surface-wave data. Phys. Earth Planet. Inter. 176, 198–212. (2009).

  37. Kumar, U. & Legendre, C. P. Crust-mantle decoupling beneath Afar revealed by Rayleigh-wave tomography. Sci. Rep. 12, 17036. (2022).

  38. Accardo, N. J. et al. Surface wave imaging of the weakly extended Malawi Rift from ambient-noise and teleseismic Rayleigh waves from onshore and lake-bottom seismometers. Geophys. J. Int. 209, 1892–1905. (2017).

  39. Ayele, A., Stuart, G. & Kendall, J.-M. Insights into rifting from shear wave splitting and receiver functions: An example from Ethiopia. Geophys. J. Int. 157, 354–362. (2004).

  40. Kendall, J.-M. et al. Mantle upwellings, melt migration and the rifting of Africa: Insights from seismic anisotropy. Geol. Soc. Lond. Spec. Publ. 259, 55–72. (2006).

  41. Al-Lazki, A. et al. Upper mantle anisotropy of southeast Arabia passive margin [Gulf of Aden northern conjugate margin], Oman. Arab. J. Geosci. 5, 925–934. (2012).

  42. Kaviani, A. et al. Mantle-flow diversion beneath the Iranian plateau induced by Zagros’ lithospheric keel. Sci. Rep. 11, 2848. (2021).

  43. Beyreuther, M. et al. ObsPy: A Python toolbox for seismology. Seismol. Res. Lett. 81, 530–533. (2010).

  44. Ligorría, J. P. & Ammon, C. J. Iterative deconvolution and receiver-function estimation. Bull. Seismol. Soc. Am. 89, 1395–1400. (1999).

  45. Zhu, L. & Kanamori, H. Moho depth variation in southern California from teleseismic receiver functions. J. Geophys. Res. Solid Earth 105, 2969–2980. (2000).

  46. Christensen, N. I. & Mooney, W. D. Seismic velocity structure and composition of the continental crust: A global view. J. Geophys. Res. Solid Earth 100, 9761–9788. (1995).

  47. Efron, B. & Tibshirani, R. Bootstrap methods for standard errors, confidence intervals, and other measures of statistical accuracy. Stat. Sci. 1, 54–75. (1986).

  48. Xu, M., Huang, H., Huang, Z. & Wang, L. SplitRFLab: A MATLAB GUI toolbox for receiver function analysis based on SplitLab. Earthq. Sci. 29, 17–26. (2016).

  49. Dreiling, J. & Tilmann, F. BayHunter: McMC transdimensional Bayesian inversion of receiver functions and surface wave dispersion. GFZ Data Serv. (2019).

  50. Rümpker, G., Kaviani, A. & Latifi, K. Ps-splitting analysis for multilayered anisotropic media by azimuthal stacking and layer stripping. Geophys. J. Int. 199, 146–163. (2014).

  51. Latifi, K. et al. The effect of crustal anisotropy on SKS splitting analysis—Synthetic models and real-data observations. Geophys. J. Int. 213, 1426–1447. (2018).

  52. Silver, P. G. & Chan, W. W. Shear wave splitting and subcontinental mantle deformation. J. Geophys. Res. Solid Earth 96, 16429–16454. (1991).

  53. Liu, H. & Niu, F. Estimating crustal seismic anisotropy with a joint analysis of radial and transverse receiver function data. Geophys. J. Int. 188, 144–164. (2012).

  54. Silver, P. G. & Chan, W. W. Implications for continental structure and evolution from seismic anisotropy. Nature 335, 34–39. (1988).

  55. Bowman, J. R. & Ando, M. Shear-wave splitting in the upper-mantle wedge above the Tonga subduction zone. Geophys. J. Int. 88, 25–41. (1987).

  56. Wüstefeld, A., Bokelmann, G., Zaroli, C. & Barruol, G. SplitLab: A shear-wave splitting environment in Matlab. Comput. Geosci. 34, 515–528. (2008).

  57. Kennett, B. L. N. & Engdahl, E. R. Traveltimes for global earthquake location and phase identification. Geophys. J. Int. 105, 429–465. (1991).

  58. Zheng, T., Chen, L., Zhao, L. & Zhu, R. Crustal structure across the Yanshan belt at the northern margin of the North China Craton. Phys. Earth Planet. Inter. 161, 36–49. (2007).

  59. Wessel, P. et al. The generic mapping tools version 6. Geochem. Geophys. Geosyst. 20, 5556–5564. (2019).

Download references


We are grateful to the Somaliland Ministry of Minerals and Energy for approving the deployment of the seismometers and acquiring gravity data. Special thanks go to Abdirahman Jama Mohammed and Abdillahi Ismail Hasan Mataan for assisting with the fieldwork. The software packages Obspy, SplitRFlab, BayHunter, Seispy, and Generic Mapping Tools (GMT) were used to perform the analysis and prepare geographic plots 21,43,48,59.

Author information

Authors and Affiliations


M.Y.A.: Conceived the project, funded and supervised the research, acquisition of the data, data processing, results validation, coordination between the co-authors, drafting the paper, and final submission. M.I.: Data analysis, figure preparation, writing the first draft of the paper. I.M.Y.: Acquisition of the data, data analysis, assisting in writing the manuscript. A.K.: Advising during the research, and assisting in writing the manuscript.

Corresponding author

Correspondence to
Mohammed Y. Ali.

Ethics declarations

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

About Scientific Reports

Scientific Reports is an open-access journal publishing original research from across all areas of the natural sciences, psychology, medicine, and engineering.

This site uses Akismet to reduce spam. Learn how your comment data is processed.

error: Content is protected !!