The key role of conjugate fault system in importing earthquakes into the eastern flank of the Red Sea

This study aims to synthesize seismic observations with gravity and magnetic data and to suggest a new scenario on the development of the Harrat Lunayyir (HL) tectonic system on the eastern Red Sea coastline, Saudi Arabia. Gravity and aeromagnetic anomalies distinctly mapped the NE and NW trends, while the InSAR data depict a small NW–SE graben and an NW–SE dyke. High-resolution relocations, which are well-consistent with the focal mechanism solutions for events with magnitudes greater than 3.0, admit two distinctly fault styles of different orientations. Thus, leading to the NE and NW fault planes’ reactivation related to the Precambrian basement faults and the Red Sea rift system, respectively. The spatiotemporal distributions of epicenters and focal mechanism solutions suggest a new seismic deformation scenario of the 2009 earthquake seismic activity. The low static frictions of 0.2–0.35 obtained from the stress inversion indicates reactivation of preexisting faults in the respective seismogenic zones. The obtained results give rise to a swarm-like sequence of tectonic implications, two activated fault styles differently oriented, and an NE conjugate fault system inherited in the region, which plays a vital role in transferring the ambient stress regime into the Red Sea’s eastern flank.


Introduction
The current geodynamic process acting in the Red Sea is characterized by the extensional tectonic stress regime as accommodated by the Arabian plate's divergence movement away from the African plate (Stoeser and Camp 1985;Abdelfattah et al. 2020). Within the Red Sea's eastern flank, the stress regime is presumably imported by a shallow active asthenosphere (Abdelfattah et al. 2017b); forming a complex structure consisting of several fault styles differently orientated, such as the NW-SE dykes and the NNW normal faults. Seismological studies and geodetic observations revealed that the Arabian Shield is undergoing an NE-SW extensional stress regime, presumably prompted by the Red Sea rift system , 2012ArRajehi et al. 2010). However, recent studies of Abdelfattah et al. (2017bAbdelfattah et al. ( , 2020 showed incompatibility between the large-scale regional tectonic stress in the Red Sea opening and suggested that the stress is transferring into the Red Sea's eastern flank across the conjugate fault system in the region. The active continental lithosphere of the northernmost part of the Red Sea's flanks is tectonically manifested by inland earthquakes occurring in several seismic dislocation zones (Fig. 1). It is worth knowing that relatively low temperatures characterize the earthquake occurring in the brittle volume of crust, which may not be valid in the volcanic regions, where arguments arise due to the coupling between the tectonic and magmatic processes. Owing to the fault geometries can be adapted to reactivate in tectonic stress status, it is crucial to know the reactivation potential of tectonic faults in the region with tectonic and magmatic processes.
One of these regions is the Harrat Lunayyir (HL) enriched by the active volcanic field. The HL is a part of the Arabian Shield that comprises the Red Sea's eastern flank. It attracts global interest because of the 2009 earthquake episode of the largest event with magnitude Mw 5.4, preceded by thousands of events. Since this seismic event, scientists debate whether the HL's earthquakes were due to tectonic or magma upwelling through the crust. The essential features of HL include: (1) the presence of Cenozoic volcanic lavas over the Precambrian crystalline rocks; (2) the earthquake episode occurred in 2009 inducing more than 30,000 events within 3 months (April-June 2009) before and after the mainshock occurred on May 19 with magnitude M5.4 (SGS catalog). This seismic activity is called Al-Ays earthquake, since Al-Ays is the town's name located 40 km southeast of the mainshock's epicenter. Al-Ays' town was affected by minor damage evaluated with VI degrees' intensity in the Mercalli scale (Al-Amri and Fnais 2009). The seismic activity at HL attracts global interest, debating its origin. In general, the seismic activity could be tectonic, magmatic, or both of them. For HL, most of the studies published after the earthquake of 2009 stated that the earthquake swarm of HL is controlled mainly by a magma rising and emplacement of shallow depth dykes (Abdelfattah et al. 2017a;Al-Amri et al. 2012;Al-Zahrani et al. 2012Koulakov et al. 2014;Mukhopadhyay et al. 2013;Pallister et al. 2010;Zobin et al. 2013). Such a scenario is proven in other regions, such as the Afar area (East African rift), where seismicity is associated with active dykes' intrusion with clear evidence. However, there is no evidence of an imminent volcanic eruption in the HL area after 10 years of the main-shock of 2009. The last eruption in HL occurred in the tenth century (Siebert et al. 2011). There has been a significant doubt in the primary model stating that a magma upwelling caused HL 's seismicity through the crust. The key question on the origin of the seismic activity, that occurred in 2009 in the area, remains. Recently, Rashed et al. (2020) re-examined and discussed evidence of the dyke intrusion model adopted by several authors. They favor the tectonic origin. This study contributes to this debate on the cause of the HL area's seismic activity.
To explore the change or redistribution of the tectonic stress conditions in the HL active inland zone, we re-investigate the evidence associated with the 2009 HL seismic activity based on gravity data analysis. We use edge detector techniques combined with the spatiotemporal distribution of epicenters, focal mechanism solutions, and the result of aeromagnetic anomalies. This study provides an impetus to characterize adequate delicate fault structures and the current status of tectonic stress controlling the earthquake generations and adapting the fault geometries that might be reactivated due to the tectonic stress redistribution accumulations. In this study, we attempt to decipher the following: 1 The controversy on whether the 2009 earthquake episode originated due to stresses driven tectonic or magmatic processes. 2 The activated delicate fault structures.
3 The key role of the conjugate fault system to importing earthquake activities in the Arabian Shield.

Geological setting
The Red Sea's eastern margin is represented by large Precambrian basement rocks (600-900 Ma) exposed on the surface with little or no sediment cover known by the Arabian Shield. The Arabian Shield is considered as a Juvenile Neoproterozoic crust. Its structures arise from the convergence, and collisions between continental fragments occurred between 630 and 550 Ma (Stern and Johnson 2010). During Cenozoic volcanic, the Arabian Shield is tectonically accommodated by a divergence between the Arabian and African plates that resulted in the Red Sea's opening. There is evidence that the Precambrian rocks of the Arabian Shield have been under several tectonic phases. Among the significant structures in the Arabian Shield, the NW-SE normal faults and dykes, are recognized to be parallel to the Red Sea and NE-SW strike-slip faults, respectively. In HL and its surrounding, NW-SE faults and dykes are described by Bosworth (2015). The dykes and normal faults are emplaced during the early stages of the Red Sea rift (Stern and Johnson 2019). Moreover, the Arabian Shield is scattered by Cenozoic volcanic rocks known as Harrats, covering a vast area (about 180,000 km 2 ). This Cenozoic volcanic activity is related to the regional extension tectonic event that started 30 Ma ago and continues up to the present, resulting in the Red Sea rifting (Coleman et al. 1983). The young Harrats in Arabian Shield includes Harrat Bouqum, Kishb, Rahat, Khaybar, and HL. The HL, as a case study, is the smallest volcanic field. It is 60 km from the Red Sea coast and is considered as part of the eastern margin of the northern Red Sea. The HL shows several cones trending along NW-SE (Al-Amri et al. 2012), with the latest eruption occurred about 1000 years ago. Figure 2 shows the geologic map of HL and the surrounding area. The figure shows different igneous rocks that compose the Precambrian basement exposed at the surface with remarkable NW-SE and NE-SW trends. Several surface faults are mapped by Jónsson (2012), and most of them are trending along NW-SE direction, parallel to the Red Sea rift. They are considered normal and tensional fractures. The geologic map also shows the recent volcanic rocks observations related to the latest eruption in the tenth century. The HL field has a random shape without any preferential trend. All Cenozoic volcanic rocks in the Arabian Shield, including HL, are assumed to be caused by hot material in the upper mantle (Camp and Roobol 1992).

Gravity anomaly analysis
Bouguer gravity anomalies We used the high resolution (1' × 1') global gravity model (Sandwell et al. 2014). The global model is highly useful in revealing regional crustal structures. Figure 3 shows the complete Bouguer anomaly map of the HL area. The map shows long and short wavelengths due to different subsurface densities with a general increase from − 70 to 50 mGal toward the West. In general, the long-wavelength components represent deep sources. In our case, it corresponds to the crust/mantle discontinuity. Indeed, the general increase toward the West fits the crust structure of passive margin, revealed by the thinning of the crust toward the Red Sea. The crustal thinning is described by Al-Damegh et al. (2005) and Hansen et al. (2007). To efficiently observe the crustal heterogeneities' effect, we extract a regional anomaly representing the effect of crust/mantle discontinuity from the observed anomaly. The regional anomaly representing the crust/ mantle discontinuity must be a smooth surface obtained from the polynomial method, since the discontinuity is deep (about 30 km). The extent of the study area is limited to some hundred's kilometers. We calculate different polynomial surfaces using the observed Bouguer values. The degree of the best polynomial surface is determined using the method proposed by Gabtni and Jallouli (2017). The best polynomial surface, fitting the long-wavelength considered a regional, is the third-degree polynomial surface (Fig. 4a). This regional anomaly is subtracted from the observed complete Bouguer anomalies to obtain a residual. Figure 4b shows the residual anomaly map obtained by subtracting the regional anomaly from the observed Bouguer anomaly. It represents the effect of different densities within the uppermost crust of the study area, including density variation within the Precambrian basement related to crustal heterogeneities, dykes, or faults, To efficiently enhance these anomalies, we calculate the vertical gradient of the residual anomalies. This process amplifies all anomalies, especially anomalies of high frequencies. It also aids in the separation of anomalies. The increase in frequency leads to an increase in the amplification using this process. The vertical gradient is calculated in spectral domain using the Generic Mapping Tools developed by Wessel and Smith (1998). It is obtained by transforming the field to frequency domain, multiplying the field in frequency domain by the radial wave number (k r ) and then transforming back the obtained signal to space domain. Figure 5 shows the resultant vertical gradient map. It demonstrates prominent NE-SW anomalies interfered with short NW-SE wavelength anomalies. These latter trends were not distinguishable in the residual map, and they appear in the vertical gradient map. The anomalies require more accurate and higher sample frequency of gravity data, since they are of short wavelengths. These NW-SE short-wavelength anomalies correspond to the known NW-SE faults and dykes, parallel to the Red Sea rift. The structures are well described and evidenced by aeromagnetic data Rashed et al. 2020). However, the repetitive of prominent positive and negative gravity anomalies trending in the NE-SW direction is not common. These anomalies have a width of about 30 km and an amplitude of 10 to 20 mGal. This is why the geological map ( Fig. 2) shows structures and trends represented by dykes, normal faults, or thin sediment layers trending in NW-SE direction related to the Cenozoic extension event and the Red Sea's opening. Gravity anomalies hardly express these Cenozoic structures. In addition, no anomaly fits the recent HL volcanic rocks. The NE-SW gravity anomalies are dominant in the residual and vertical gradient maps. They may override the gravity signatures of the known NW-SE structures.

Residual anomaly analysis
A key question arises here on the cause of these NE-SW prominent gravity anomalies (Figs. 4b and 5). These anomalies are related to the Precambrian basement of the Arabian Shield based on the wavelength (10 km), amplitudes (tens of mGals), extent (50 to 100 km), and the geological context of the study area. The Arabian Shield has been under several tectonic phases. Its structure occurs as a result of collisions between continental fragments. The so-called basement implicates many heterogeneities. Therefore, we consider that the observed NE-SW prominent gravity anomalies are caused by changes in composition and densities within the Precambrian basement heterogeneities.

Boundary analysis
The observed NE-SW residual anomalies associated with crustal heterogeneities are significant. These heterogeneities correspond to an arrangement of different basement blocks with different compositions and densities. Based on the geological context of the region, the boundaries between blocks should be inherited faults, initiated by collisions of continental fragments of crust that occurred between 630 and 550 Ma (Stern and Johnson 2010). We applied edge detector filters based on gradients of potential field anomalies to detect basement blocks' boundaries with different densities. These filters aid the location of the edge of sources, faults, dykes, densities, or change of magnetic susceptibility using gravity and magnetic anomalies. Cordell and Grauch (1985) initiated these techniques. These techniques were developed by Blakely and Simpson (1986) and applied by several authors (Jallouli and Mickus 2000;Jallouli et al. 2013). In this study, we calculate the Magnitude of the Horizontal Gradient (MHG) to locate the edges of the crustal block: where ( ∂U ∂x ) 2 and ( ∂U ∂y ) 2 are the horizontal gradients of the potential field (U) along x and y directions, respectively. Figure 6a shows the MHG of the residual gravity anomalies. Peaks or Maxima of MGH correspond to the edges of the crustal blocks with different densities. This map shows that most of the gravity trends are in the NE-SW to NNE-SSW direction. We apply a directional gradient to better illustrate these trends. Figure 6b shows the horizontal gradient along the NW-SE direction (N135 o E) almost perpendicular to the main trends. This directional horizontal gradient allows enhancing NE-SW gravity trends corresponding to the edge of the crustal blocks. In other words, it will enable better detection of the NE-SW and NNE-SSW boundaries in the basement, indicated by maxima and minima of the obtained directional horizontal gradient (Fig. 6b).

Aeromagnetic anomalies analysis
The Aeromagnetic data cover the HL and its surrounding areas. Several authors described and analyzed the magnetic anomalies over the HL area Rashed et al. 2020). These data allow us to map anomalies of high frequencies, higher than those expressed by Bouguer data. The aeromagnetic maps show different narrow NW-SE trends (Fig. 7a). These trends correspond to a set of faults and thin dykes. We remind that these dykes evidenced by aeromagnetic data are not or barely shown in gravity maps. There is a complementarity between available gravity and magnetic findings.
Moreover, all magnetic analyses illustrate the clear magnetic signature of the recent HL lava sheets that emplaced in the tenth century. Anomalies of high frequencies and amplitudes express these recent lavas of HL. The aeromagnetic map (Fig. 7a) shows also NW-SE and NE-SW trends. Based on the aeromagnetic anomalies of HL area (Fig. 7a), we mapped the main magnetic trends. The obtained map (Fig. 7b) shows a predominance of NW-SE trends corresponding to faults and dykes. These dykes are aged 20-24 Ma, and they are associated with the beginning of the Red Sea rift (Stern and Johnson 2019).

