THE SOUTH ALBORAN BASIN SEISMIC SWARM IN 2020-2021 AND ITS POSSIBLE RELATION TO THE SOUTHERN DEVELOPMENT OF THE AL-IDRISSI FAULT SYSTEM

1. INTRODUCTION

Earthquake swarms represent a distinct form of seismic activity that diverges from the typical mainshock-aftershock pattern or Omori’s law [Sykes, 1970; Hainzl, 2003]. They are often associated with volcanic processes [Aoki et al., 1999; Hayashi, Morita, 2003]. Conversely, non-volcanic earthquake swarms have been observed at plate boundaries [Shimamoto, 1993; Tryggvason, 1973] and within tectonic plates [Lopes et al., 2010]. To comprehend the seismic sources and spatial distribution of seismic events, it is imperative to map seismic swarms. The proposed methodology [Sharon et al., 2020] stands out as highly effective for earthquake modeling, facilitating the quantitative characterization of regional seismicity by generating both earthquake kernel density maps and seismic moment (M₀) kernel density maps [Deif et al., 2023; Sharon et al., 2020].

The recently increased seismic activity has been observed in the Alboran Sea, particularly within the South Alboran Basin, marked by a notable increase in seismic events, including aftershock sequences and clustered seismic swarms. An overview of historical and instrumental seismicity in the Alboran Sea reveals significant seismic occurrences such as the Adra earthquakes [Peláez et al., 2007] in 1804 and 1910 (MSK Intensity VIII–X), the Al Hoceima earthquakes in 1994 and 2004 with magnitudes Mw=6 and Mw=6.3, respectively [Cherkaoui, El Hassani, 2012], as well as the 2016 Mw=6.4 earthquake off the coast of Al Hoceima [Buforn et al., 2017; Medina, Cherkaoui, 2017], considered the most intense seismic event in the region. Positioned along the convergent boundary of the African and Eurasian plates, the Al Hoceima region is considered as the most seismically active area in the southern Alboran Sea [Hamdache et al., 2022]. Recent data from the National Center for Scientific and Technical Research (CNRST) indicates a prolonged seismic swarm in the South Alboran Basin spanning several months in 2020–2021. This study aims to delineate the seismic swarm to discern the spatial clustering of seismic events both at the surface and in depth. Employing the methodology proposed by [Sharon et al., 2020], we quantify the density of seismic events per square kilometer during the swarm period, encompassing approximately two years (2020 & 2021). Additionally, to visualize the spatial distribution in depth, the hypocenters have been mapped in three dimensions to elucidate the tectonic processes in the region where the seismic activity occurred.

2. GEOLOGICAL SETTING