Seismological aspects
A relationship between the earthquake locations and nearby geological faults is essential to manifest the fault structures. Moreover, an earthquake's fault plane solutions provide information on the aligned fault structures and the current status of tectonic stress that controls earthquake generations in an earthquake-prone area. The precise locations and focal mechanism solution play a crucial role in this study. We determine the 2009 earthquake episode's accurate hypocentral relocations using a scheme of the absolute location described later. In addition, we determine the focal mechanism solutions using the moment tensor technique based on the double-couple source approximation. Figure 1 shows the seismic stations' map and the epicenters for events of magnitudes greater than three (Abdelfattah et al. 2017a). The seismographic network comprises 23 seismic stations equipped with three components, distributed around the epicentral area. Each seismic station was equipped with a 1-Hz seismometer, with continuously recorded at a sampling rate of 100 Hz, and a GPS receiver to maintain an internal clock's accuracy with an order of 1 ms.

Hypocentral locations
The data set used to relocate the hypocenters belonging to the 2009 HL earthquake swarm comprised 136 events, representing the significant events with magnitudes three at least to provide high accuracy picking of the arrival times. The arrival times of both P-and S-waves were manually picked from time-series recorded by vertical and two horizontal components, respectively, using the GSAC software program (Herrmann and Ammon 2004). The one-dimensional velocity model of the P-wave was estimated by Rodgers et al. (1999). The S-wave model  Table S1. The 2009 HL seismic activity was described as a swarm-like earthquake sequence with a pattern that began with a swarm activity denoted by the events 1-35 (triangles) that delineated the NE fault trend; the foreshock-mainshock stage represents the events 36-109 (circles) that are located over the NNW fault trend and the aftershock stage represents the events 110-136 (diamonds) that distributed over the NNW fault trend. The red dashed line represents the fault trends identified from the distribution of epicenters was assumed by multiplying a scaling factor of 1/ √ 3 by the P-wave velocities. We determined the earthquake's hypocenters using the ELOCATE code in GSAC software (Herrmann and Ammon 2004). The relocation was improved by incorporating station delays, measured using the average difference between the observed and theoretical travel times. The iteration was terminated when the root mean square (RMS) of the residual arrival times was less than 0.01 s. The average spatial errors were 137 m and 540 m in the horizontal and vertical direction, respectively. Figure 8 shows the distributions of relocated epicenters obtained from this study, which distinctly delineated two fault segments of different orientations. The figure shows that the NE striking fault segment was initially activated due to shear stress that was incompatible with the region's ambient tectonic stress. This fault segment was developed to the eastward before stopped by the prevailing NW dykes in the region. Consequently, the stress accumulation was continued and redistributed to be accommodated along the NW fault, where the main rupture occurred.

Focal mechanism solutions
One of the robust tools that stimulate the extraction of appropriate information on the aligned structures is the earthquake focal mechanisms. In the analysis, we retrieve the focal mechanism solutions using the waveform moment tensor inversion technique of Ichinose et al. (2003) based on the double-couple source point approximation. The technique is employed to avoid getting multiple focal mechanism solutions of small-size earthquakes whenever the number of seismic stations recording events is deficient and azimuthal coverages are inadequate. In applying the waveform inversion in retrieving the focal mechanism solutions, two seismic stations with six components may be adequate to resolve reasonable focal mechanism solutions.
We extract the corrected ground velocity waveforms after removing the linear trend and deconvolving the seismograms' instrumental response. The second-order bandpass filter with corner frequencies in the range of 0.03-0.1 Hz was applied to remove complexities inserting into the waveforms due to the medium heterogeneities along source-to-receiver paths, resulting in simplifying the low-frequency signals of the surface and S waves. We employ the one-dimensional velocity model of Rodgers et al. (1999) to calculate the Green's functions based on the fast reflectivity and f-k summation computer code of Zeng and Anderson (1995). The match between the observed waveforms and optimum synthetic ones was measured based on the variance reduction's function. Figure 9 shows an example of waveform inversion representing the comparison between the observed and synthetic seismograms computed for the optimum solutions of double-couple approximation. The focal mechanism solutions of the analyzed events are shown in Fig. 10 and listed in Additional file 1: Table S1.
Additional file 1: Figure S1 shows the distribution of plunge angles for compression and tension axes based on the focal mechanism solutions obtained from the waveform inversion. The majority of earthquakes suggest that the prevalent mechanisms are normal to oblique faults, reflecting two fault segments consistent with that delineated from epicenters' distribution. Additional file 1: Figure  S2 shows a histogram showing a noticeable trend of the T-axis orientations for an extensional stress regime, which is oriented on the NE striking, parallels to the ambient tectonic stress acting in the Red Sea. On the other hand, a wide range distribution of the P-axis illustrates various focal mechanisms ranging from strike-slip to normal faulting (Additional file 1: Fig. S2). The second extensional stress regime is accommodated on the NW fault plane that reactivates the oblique faults of the region's sinistral strike-slip movement. Accordingly, the stress accumulation mechanism that led to the 2009 seismic activity indicates that the earthquake generations are regulated by the regional tectonic stress acting due to the Red Sea rift.
The iterative joint inversion for stress and fault orientations of Vavryčuk (2014) was applied to compute the tectonic stress in the HL seismogenic zone using focal mechanisms of 136 events listed in Additional file 1: Table S1. Three angles of the orientations of the principal stresses (σ 1 , σ 2 , and σ 3 ) and the shape (stress) ratio (R) are estimated using the following equation; The inversion of Vavryčuk (2014) is based on evaluating the fault instability (I) that is measured along the range of 0 to 1 to resolve the fault plane ambiguities using the following equation: where τ and σ are the normalized shear and normal stresses and μ is the friction coefficient. To evaluate the uncertainties in the principal stress directions, shape ratio and in the static friction, we repeated the stress inversion of focal mechanisms contaminated by artificial noise. The number of random noise focal mechanisms and the mean deviation were set to be 100 realizations and 5°, respectively. To evaluate the fault instability in Eq. 2, a value of friction is required. To resolve the value of friction that often ranges between 0.2 and 0.8, we used several values of friction range between 0.1 and 1.0 with a step of 0.05 and adopt the value which produces the highest overall fault instability for the data inverted. In Mohr's circles, the inclination of the line depends on the value of friction, where the angle is given by tan −1 (μ). The results of the stress inversion are displayed in Fig. 11a-c, depicting the stereogram of the principal stresses with quite wide-spread of uncertainties, where the confidence region is sensitive to the accuracy of focal mechanism solutions used in the inversion. The distribution of principle stresses in Fig. 11a show a wide range of the plunges of σ 1 and σ 2 stress axes. The positions of fault planes in Mohr's circle diagram (Fig. 11a) are characterized by high shear stress. The Mohr's diagram displayed Fig. 9 Plot showing the focal mechanism solutions and a comparison between the observed and synthetic seismograms computed for the deviatoric solution as derived from the waveform inversion Fig. 10 Fault plane solutions derived by the waveform inversion for the events analyzed in this study. The size of beach-balls is scaled in the corresponding magnitudes. The numbers above each beach-ball denote to the successive event number and to the corresponding magnitudes as listed in Additional file 1: Table S1. The spatial distribution of focal mechanisms shows the consistency with the alignments of hypocenters and the nodal planes that showed sinistral strike-slip movement on the NE fault trend and normal faulting mechanism on the NW fault trends Fig. 11 Iterative joint inversion and fault orientations for the event listed in Additional file 1: Table S1 (a) for all events 1-136, (b) for events 1-35 that represent the activity on the NE fault trend and the foreshock stage, and (c) for events 36-136 that represent the activity on the NW fault trend and the foreshock-mainshock and aftershock stages. N is the number of realizations of random noise for accuracy estimates and R is the shape (stress) ratio, respectively, as shown in the shape ratio plot. Parameters σ 1 , σ 2 , and σ 3 are of the principal stresses as shown in stress and Mohr plots. In Mohr's circles the red lines are inclined depending on the value of friction, representing by tan −1 (μ) that the fault planes are distributed in the upper and lower half-planes, implying that the fault planes and its conjugate ones were reactivated. Furthermore, the low value of the shape ratio of 0.1 ± 0.3 (Fig. 11a) indicates significant variability in the focal mechanism solutions, where the σ 1 and σ 2 axes are not well constrained in the plane perpendicular to the σ 3 axis, showing that the division into two groups is necessary to explain the data. The orientation of principal stress axes σ 1 , σ 2 , and σ 3 is retrieved to be 286°/69° ± 24°, 153°/15° ± 25°, and 59°/15° ± 2° for azimuth and plunge, respectively. For the foreshock stage, Fig. 11b demonstrates the stress inversion for the events 1-35 that revealed sinistral strike-slip focal mechanisms. The orientation of principal stresses reflects an oblique focal mechanism of major compressive stress (σ 1 ) of 330° ± 4° trend and an approximately shallow plunge of 27° ± 4° and minor compressive stress (σ 3 ) with shallow plunge of 13° ± 3° and 66° ± 3° NE trend, respectively. For the foreshock-mainshock and aftershock stages, Fig. 11c shows the stress inversion results for the events 36-136 that exhibited normal fault plane solutions. The orientation of principal stresses reflects normal focal mechanism that characterized by major compressive stress (σ 1 ) of 87° ± 3° plunge and 229° ± 3° SW trend and minor compressive stress (σ 3 ) with shallow plunge of 3° ± 3° and 60° ° ± 3° NW trend, respectively. The most conspicuous feature is the low value of the static friction of 0.2 ~ 0.35 obtained in the stress inversion for the foreshock stage and the foreshock-mainshok stage as well as the aftershock stage. The low value of static friction indicates a reactivation of a preexisting faults that are widely distributed in the seismogenic zone of HL.