The Alboran Sea, situated in the western Mediterranean region between the Iberian Peninsula and northern Morocco, lies at the boundary the Eurasian and African plates. Bathymetric data sourced from the General Bathymetric Chart of the Oceans (https://www.gebco.net/) illustrate the intricate morphology of the sea floor, accentuated by the presence of mountainous landforms known as the Alboran ridges (AR). To facilitate study, many researchers [Gràcia et al., 2019; Ercilla et al., 2021; d’Acremont et al., 2022] have partitioned the basin into three distinct segments (Fig. 1): the South Alboran Basin (SAB), the East Alboran Basin (EAB), and the West Alboran Basin (WAB). Present-day crustal deformations result from fault systems within the broader framework of NW-SE to NNW-SSE trending convergence (4.5–5.6 mm/yr) between the Nubian and Eurasian plates [DeMets et al., 2015], serving as the primary driver of seismic activity.

The sea floor of the Alboran Sea is characterized by a network of normal and strike-slip faults. The Youssuf fault (YF) network, representing the dextral strike-slip system [Stich et al., 2020], contributes to numerous seismic events in northern Algeria. Additionally, fault systems such as Carboneras (CF), Nekor (NK), and Al Idrissi (AIFS) correspond to sinistral strike-slip faults, with the Al Idrissi fault system obviously associated with significant earthquakes in northern Morocco [Gràcia et al., 2019; Stich et al., 2020].

3. DATASET

Given the large scale of the investigation area, precise earthquake coordinates are essential for effective spatial data processing. While most seismic stations in Morocco are concentrated in the northern region, Spain is equipped with a dense network of seismic stations (https://visualizadores.ign.es/estaciones_sismicas/). For this study, we utilized data from both the CNRST and the Instituto Geográfico Nacional (IGN) to handle spatial data processing and mapping. The seismic events under analysis were located within the coordinates 35–37 °N and 2–8 °W, covering the period from January 1, 2020, to December 31, 2021. Over this timeframe, there were recorded 6.354 seismic events with magnitudes ranging from 1.5 to 5.1 (1.5<M≤5.1) and depths varying from 1 to 200 kilometers. The maximum location uncertainties for these events are less than 10 km horizontally. We also referenced seismogenic zones defined by the Geological and Mining Institute of Spain (ZEIS Database, https://info.igme.es/zesis/) to identify precisely the locations of seismic swarms. While both the CNRST and IGN catalogues were employed for spatial data processing and mapping, the statistical analysis was based solely on the CNRST catalogue, which reports local magnitudes. Detailed descriptions of the methodology and statistical approaches will be provided in the methods section.

Fig. 2 depicts the spatial distribution of seismic events in the Alboran Sea, highlighting clusters of heightened seismic activity concentrated within specific seismogenic zones such as the Alboran Ridge – Central Rif (AR – CR), the Granada Basin (GB), and the Gulf of Cadiz North (GCN).

The focal mechanisms, as shown in Table.1, for earthquakes with magnitudes greater than 3, were acquired from the IGN open access data.

4. SPATIAL DATA PROCESSING

To quantify regional seismicity, we utilized two parameters: earthquake kernel density and seismic moment (M₀) kernel density. The subsequent spatial data processing steps outline the calculation of both parameters. A regional analysis was conducted in horizontal coordinates on a 2D grid with 1 km intervals. For each grid point, both parameters were computed using all recorded events within a 6 km radius. Kernel density estimation was employed as a method to ascertain the spatial distribution via a probability density function, where the distance from a reference point was utilized to weight each event.

The earthquake kernel density parameter, denoted as N_k, was computed by summing the weighted events within a 6 km radius of each grid point (Fig. 3), dividing their total by the sampling area (πr²), and normalizing by the earthquake catalogue duration (T) [Sharon et al., 2020]:

 (1)

where N represents the number of earthquakes within a circle of radius r, d(n) denotes the distance between an earthquake event (n) and the center of the circle, σ represents the standard deviation of the Gaussian function, and T signifies the duration of the earthquake catalogue. The resulting values are expressed in events per square kilometer per year (events·km^–2·year^–1). The kernel density method allows for a more sophisticated visualization of seismic activity by smoothing the distribution of seismic events, which helps to reveal patterns and trends that might not be evident in a conventional epicenters map.

Our study encompasses a large area spanning over 33000 km². The process of calculating the earthquake kernel density parameter is complex, particularly when considering the duration of the earthquake catalogue. Therefore, the utilization of computational tools is recommended to facilitate this calculation. For our investigation, we utilized the code (it's available at https://fr.mathworks.com/) and implemented several modifications to enhance the accuracy of our results.

To compute the kernel density of the seismic moment, denoted as ρ_Mok, the seismic moment released by each seismic event must be determined using the following relationship proposed by [Hanks, Kanamori, 1979]:

 (2)

Each energy quantity is weighted according to the distance of the seismic event from the center of the circle, just like the earthquake kernel density estimation. The obtained results are expressed in J·km^–2·year^–1.

 (3)

The open-source software Zmap [Wiemer, 2001] was used for statistical analysis of seismic swarm event sizes and spatial distributions in order to deduce changes in local stresses. We used the Gutenberg Richter law to determine the b-value of the seismic swarm sequence, which is expressed as [Gutenberg, Richter, 1944]:

 (4)

where N is the number of events with magnitude larger than M_th. In this study, the b-value of the Gutenberg – Richter relation is used as an indicator of the stress regime [Schorlemmer et al., 2005]. The maximum curvature technique is used to analyze the frequency-magnitude data, and the maximum likelihood approach is used to estimate the b-value [Aki, 1965]:

 (5)

where  is the average magnitude value, M is the lower limit of the magnitude in the catalog, and ΔM is the catalog binning width. The following relationship describes the uncertainty in estimating the b-value [Shi, Bolt, 1982]:

 (6)

where n_с is the number of earthquakes in the given sample.

The vertical stress is minimal for b values less than 1, indicating a compressive regime [Schorlemmer et al., 2005]. If the vertical stress is intermediate, the value of b is likely to be close to the universal value of unity, indicating a slip fault. If it exceeds 1, the stress regime is extensional [Schorlemmer et al., 2005].

5. RESULTS AND DISCUSSION

As anticipated, Fig. 4 illustrates that the majority of earthquakes are concentrated along the AR – CR seismogenic zone, which includes the Al Idrissi fault system, with a lesser occurrence along the GB seismogenic zone. Moreover, notable clusters are observed in the GCN seismogenic zone, albeit with lower density values compared to the AR and GB zones. The seismic activity shows a more dispersed pattern at the junctions of the CIB and AR seismogenic zones, indicating the intersection of the AIFS and CF fault systems.

These findings reveal the occurrence of two seismic swarms in the Alboran Sea between 2020 and 2021. The most intense swarm, in terms of earthquake frequency, is associated with the Central Rif – Alboran Ridge, while the least intense swarm is linked to the Granada Basin.

According to the findings depicted in Fig. 5, the distribution pattern of seismic moment kernel density in the Alboran Sea closely mirrors that of earthquake kernel density, highlighting pronounced seismic activity zones along the AR and GB. Building upon the preceding observations, it is noteworthy that the locations of the two seismic swarms directly oppose each other. This prompts speculation regarding a potential correlation between these distinct seismic events.

Fig. 6 presents the temporal distribution of seismic events in the GB and AR – CR swarms. The GB swarm shows a sharp increase in seismicity around t=3∙10⁷ seconds, followed by a decline in both the frequency and magnitude of events, indicating a rapid release of stress. In contrast, the AR – CR swarm exhibits continuous seismic activity over time, with frequent moderate to large events. This continuous seismicity suggests a prolonged period of stress release, distinct from the more abrupt behavior observed in the GB swarm.

The presence of two seismic swarms, namely the AR – CR and the GB, positioned opposite each other (face-to-face), may either be coincidental or suggestive of a deeper correlation. A review of literature concerning the Alboran Sea accentuates the intricate nature of the Betic-Rif orogenic system, fueling ongoing debates surrounding diverse geodynamic reconstructions. Prominent models propose subduction processes involving slab roll-back, delamination, slab break-off, or slab tearing, as postulated by [Royden, 1993; Lonergan, White, 1997; Spakman, Wortel, 2004; Jolivet et al., 2008; Vergés, Fernàndez, 2012]. Multiple tomographic studies have delineated a high-velocity structure pervading the upper mantle beneath the Betic Arc and Gibraltar Arc, as highlighted by [Gutscher et al., 2002; Spakman, Wortel, 2004; Bezada et al., 2013]. This anomalous structure displays characteristics resembling a slab and is likely indicative of subducted oceanic lithosphere (Fig. 7).

Seismic activity in the Alboran Sea primarily consists of shallow earthquakes below 40km, as reported by [Buforn et al., 1995; Morales et al., 1997]. Furthermore, intermediate earthquakes (40–150 km of depth) are frequent, predominantly associated with the subducting plate [Santos-Bueno et al., 2019; López-Sánchez et al., 2022]. Deep events, occurring rarely at depths of around 600–650 km, have also been documented, as observed by [Buforn et al., 2011; Chung, Kanamori, 1976]. Given this context, we are unable to offer a definitive explanation for the juxtaposition of these two seismic clusters, AR – CR and GB, in this instance. Consequently, it is advisable to construct a three-dimensional sketch illustrating the distribution of earthquakes at depth to ascertain the most appropriate interpretation.

Fig. 8, a, illustrates the three-dimensional distribution of earthquakes spanning from the Central Rif to the Betic Zone. Examining this diagram reveals a distinct orientation of seismic depth events extending eastward from the Central Rif. Conversely, within the Betic Zone, the structural pattern of earthquakes at depth appears less discernible, likely attributed to a smaller quantity of earthquakes compared to those located in the northern region of Morocco. The focal mechanism of the significant seismic event within the seismic swarm in the AR – CR indicates a left-lateral slip (IGN open access data), with a predominant nodal plane orientation of 14°/82°/–13° (strike/dip/rake). [Lozano et al., 2022] investigated the seismic sequence spanning from 2020 to 2021 in the GB and observed that the hypocenters follow a near-vertical pattern associated with normal faults. [Madarieta-Txurruka et al., 2022] conducted a study on the 2021 seismic sequence in the GB, revealing a concentrated "Chimney-shape" seismic activity resulting from the activation of a normal fault.

Based on the dissimilar orientation of seismic activity at depth between the AR – CR and the GB, along with differences in the tectonic regime of associated faults in these two seismogenic zones, it can be inferred that the positioning of these two seismic swarms opposite each other was coincidental, likely attributed to the high seismicity rate in these two regions. Since the two seismic swarm sequences are unrelated, we have turned our attention to studying the seismic activity in the AR – CR area. Here, we observed that the seismic activity is slanting towards the east. At first glance, this appeared to be related to subduction processes (Fig. 8, a, b), However, existing research indicates that the subduction zone, where intermediate earthquakes typically occur, is situated between longitudes 6–5 W° [Santos-Bueno et al., 2019; López-Sánchez et al., 2022], whereas the observed seismicity is concentrated between longitudes 3–4 W° (Fig. 8, b). These findings are intriguing and highlight the need for further exploration of this seismic swarm to determine if there are nearby structures, possibly associated with the AIFS system, influencing the seismic activity in the eastward direction.

The seismic activity in the CR – AR region seems to be focused mainly on the southern part of the AIFS, as shown in Figs. 4, 5, and 8. Many studies [El Alami et al., 1998; Medina, 2015; Medina, Cherkaoui, 2017; Gràcia et al., 2019] indicate that this specific segment of the fault system experienced significant seismic events during the 2016 crisis, as well as in previous occurrences in 1994 and 2004. Interestingly, there is a lack of strong seismic events in the central and northern segments of the AIFS. This suggests that these segments might be inactive or undergoing slow movement [Gràcia et al., 2019]. Additionally, research suggests that the AIFS has the potential to cause multiple segment ruptures [Gràcia et al., 2019], a phenomenon closely linked to the seismic activity observed in the South Alboran Sea in 2020 and 2021. Density distribution analyses highlight a high concentration of seismic events at the southern end of the AIFS, with significant energy release. To further investigate the cause of the South Alboran seismic swarm, changes in the tension regime during the swarm were monitored by observing the evolution of the b-value.

The b-value calculated for the seismic cluster in the South Alboran Sea is 1.16±0.04 (Fig. 9), which closely aligns with the b-values reported by [Hamdache et al., 2022] for the aftershock sequences following the 1994 event (b=1.01±0.07) and the 2004 event (b=1.04±0.05) by [Telesca et al., 2009]. Since these values are close to one, they suggest compatibility with a slip stress regime, as proposed by [Grob, van der Baan, 2011].

To ensure the accurate estimation of the b-value, we employed moving windows that maintain a consistent number of epicenters. This approach, as proposed by [De Rubeis et al., 1993], tracks the variation of spatial correlation dimension over time. Following numerous tests, we determined that a fixed number of 150 epicenters per window provides precise estimations. Analysis of the b-value variation over time reveals distinct stress regimes (Fig. 10): an initial strike-slip regime with a b-value approximating 1 transitions abruptly to an extensional regime [Schorlemmer et al., 2005]. This abrupt increase in b-value suggests the opening of fractures, thereby possibly supporting the concept of the AIFS southern development proposed by [Gràcia et al., 2019]. Subsequently, these two regimes have alternated. Comparable patterns have been observed in seismic swarms, such as the 2012–2013 seismic swarm in the Eastern Guadalquivir Basin (Southern Spain), where b-values indicate alternating extensional and strike-slip regimes [Peláez et al., 2015].

However, the b-value does not definitively establish the stress regime due to its sensitivity to various influencing factors. Although we utilize the b-value as an indicator of stress variation, a comprehensive analysis, including focal mechanism inversion, is necessary for a more accurate identification of the stress regime.

The focal mechanism data reveal a sequence of seismic events dominated by strike-slip faulting, with occasional evidence of an extensional component. The majority of events are characterized by steep dip angles and rake values near 0° or ±180°, indicating a predominantly lateral motion. However, several events show rake values between –20° and –80°, suggesting oblique mechanisms with a normal (extensional) component (Table 1). This alternation between pure strike-slip and strike-slip with normal faulting may indicate variations in the local stress field, potentially associated with transient stress changes or complex fault geometries (Fig. 11). The consistent presence of these strike-slip mechanisms, with intermittent extensional elements, implies a tectonic regime primarily governed by horizontal shear, with occasional influences from extensional forces likely due to local fault interactions or episodic changes in stress conditions.

The correlation between focal mechanism results and b-value trends appears complex rather than straightforward. The b-value graph shows periods of decreasing b-values, followed by sharp increases. This pattern suggests changes in stress regimes, with drops in b-value potentially indicating increased stress or fault locking, and rises pointing to possible fracture opening or release of stress. Meanwhile, the focal mechanisms show an alternation between strike-slip and strike-slip with normal components, indicative of a transtensional regime. This type of faulting suggests both lateral motion and some degree of extensional force, which aligns with certain rises in the b-value, as these could correspond to stress release events associated with opening fractures or fault segments.

However, the alternation in focal mechanisms does not perfectly mirror the exact timing or scale of b-value fluctuations. Instead, both sets of data collectively point to a variable stress regime but may reflect slightly different aspects of the tectonic process. The b-value is a statistical measure of earthquake magnitude distribution and tends to capture broader stress trends, while focal mechanisms provide a more direct view of the specific faulting style for individual events. Together, these analyses Alboran Ridge the region’s complex tectonic environment, where the interplay of strike-slip and transtensional forces may contribute to the ongoing southern development of the AIFS.

6. CONCLUSION

This study investigates the spatial distribution of seismic events in the Alboran Sea during 2020 and 2021, revealing a heterogeneous pattern with distinct clusters of heightened seismic activity concentrated in key seismogenic zones: the Alboran Ridge – Central Rif (AR – CR), Granada Basin (GB), and the Gulf of Cadiz (GCN). Magnitudes recorded range between 1.5 to 5.1. A seismicity-based approach is employed to pinpoint these clusters, utilizing two parameters: earthquake kernel density and seismic moment kernel density. The analysis identifies two seismic clusters: one situated north of Al Hoceima, specifically in the southern segment of the AIFS, and another in the Granada Basin. Furthermore, it reveals that these clusters exhibit substantial energy release, particularly in the southern part of the AIFS, with values reaching up to 1.6·10¹⁴ J/km^–2·year^–1. 3D seismic modeling conducted in the region spanning from the central Rif to the Betic zone indicates that the predominant seismic activity occurs at shallow depths, typically less than 40 kilometers, while fewer earthquakes are observed at intermediate depths, ranging between 40 and 200 kilometers. The concentration of seismic events in northern Morocco exhibits a depth-oriented eastward pattern, prompting speculation regarding the presence of adjacent geological structures, potentially linked to the AIFS system.

The main seismic event within the seismic swarm had a magnitude of 5.1 and was centered at the southern edge of the AIFS. Analysis of the focal mechanism of this event reveals a left-lateral slip stress regime, with a dominant nodal plane orientation of 14°/82°/–13° (strike/slip/rake). Examining the sequence of seismic activity in the southern Alboran Basin, a b-value close to 1, as derived from the Gutenberg – Richter law, suggests a slip-dominated regime. Further investigation into the variation of the b-value over time exposes distinct stress regimes: an initial decrease in stress with a b-value approximating one, succeeded by a notable rise in the b-value indicating an extensional regime, likely attributed to fracture opening. Subsequently, these regimes alternate. The focal mechanism data reveal a sequence of seismic events dominated by strike-slip faulting, with occasional evidence of an extensional component. The temporal variation of focal mechanisms throughout the seismic swarm shows a shift between strike-slip and strike-slip with normal components, suggesting a transtensional regime. This type of faulting suggests both lateral motion and some degree of extensional force, which aligns with certain rises in the b-value, as these could correspond to stress release events associated with opening fractures or fault segments, although the timing does not perfectly match b-value fluctuations. Together, these findings highlight a complex tectonic environment marked by strike-slip and transtensional forces, supporting the ongoing southern development of the AIFS, as suggested by several researchers [Buforn et al., 2017; Gràcia et al., 2019].

The applied approach proposed by [Sharon et al., 2020] has been used to map and localize the seismic swarm. While this approach has demonstrated success, it remains insufficient to fully comprehend this intricate phenomenon. Thus, for a comprehensive analysis of the seismic swarm sequence, it is highly advisable to employ fractal theory [Mandelbrot, 1989] to quantify both the spatial and temporal clustering of the earthquakes.