Discussion
The integration of gravity, aeromagnetic, and seismic data gives an impetus to map the tertiary structures in the HL area represented by the NW dykes and the NE fault trends. The NW trends presumably correspond to dykes and normal faults parallel to the Red Sea axial rift. However, the NE trends are interpreted as inherited conjugate faults accommodated by convergences and collisions during Precambrian.
The NE trends are well delineated by gravity anomalies, implying a heterogeneous distribution of distinctive densities in the upper-most crust and indistinctive magnetic susceptibility contrasts. On the other hand, the NW trends are well detectable in the magnetic anomaly map (Fig. 7a). They are not visible in the gravity anomaly map (Fig. 4b), implying a high amplitude of magnetic anomalies and low differentiation in density contrasts. A close examination of the gravity map shows some rough NW trends related to dykes. Such trends could be better expressed using higher resolution gravity data. A distance of one kilometer between measurement stations should well-express these dykes and faults. There is a complementarity between gravity and magnetic findings based on the available data. Gravity anomalies provide evidence of crustal heterogeneities, and NE inherited faults. In contrast, magnetic anomalies provide NW trends evidence due to dykes and faults related to the Red Sea's opening. This is a good example, where gravity and magnetic finding are different and complementary.
We precisely relocated earthquakes for events having magnitudes greater than three. The change of the seismic activity between the NW and NE faults is the key feature of the seismicity pattern. The spatiotemporal analysis incorporated with the focal mechanism solutions illustrates three stages. The first stage represents the seismic activity preceded 6 days before the maximum magnitude event that mostly depicts the NE fault trend mapped with the events 1-35, listed in Additional file 1: Table S1, which are characterized by focal mechanism solutions of strike-slip faulting. The second stage of seismicity represented the central rupture zone and was delineated on the NW fault trend of normal fault plane solutions representing by the events 36-109. The rest of the seismicity represented the last stage and was distributed on the NW fault trend assigned by the events 110-136 that displayed normal fault plane solutions. The HL seismic activity can be explained based on the spatiotemporal distribution of earthquakes, fault structures, and stress regimes. First, the seismic deformation began with a tectonic earthquake swarm activated by a shear failure along an oblique fault, relating to the conjugate fault system inherited in the region. Second, the shear failure produced discontinuity due to one of the preexisting NW dykes distributed widely in the region that reactivated by foreshock-mainshock stage. Finally, we identified the rest of the seismic activity near the maximum magnitude event as aftershocks. The seismic activity represents a dyke cluster within the brittle volume, where magmafilled incubated in the upper crust from the last eruption. The current results revealed strong compatibilities based on the distribution of epicenters, the focal mechanism solutions, the distribution of compressional and tensional axes' plunge angles, and the stress inversion. The findings showed that during the 2009 seismic activity in the HL, the stress process was initially accommodated along the NE fault plane. The NE fault plane's development was prevented by one of the pre-existing NW dykes prevailing in the region. Consequently, the stress accumulation mechanism was transferred to reactivate the NW fault plane by the foreshock-mainshock stage. The NW fault is evidenced by the InSAR data analysis (Baer and Hamiel 2010). The feature of seismicity pattern in the HL is likely related to the model of Hill (1977) that suggested three stages, which initiates by a series of shear failures occurring on conjugate fault planes joining en-echelon offset dykes at oblique angles, representing a tectonic earthquake swarm related to ambient stress regime within brittle crusts. The second stage suggests that the seismic activity along the NNW dyke is oriented parallel to ambient maximum compressive stress, representing an earthquake sequence comprised of foreshocks, mainshock, and aftershocks.
Tectonically, the occurrences of normal strike-slip fault mechanisms are expected in the region due to the tectonic processes associated with the Red Sea rift and the Arabian Shield's counterclockwise rotation relative to the African plate. The sinistral displacements on the fault planes are recognized by the shear stress regime in the Gulf of Aqaba-Dead Sea fault system and the East Anatolian fault. It is noteworthy that there are offsets along the Red Sea axial rift (see Fig. 1), identified as inherited conjugate faults. The NE faults mapped by gravity anomalies can be related to one of these faults. On the other hand, the NW fault plane is related to the Cenozoic faults parallel to the Red Sea axial rift. The gravitational-derived faults, magnetic data, and InSAR data provide evidence for both the faults structures emanating from the Red Sea rift and counterclockwise rotation shear deformation due to the relative motion between the Arabian and African plates. It is worth-knowing that the conjugate fault system, evidenced by the gravity data analysis and inherited in the region, plays a vital role in Arabian Shield. This finding is consistent with the results obtained from recent studies of Abdelfattah et al. (2017bAbdelfattah et al. ( , 2020 for earthquake sequences in the Jizan and Namas seismogenic zones, located in the Arabian Shied southward of the HL seismogenic zone.

Conclusion
This study focused on the HL seismogenic and volcanic zone in the western shield of Saudi Arabia, which has a complex structure comprises dykes, recent volcanic eruptions, and fault segments of various orientations. We performed an integrated analysis based on gravity and seismological data, considering the results obtained from previous studies, such as magnetic anomalies, the model of InSAR data, and earthquake spectral characteristics. Gravity and magnetic data analysis revealed a set of inherited NE to NNE conjugate faults and NW dykes parallel and related to Red Sea rift. The model of InSAR data shows a small NW-SE graben and an NW-SE dyke. The hypocenters and focal mechanism solutions were determined for earthquakes with a magnitude greater than three among the earthquake swarm. Independently, the orientation of the fault segments was shown by the spatiotemporal distribution of epicenters that is consistent with a nodal plane determined from the focal mechanism solutions. The results of potential and seismic data analysis suggest a new scenario for the genesis and the repartition of HL seismicity, implying a vital role of the tectonic and the ancient conjugate faults affecting Arabian Shield's Precambrian basement. We obtained a stress accumulation based on the Red Sea rift, restoring the HL's inherited conjugate NE faults. Consequently, a dissipation of the stress along NW faults and dykes emplaced during rifting. Thus, we obtained the following conclusion: 1. Three stages are assigned to the feature of seismicity pattern in the HL; foreshock stage along the NE fault, foreshock-mainshock stage with the activity on a patch of the NW fault, and aftershock stage on the NW fault. The feature of this pattern implies that the ambient tectonic stress is presumably transferred into the Red Sea's eastern flank through the inherited conjugate fault system with strike-slip faulting mechanism. These faults play a vital role in generating earthquake activities in the Arabian Shield. 2. Consistent with high-frequency contents that characterize the HL seismogenic zone events, the brittle-failure earthquakes provide a powerful means of explaining the 2009 earthquake episode. The earthquake activity is described as a tectonic earthquake swarm, a swarm-like sequence with an exceedingly large mainshock. An overlapping type of the swarm events followed by a typical earthquake sequence consists of foreshocks, mainshock, and aftershocks. The low-frequency events explained a source effect rather than a propagation effect, where the mainshock released most of the stress around the fault plane. The resulting low stress generated only events dominated by the low frequencies just for a while.