Slip Partitioning in the 2016 Alboran Sea Earthquake Sequence (Western Mediterranean)

Introduction

Most earthquakes worldwide show directions of fault displacement close to horizontal or close to vertical, while oblique slip directions are less common (Célérier, 2008). This observation is usually thought to represent the influence of the Earth’s free surface on crustal stresses. Vanishing shear stress at a free surface requires that one principal stress is oriented vertically, thus complicating the occurrence of general, oblique directions of shear, at least on shallow faults (Anderson, 1905). The predominance of strike-slip and dip-slip faulting still holds in many scenarios where regional deformation is actually oblique to the trend of the principal faults or belts of tectonic deformation. Such a behavior was first recognized at the convergent margin of Sumatra, where oblique motion between the Australian and Sunda plates is partitioned into shortening normal to the plate boundary along the subduction interface, and strike slip faulting parallel to the plate boundary in the overriding plate (Fitch, 1972). The case of Sumatra is no exception, and slip partitioning has been analyzed on many plate boundaries since then (e.g., Yu et al., 1993; McCaffrey, 1996). Slip partitioning can also be active during single earthquakes, as demonstrated in an exemplary manner by the complex rupture cascade of the M_W 7.8, 2016 Kaikoura, New Zealand earthquake (Xu et al., 2018; Ulrich et al., 2019), or, previously, by concurrent strike-slip and reverse faulting during the M_W 7.9, 2002 Denali earthquake in Alaska (Eberhart-Phillips et al., 2003; Bemis et al., 2015). Even for very large events, partitioning of slip may occur over very close distances of order of km, as observed along the 400 km long surface rupture from the M_W 7.8, 2011 earthquake on the Kunlun fault in Tibet (Van der Woerd et al., 2002; King et al., 2005).

Slip partitioning for crustal deformation occurs at different length scales, from the order of plate-boundaries down to local structures. Case histories for slip partitioning in smaller earthquakes include the 2014 Northern Nagano earthquake (Japan, M_W 6.2) or the 2018 Hualien earthquake (Taiwan, M_W 6.4), with both events involving strike-slip and reverse motion in settings where fault curvature imposes local variation of the slip direction (Kobayashi et al., 2017; Lo et al., 2019). Local slip partitioning can be frequently inferred from geological observations, for example along the Dead Sea Fault System (Bartov and Sagy, 2004; Gomez et al., 2007), and also has been confirmed from the locations and mechanisms of microearthquakes in many different settings, (e.g., Yamini-Fard et al., 2006; Wenzheng and Hauksson, 2011). By definition, strike-slip faulting is an integral part of slip partitioning, and in general, local slip partitioning is linked to a nonplanar geometry of the faults. Irregularities of the fault trace are common for strike-slip faults, and may arise from the evolution of disconnected faults over time, compositional heterogeneity or fault interaction, among others, (e.g., Mann, 2007). If step-overs or bends are present along a fault, strike-slip motion has to deal with these obstacles, and local deformation will deviate from simple shear by a component of shortening or extension (Fossen and Tikoff, 1998; McClay and Bonora, 2001; Nabavi et al., 2017). Transtensional settings develop around fault releasing bends and extensional step overs, which occur where the sense of fault slip and the sense of the fault offset are the same, like for a right-stepping pull apart on a right-lateral fault. In turn, transpersonal settings develop around fault restraining bends and contractional step overs, which occur where the sense of fault slip and the sense of the offset are opposite to each other.

Here, we analyze the kinematics of an intense earthquake sequence that occurred in 2016 in the southern Alboran Sea (Buforn et al., 2017; Medina and Cherkaoui, 2017; Galindo-Zaldivar et al., 2018; Kariche et al., 2018; Gràcia et al., 2019). Seismicity started with a foreshock, which was a moderate magnitude, M_W 5.1 earthquake on January 21st. This initial event was followed by a M_W 6.3 earthquake at 04:22 UTC on January 25th. Despite its offshore location, about 40 km from the African coast and 125 km from Iberia (Figure 1), this earthquake caused damage along the north Moroccan coast (I = VI) and was widely felt in southern Spain (I = IV, IGN, 2016). During the following months, aftershock activity was intense, including numerous felt, moderate magnitude earthquakes. Moment tensor solutions for this sequence are available from near real-time moment tensor determination at the Instituto Geográfico Nacional (Rueda and Mezcua, 2005) and dedicated studies (Buforn et al., 2017; Gràcia et al., 2019). According to published results, foreshock and mainshock show similar strike slip faulting solutions, while the aftershock series shows more complex characteristics of faulting. In addition to further strike-slip events, pure reverse faulting mechanisms occur, including many of the largest aftershocks. We investigate this case of slip partitioning through relocation and regional moment tensor inversion for further small-to-moderate magnitude aftershocks. We also use moderate aftershocks as empirical Greens functions to evaluate the rupture process of the mainshock and establish the spatial relationship between mainshock and aftershocks. We interpret our results in light of the local faulting geometry, originating transpression during the 2016 earthquake series, and in the regional tectonic context, characterized by present-day transtensional conditions in the Alboran basin.

FIGURE 1. Earthquake locations in the Alboran region from 1990 to 2020 (IGN catalog, www.ign.es), distinguishing between shallow events (yellow dots), intermediate depth events (deeper 40 km, blue dots), and shallow seismicity during 2016 (red dots), including the sequence in the southern Alboran Sea. The principal strike slip faults are labeled (left-lateral Al-Idrisi and Carboneras faults, right-lateral Yusuf fault; all faults from Gràcia et al., 2019). Stars mark the epicenters of the 2016 mainshock and the largest instrumental earthquakes in the Alboran region, including the 1910, M_W 6.1 Adra earthquake (Stich et al., 2003) and the 1994, M_W 5.9 and 2004, M_W 6.3 earthquakes near the city of Al Hoceima (El Alami et al., 1998; Stich et al., 2005). Bathymetry is taken from the EMODnet Bathymetry Consortium (2018). Squares mark selected coastal cities and towns around the Alboran Sea.

Formation and Seismotectonics of the Alboran Basin

The Alboran Sea is the western termination of the Mediterranean and is located on the inner side of the Gibraltar Arc, formed by the Betic and Rif mountain belts in southern Spain and northern Morocco. The relationship between the extensional Alboran Basin and the tightly curved orogenic arc has been analyzed by means of diverse geodynamic models. Besides explanations in terms of subduction-related processes, (e.g., Royden, 1993; Lonergan and White, 1997; Gutscher et al., 2002; Spakman and Wortel, 2004; Gutscher et al., 2012; Faccenna et al., 2014), the configuration of the Alboran-Betic-Rif system also inspired less common geodynamic interpretations, like for example a radial extensional collapse (Platt and Vissers, 1989; Molnar and Houseman, 2004). Recent results from seismic tomography, based on dense station coverage in the Iberia–Maghreb region during several temporal deployments, show clearly the position and geometry of the lithospheric slab in the upper mantle beneath the Alboran region, (e.g., Bezada et al., 2013; Palomeras et al., 2014; Fichtner and Villaseñor, 2015). These images substantially strengthen geodynamic interpretations in terms of subduction related processes. The Alboran slab forms part of the Western Mediterranean subduction system that is characterized by fragmentation of the slab, generalized trench retreat and nearly complete consumption of the remaining Tethys lithosphere since the slowing down of African–Eurasian convergence about 25 Myr ago, (e.g., Faccenna et al., 2014). In particular, the formation of the Alboran Basin is tied to the rapid rollback of the Alboran slab during the Miocene. Receiver functions have imaged a narrow corridor of oceanic lithosphere that is still preserved below the Alboran units at the vertex of the Gibraltar arc, representing the Mesozoic connection between the spreading centers of the Alpine Tethys and the Central Atlantic (Molina-Aguilera et al., 2019; Santos-Bueno et al., 2019).

According to stratigraphy and structural evidence, the most significant extension recorded in the Alboran crust started around the beginning of the Miocene (23 Myr) and lasted until ∼8 Myr in the middle Tortonian, (e.g., Bourgois et al., 1992; Comas et al., 1992; Billi et al., 2011; Giaconia et al., 2015). Kinematic reconstructions suggest that this extensional phase coincides with fast displacement of the Alboran domain, following the retreating Alboran slab for about 500–700 km from a position south of the Balearic Islands at the beginning of the Miocene, to nearly the present position at 8 Myr (Rosenbaum et al., 2002; Vergés and Fernàndez, 2012; van Hinsbergen et al., 2014). The reconstructions corroborate a WSW direction of slab rollback as a requirement to explain the length (>600 km) of the subducted Alboran slab seen in tomographic images. Slab rollback toward Gibraltar and the docking of the allochthonous Alboran units onto the passive continental margins of Iberia and Morocco is also recorded by a change in the geochemical signatures of late Miocene volcanic rocks (12–6 Myr), showing a spatial and temporal transition from subduction-related to intraplate type magmatism (Duggen et al., 2004; Duggen et al., 2005). Since the Tortonian, the Alboran domain slowed down significantly and became affected by slow, oblique convergence between the Iberian and Moroccan margins (e.g., Billi et al., 2011; Giaconia et al., 2015). This fundamental change led to shortening and partial basin inversion (e.g., in the southern Alboran Sea, Martínez-García et al., 2013; Martínez-García et al., 2017; Lafosse et al., 2020). Despite ongoing plate convergence and geological evidence for contraction since the end of the Miocene, the general present-day stress regime in the Alboran domain is transtensional, according to observed geodetic velocities and earthquake focal mechanisms (e.g., Stich et al., 2006; Serpelloni et al., 2007).

Moment tensor focal mechanisms for earthquakes in the Alboran domain show a clear predominance of strike-slip faulting, occasionally including components of normal faulting. The characteristic orientation of tension axes is near ENE–WSW (Figure 2, e.g., Martín et al., 2015; Custódio et al., 2016; Stich et al., 2019). Strike-slip faulting with optional normal components has also been observed for the so far largest instrumental earthquakes in the Alboran region, including the 1910, M_W 6.1 Adra earthquake (Stich et al., 2003a; Stich et al., 2003b) and the earthquakes in 1994 (M_W 5.9) and the 2004 (M_W 6.3) near Al Hoceima (Calvert et al., 1997; Bezzeghoud and Buforn, 1999; Stich et al., 2005; Biggs et al., 2006). In agreement with the dominant faulting style, overall transtensional stress conditions have been inferred for the Alboran region, with stress tensors closer to uniaxial extension in the northern part, and closer to strike-slip condition in the southern part of the Alboran domain (Stich et al., 2006; De Vicente et al., 2008; Olaiz et al., 2009; Soumaya et al., 2018). Compressional principal stress axes show azimuths around N345°E and tensional principal axes around N75°E. Recent transtension is also documented geologically from active faults at the Moroccan margin (d’Acremont et al., 2014; Lafosse et al., 2017; Lafosse et al., 2020). The tectonic regime in the Alboran region is in a marked contrast with the adjacent sections of the plate boundary zone, where transpersonal environments with reverse and strike-slip focal mechanisms are characteristic for northern Algeria as well as for the SW-Iberian margin (e.g., Bezzeghoud and Buforn, 1999; Stich et al., 2003a; Stich et al., 2003b).

FIGURE 2. Regional moment tensor solutions along the Iberia-Maghreb plate boundary zone until 2014 (Martín et al., 2015). We show double couple components of the mechanisms in lower hemisphere projection. The size of the beach-balls is scaled with magnitude (see bottom right). We use EMODnet Bathymetry, as well as source and representation of seismicity like in Figure 1. Strike slip and normal faulting earthquakes are dominant in the plate boundary zone between about 2°W and 5°W.

Extension in the central Betics and Alboran Sea is not straightforward to understand in the context of plate convergence between the Iberian and Maghrebian margins (e.g., Cunha et al., 2012). In fact, the Global Navigation Satellite System (GNSS) velocity field suggests that deformation in the Iberia-Maghreb region is a result of the superposition of two independent processes (e.g., Fadil et al., 2006; Stich et al., 2006; Serpelloni et al., 2007; Koulali et al., 2011; Mancilla et al., 2013; Palano et al., 2013). On regional scale, the Iberia-Maghreb region is subject to oblique Eurasia-Nubia plate convergence in WNW–ESE direction at a rate of ∼4–5 mm/yr. The most striking anomaly in this context is the peculiar SW motion of the Rif mountain range in Morocco with respect to stable Nubia at a rate of ∼3–4 mm/yr (Fadil et al., 2006; Koulali et al., 2011). Motion of the Rif block is aligned with the direction of extension in the central Betics (Mancilla et al., 2013) and Alboran Basin (Spakman et al., 2018). The anomalous emplacement of the Rif block can be modeled by applying a horizontal traction in SW direction to the base of the Rif (Pérouse et al., 2010), such that a local, subcrustal process superimposed on plate convergence can be considered causative for transtensional conditions in the Alboran Sea. A plausible origin of Rif motion is the present stage of slab rollback and detachment, where the mantle resistance to the hanging slab produces traction in direction opposite to the NW-directed absolute plate motion of the Alboran region with respect to the underlying mantle (Spakman et al., 2018; Capella et al., 2020). Strain tensors calculated from GNSS velocities around the Alboran Sea corroborate that a transtensional regime is characteristic of the Alboran Basin today (Serpelloni et al., 2007; Billi et al., 2011; Palano et al., 2013). The Alboran transtensional regime is characterized by ∼2.5 mm/yr of extension in direction of present-day Alboran stretching (ENE–WSW), or ∼3 mm/yr of left-lateral slip along the NNE–SSW faulting directions inferred from strike-slip moment tensors in the Basin (Stich et al., 2006).

The interplay of large-scale plate motion with local subcrustal processes is mapped in the regional fault network by changes in the seismotectonic pattern along the Maghrebian margin from east to west. While the Alboran-Rif domain is characterized by conjugate strike-slip faulting with normal component, seismotectonics along the coast of Algeria and Tunisia is dominated by reverse faulting and right-lateral strike-slip, (e.g., Bezzeghoud and Buforn, 1999; Braunmiller and Bernardi, 2005; Stich et al., 2006; Soumaya et al., 2018). Over the recent decade, marine multibeam and seismic reflection imaging in the Alboran basin was able to characterize the principal active faults in the offshore part of the Alboran domain, (e.g., Gràcia et al., 2006; Gràcia et al., 2012; Martínez-García et al., 2013; Perea et al., 2018). The dominant tectonic feature of the Alboran is the complex fault network of a major left-lateral shear zone that extends for more than 400 km in roughly NE–SW from SE Spain across the Alboran Sea to the eastern Rif Mountains near Al Hoceima (Trans-Alboran Shear Zone, Bousquet, 1979; De Larouzière et al., 1988; Stich et al., 2006). The two principal left-lateral lineaments that cross the Alboran Sea are the offshore part of the Carboneras fault system and the Al-Idrisi fault system, both with ∼90 km length (Figure 1). The right-lateral Yusuf fault connects the Trans-Alboran Shear Zone with the Algerian continental margin. While in general the instrumental seismicity of the Alboran Basin shows little activity of these principal tectonic structures, the 2016 Alboran Sea earthquake is an exception from this rule. The 2016 earthquake sequence can been assigned to the Al-Idrisi fault, and the mainshock appears to nucleate right at a pronounced bend in the fault system (Buforn et al., 2017; Medina and Cherkaoui, 2017; Kariche et al., 2018; Gràcia et al., 2019).

Relocation and Regional Moment Tensor Inversion for the 2016 Alboran Sea Earthquakes

Source properties of small to moderate earthquakes have to be inferred from local and regional observations, because at far distances the signal of the elastic wavefield falls below noise level. We collect the recordings from 35 seismic broadband stations in Spain, Morocco and Portugal, including 30 near regional stations at distances less than 350 km, as well as five intermediate distance stations (350–500 km) to complement the dataset with additional observations in southwest, northwest and north directions from the earthquakes (Figure 3). In order to understand faulting processes during the 2016 Alboran Sea earthquake sequence, we perform time domain moment tensor inversion. Before moment tensor inversion, we relocate the mainshock, foreshock and major aftershocks. Relocation is motivated mainly by reported catalog depths, in excess of 20 km for many events of the sequence (e.g., IAG, iagpds. ugr.es, IGN, www.ign.es), which appears unrealistic in the extended continental crust of the Alboran Basin (e.g., Comas et al., 1999; Torne et al., 2000; Booth-Rea et al., 2007). Shallow source depths are expected for all events according to previous estimates for the seismogenic layer thickness (Fernández-Ibáñez and Soto, 2008; Grevemeyer et al., 2015). Relocation is performed by a grid-search algorithm (NonLinLoc, Lomax et al., 2000), involving a preliminary 3D Earth model for the region (Martín et al., unpublished manuscript). The Earth model was mainly built from available studies of crustal thickness, (e.g., Mancilla and Díaz, 2015), sediment thickness (e.g., Soto et al., 2008) and crustal structure inferred from seismic transects, (e.g., Gil et al., 2014).

FIGURE 3. (A) Distribution of 35 regional seismic broadband stations (triangles) used for waveform analysis of the 2016 Alboran sequence (star), (B) depth of the moment centroid from regional moment tensor inversion, (C) example of waveform fitting from moment tensor inversion at intermediate periods (20–50 s), showing results for the January 25 mainshock. Station panels show radial, transverse and vertical displacement waveforms from top to bottom, black for observations and red for predictions (time in s, amplitudes in m). The best-fitting solution was obtained at 6 km depth. (D) Example of waveform fitting (as in C) for the M_W 5.1 reverse faulting aftershock occurred on 2016/02/22, 03:46:03 UTC, at identical centroid depth (6 km).

Our relocation yields shallow hypocentral depths as expected, mostly below 15 km. It places the 2016 mainshock at the bend of the Al-Idrisi fault, suggesting a relationship between this major earthquake and the major strike-slip fault. This agrees with most previous interpretations (Buforn et al., 2017; Medina and Cherkaoui, 2017; Kariche et al., 2018; Gràcia et al., 2019), with the exception of Galindo-Zaldivar et al. (2018), who suggest that the earthquake could be related to newly formed faults like the 1994 and 2004 Al Hoceima earthquakes. Relocation resolves four separate groups of epicenters within the seismic sequence (Figure 4). The first group, according to chronological order, is situated in the center of the 2016 earthquake series, and is formed solely by the epicenters of the January 21st foreshock, the January 25th main event and two aftershocks on January 25th and 31st (group A). This group is surrounded by a zone where no moderate aftershocks occurred. A second cluster toward the South (group B) is located ∼5 km west of the trace of the Al-Idrisi fault, and contains large part of the early aftershocks. A third group toward the East (group C) concentrates most of the later aftershocks from February to March. Finally, a fourth group consists of more scattered events that extend from the southern cluster toward the Moroccan coastline (group D). Besides the fore- and main shocks, also the southern aftershock cluster could be associated with the Al-Idrisi fault within reasonable location errors. Overall, the relocated sequence is shifted eastward by ∼5 km compared to routine catalog locations by the Instituto Geográfico Nacional (IGN). Results agree with other studies that use 3D velocity models (Buforn et al., 2017; Gràcia et al., 2019), and also with studies that keep a 1D model but incorporate additional stations in northern Morocco (Medina and Cherkaoui, 2017). A key role of azimuthal coverage for locating earthquakes beneath the Alboran Sea has been suggested previously (Santos-Bueno et al., 2019).

FIGURE 4. Regional moment tensor solutions for earthquakes in the south-central Alboran Sea during 2016 (double couple components in lower hemisphere projection). Solutions are labeled with the year, month and day of the event (Table 1) and, where is necessary, a capital letter (A–J) is indexed for events in the same day, according to temporal order. Mechanisms are offset from the epicenter position for visibility. Epicenters are obtained from relocation in a 3D Earth model, and represented in different colors according to the temporal order of events (earliest events in red, latest events in blue). The black ellipses (labeled A–D) mark four separate groups of events that are referred to in the text. Bathymetry from EMODnet Bathymetry Consortium (2018); Al-Idrisi Fault from Gràcia et al. (2019).

The first-order seismic moment tensor provides a general description of seismic wave radiation from a point source. Full waveforms are the preferred observable in moment tensor inversion, leading to a linear inversion scheme and a rather complete usage of the available information (Stich et al., 2003a; Stich et al., 2003b). We perform time-domain inversion for the deviatoric moment tensors, addressing small and moderate events of the 2016 earthquake sequence. We use three component displacement waveforms and filter in an intermediate period band from 15 to 35 s. Seismograms are rotated to a cylindrical coordinate system to approach a separation of radial and transverse components of the wavefield. Green’s functions for waveform inversion are computed with a reflectivity code for an average regional 1D Earth model (Stich et al., 2003a; Stich et al., 2003b; Stich et al., 2005; Martín et al., 2015). To incorporate the non-linear dependence of Green’s functions on source depth, we combine linear moment tensor inversion with a grid search among different trial depths and retain the best fitting combination of moment tensor and depth. We try centroid depths from 2 to 30 km, with 2 km increment. In this way, we also may assess the stability of the moment tensor solutions for different trial depths, and obtain an independent depth estimate for the seismotectonic interpretation of the sequence. The inversion involves careful manual weighting of individual waveforms to improve the overall waveform matches and explore the stability of the solutions (Figure 3). For the mainshock, an ad hoc centroid time shift of 3 s was assumed, to account for finite rupture duration, a choice that will be justified a posteriori through empirical Green function (EGF) analysis.

We were able to obtain stable moment tensor solutions for the foreshock, mainshock and 48 small and moderate aftershocks of the 2016 Alboran Sea sequence (Figure 4; Table 1). This includes the solutions that achieve a satisfactory match between observed and modeled waveforms, and do not show erratic changes for minor variations of seismogram weighting factors and centroid depths. Our solution for the mainshock shows strike-slip faulting with fault angle parameters of N39°E, 88°, −38° (strike, dip, rake) for the left-lateral nodal plane or N131°E, 51°, 183° for the right-lateral plane. The left lateral plane agrees with the kinematics and the local strike of the Al-Idrisi fault north of the fault bend, and can be likely considered the rupture plane (Buforn et al., 2017; Gràcia et al., 2019). The mainshock moment tensor shows a double-couple source with negligible compensated linear vector dipole contribution of 4%, and shallow centroid depth of 6 km. The estimated scalar moment of M₀ = 3.5 × 10¹⁸ Nm reveals the 2016 earthquake as the largest instrumental event at this part of the Africa–Eurasia plate boundary zone, larger than the 2004 Al Hoceima earthquake (M_W = 6.3, M₀ = 2.9 × 10¹⁸ Nm, using the same methodology; Stich et al., 2005) and the 1910 M_W 6.1 Adra earthquake (Stich et al., 2003a; Stich et al., 2003b). The collocated foreshock shows a very similar mechanism, with fault angle parameters of N28°E, 73°, −22° (strike, dip, rake) for the left-lateral nodal plane. Among the aftershock moment tensors, we obtain 28 solutions with magnitude M_W 4.0 and up, highlighting the importance of this sequence for the regional context. The remaining 20 moment tensor mechanisms correspond to minor aftershocks with magnitudes M_W between 3.5 and 3.9.

TABLE 1. Results from moment tensor inversion for 50 earthquakes of the 2016 sequence (compare Figure 4).

Compared to previous studies, we obtained a larger number of aftershock moment tensors [48, compared to 11 and 20 solutions, respectively, in Buforn et al. (2017) and Gràcia et al. (2019)], which is possibly a benefit from the denser station coverage in this study. Our mechanism for the mainshock is comparable to moment tensor solutions by Buforn et al. (2017) and Gràcia et al. (2019), showing similar seismic moment, although relatively minor double-couple rotation angles of 15° and 34°, respectively (Kagan, 2007), but our mechanism is different from the IGN automated moment tensor inversion showing a rotation angle of 60°. For the foreshock, we observe moderate rotations of 36° and 39° between our mechanism and solutions by IGN and Buforn et al. (2017), introduced mainly by the lower fault dip proposed in these studies for the NNE–SSW nodal plane (40°–46°). Our dip value of 73° appears more consistent with the steep dip angles of the mainshock mechanism and the Al-Idrisi fault. Instead, Gràcia et al. (2019) proposed for the foreshock a nearly opposite moment tensor mechanism (rotation angle of 93° with respect to our solution), showing right-lateral faulting along the NNE–SSW nodal plane. This orientation is not supported by the recorded waveforms, showing high similarity and identical polarity for fore- and main shocks, and also appears incompatible with the present-day stress field in the Alboran Basin and the left-lateral kinematics of the Al-Idrisi fault. From waveform and tectonic features, we think that faulting mechanisms for fore- and main shocks are likely similar. Our mechanisms are supporting a scenario of activity on one single fault segment at the beginning of the sequence, involved in the rupture of the foreshock, as well as the rupture nucleation for the mainshock four days later.

Centroid depths (Figure 3; Table 1) concentrate between 4 and 6 km and indicate shallow seismicity throughout (2–18 km), in agreement with relocated hypocenters and the shallow brittle-to-ductile transition in the area. These depth estimates from waveform inversion, although not tightly constrained, implicitly take into account the amplitude ratio between body and surface waves, and therefore are relevant support when near-source observations for shallow offshore earthquakes are lacking. Aftershock moment tensors reproduce the long duration of the 2016 Alboran Sea sequence, with 24 solutions over the first month, 18 more solutions during March, and six further moment tensor from April to June (Table 1). Our main focus in this study is the faulting style. The southern aftershock cluster (group B) is dominated by strike-slip faulting. These events tend to have different orientation compared to foreshock and mainshock; a majority of 12 solutions show strike values clustering between N5°E and N20°E for the left-lateral nodal plane, which seems more compatible with the trend of the southern Al-Idrisi fault (Figure 4). On the other hand, moment tensor inversion indicates reverse faulting for 21 aftershocks, and especially for the majority of the larger events (M_W ≥ 4.5, Table 1). Reverse aftershocks are predominant in the eastern cluster (group C, Figure 4), showing ∼NW–SE orientation of moment tensor P-axes, which rather resembles the compressional regime at the north Algerian margin or SW-Iberia (Figure 2). The cumulated seismic moment of major reverse faulting aftershocks for the 2016 earthquake is 2.3 × 10¹⁷ Nm (equivalent to a single M_W 5.5 earthquake). Overall, the combination of strike-slip faulting and reverse faulting in the 2016 earthquakes indicates a transpressional scenario, as opposed to the generally transtensional conditions inferred from previous moment tensors and the GNSS velocity field across the Alboran Sea (e.g., Stich et al., 2006; Serpelloni et al., 2007; De Vicente et al., 2008). We may characterize the 2016 sequence as a case of local slip-partitioning, where crustal deformation is taken up along subparallel strike-slip and dip-slip faults.

Finite Source Modeling for the Principal Earthquake

According to empirical source-scaling laws, a typical M_W 6.3, strike-slip earthquake has rupture length of about 22 km and rupture area of about 280 km² (Thingbaijam et al., 2017). Expected rupture areas of the largest aftershocks (M_W 5.2) are already one order of magnitude smaller. To establish the spatial relationship between the two principal clusters of aftershocks in the 2016 Alboran earthquake sequence (groups B and C in Figure 4) and the position and extension of the mainshock rupture area, we need to estimate finite source parameters for the main earthquake. In particular, we are interested in how far the rupture propagated from the nucleation zone into the restraining fault segment of the Al-Idrisi fault in one direction, and into the southern leg of the fault in the other direction.

The magnitude M_W 6.3 earthquake on January 25th is the only event of the 2016 Alboran sequence that has been recorded with sufficient quality at far distances, and teleseismic body wave recordings of this earthquake could be used to investigate the rupture process (Buforn et al., 2017; Gràcia et al., 2019). Both studies apply iterative time domain deconvolution (Kikuchi and Kanamori, 1991) to invert for the slip distribution. They use similar faulting parameters (hypocentral depth of 7 and 10 km, strike of N30°E and N214°E, and dip of 87° and 85°, respectively), and obtain overall similar slip distributions, showing maximum fault displacement of ∼1 m, total rupture length of ∼25–30 km, as well as an asymmetry of the slip distribution, with the rupture propagating predominantly northward for ∼16–20 km (Buforn et al., 2017; Gràcia et al., 2019). In both studies, slip is confined to shallow depths between 1 and 15 km. In general, teleseismic body waves are more sensitive to vertical variations of slip because these waves propagate steeply downward from the source, but the resolving power for horizontal variations may be less. To confirm the horizontal directivity of the rupture, Buforn et al. (2017) analyze spectral amplitudes of teleseismic Rayleigh waves at different azimuths, indicating a rupture propagating to N30°E azimuth, and Gràcia et al. (2019) suggest that pulse durations in regional waveforms are shorter toward north, consistent with predominant rupture propagation toward this side.

Here, we estimate apparent source time functions (ASTFs) for the main earthquake from P- and S- waveforms recorded at near-regional stations in Morocco and Spain. Source time functions are time series of moment release. The prefix “apparent” points out that these moment rate functions vary according to the direction to the recording station and the type of wave involved in the analysis, hereby encoding finite source parameters like directivity and rupture speed (e.g., Ammon et al., 2006; López-Comino et al., 2012). We attempt to isolate the ASTFs from the mainshock broadband waveforms through deconvolution of the corresponding recordings for foreshock and aftershocks. This approach assumes that waveforms from collocated small earthquakes represent EGFs (Hartzell, 1978) that contain essentially the same propagation effects as the target earthquake, but lack a significant imprint from source finiteness and rupture propagation. We perform deconvolution through spectral division, stabilized by a waterlevel at 1% of the maximum spectral amplitude, as well as a Gaussian low pass with pulse width of ∼0.75 s (at half the height of the maximum). We process time windows that contain 15 s of P-waves on radial and vertical components (except for stations EMLI and EALB, with shorter S–P times) and 10 s of S-waves on the transverse component, as well as 5 s of pre-arrival waveforms for both P- and S-waves. Several recordings saturate at the S-waves (EMLI, EALB, ANER), and amplitude information is lost. A concern with deconvolution in our case is given by the complicate and emergent waveforms recorded for shallow earthquakes at regional distances, due to multiple interactions with the crustal waveguide. We decide to enhance final ASTFs by stacking the deconvolution results for different candidate EGFs (Stich et al., 2005).

We initially include all strike-slip aftershocks (Table 1) into the analysis, without considering source proximity. The validity of EGFs and the success of individual deconvolutions are verified a posteriori according to the quality of the extracted ASTFs, validating the temporal concentration of the signal and the pre- and post signal noise level. In this step most of the deconvolved functions are discarded because of low quality, an outcome that we may relate to the characteristics of the signals, as well as the significant distances between mainshock and several trial EGF events. Used ASTFs are normalized before stacking to avoid that aftershocks with different magnitudes contribute differently to the result. We are able to retain stacked ASTFs with acceptable noise level for eight P- and eight S-waves (Figure 5). We measure apparent durations on the ASTFs from manual picks of onset and termination of the source signal. Termination is not well defined on several ASTFs where the time function fades away slowly (EADA, EBER, ALHA), in which case, we rely on nearby stations to guide the pick. We obtain apparent durations of 4.8–6.3 s for P-waves, and 4.8–8.3 s for S-waves. The longest durations for each wave type are to the southwest, and shortest durations toward north and northeast.

FIGURE 5. (A) Stacked apparent source time functions at 13 near-regional stations (blue for P-waves, red for S-waves) and apparent source durations represented in a polar plot. Duration measurements on individual apparent source time functions are annotated by vertical bars and colored segments of the waveforms. Predictions for the P and S-durations from the best fitting line source model (L1 misfit of 0.18 s) are given by blue and red curves in the polar plot. (B) Overlay of alternative solutions with acceptable L1 misfit (range of L1 between 0.18 and 0.38 s), showing the resolution of the source model (compare Table 2). (C) Best fitting line source model for a bended fault with rupture nucleation at the bend and fault azimuths in either direction incorporated as free, independent variables. The bended fault model leads to marginal improvement of the fit, and results in best azimuths of N32°E and N199°E, close to the mapped geometry of the Al-Idrisi fault bend.

We model the azimuthal variability of ASTF durations as a function of finite source parameters with a propagating line source (Cesca et al., 2011; López-Comino et al., 2012; López-Comino et al., 2016), Eq. 1:

Apparent durations depend on the angle between the fault strike φ and the azimuth of the recording station α, the rupture length L, the rupture speed v_R, the asymmetry of the rupture χ (defined between 0 for unilateral rupture and 0.5 for symmetric bilateral rupture), and the rise time t_R, which is the time that a particular point on the fault is slipping during the rupture. We solve the non-linear inverse problem through full grid search for the involved parameters in reasonable intervals and appropriate steps (Table 2). We assume local P-wave velocity of 6.1 km/s, and S-wave velocity of 3.5 km/s (Stich et al., 2005), as well as horizontal take off angles for the bulk energy of regional body waves for this shallow earthquake at regional distances. Model performance is evaluated through the average L1 misfit to P- and S-wave apparent durations at all stations.

TABLE 2. Parameters for the in line source inversion of apparent source time function durations (column 1), parameter range and step size for full grid search (column 2), best model (column 3) and range of parameters in alternative trial models that lead to tolerable misfit (column 4). Columns 4–6 summarize the setup and results for an alternative grid search that reduces trade-offs through fixed rupture speed and rise time.

The joint inversion of ASTF durations obtained from P- and S-waves yields low L1 misfit of 0.18 s (Figure 5), rupture duration of 5.2 s, rupture length of 21 km, azimuth of N24°E (corroborating left-lateral faulting) and an asymmetric pattern with 60% of the rupture propagating toward NE (Table 2). The rupture velocity is 3.0 km/s, corresponding to ∼85% of the local S-wave speed, which is a high, yet plausible value. The rise time estimate of 1.0 s agrees with the initial ascent of most seismograms and ASTFs. There are intrinsic trade-offs among rise time, length and rupture speed, therefore we use the full grid search results for an appraisal of errors. We retain alternative models that perform reasonably well applying a conservative L1 misfit threshold of 0.36 s, twice the global minimum. Variations within this set show that total duration, azimuth and asymmetry are well resolved (Table 2), while the trade-offs among rise time, length and rupture speed remain unresolved. Inversion with fixed rupture speed and rise time does lead to a well constrained length estimate of 21 ± 1 km. Finally, we generalize inversion to account for a bended rupture as suggested by the surface expression of the Al-Idrisi fault around the mainshock hypocenter. Therefore, we treat the azimuths of the long and short legs of rupture propagation as independent variables (α_S and α_L instead of α). This modification leads to a marginal improvement of global L1 fit, and estimates of N32°E for α_L and N199°E for α_S. These estimates are remarkably similar to the mapped directions of the Al-Idrisi fault on either side of the epicenter (Figure 4), supporting the hypothesis of rupture of the 2016 mainshock into the restraining fault segment.

We further attempt to use the shape of ASTFs to resolve lateral variations of fault slip. To invert for variable slip, we discretize the fault model in a regular mesh with cell dimension of 0.5 × 0.5 km², using 80 mesh elements along strike and 24 mesh elements along dip. This allows for maximum rupture length of 40 km, well above the previous estimate from apparent source durations. The hypocenter is assumed in the center of the fault model. The chosen width of the rupture zone of 12 km appears a plausible scenario for the thickness of the seismogenic layer, close to estimates for the depth of the 450°C isotherm (10 km; Soto et al., 2008), the maximum depth of faulting from Ocean Bottom Seismometer data (15 km; Grevemeyer et al., 2015), as well as including the centroid depths for more than 80% of inverted moment tensors in this study. A fault width of 12 km translates to average fault slip of 46 cm when assuming a rectangular rupture with length of 21 km, scalar seismic moment of 3.5 × 10¹⁸ Nm and rigidity of crustal rocks of 30 GPa. We test a plane fault surface (moment tensor strike of N39°E), as well as a model that incorporates a fault bend according to the mapped fault geometry, with strike of N40°E north of the hypocenter and N10°E south of the hypocenter (the hypocenter is still placed in the center of the model, and both legs are 20 km long). Synthetic ASTFs are constructed by superposition of time-lagged, triangular slip velocity functions for constant rupture speed of 3 km/s (Table 2). The approach has to deal with shortcomings of ASTFs. Besides short windows and clipping at individual stations, ASTFs show significant noise, and the complexity of several ASTFs appears mutually inconsistent, for example very abrupt terminations or the coexistence of trapezoidal and triangular shapes (Figure 5). We suspect that ASTFs can be affected by ringing artifacts and other deconvolution instabilities (López-Comino et al., 2016).

Slip inversion is solved through a search approach (López-Comino et al., 2015; López-Comino et al., 2016) based on forward modeling for pseudorandom trial distributions with prescribed statistics for the slip heterogeneity (Mai and Beroza, 2002). According to the geometry of the regional wavefield, with predominately horizontal radiation into the crustal waveguide, we cannot expect resolution for variations of slip with depth, while data are expected to be sensitive to horizontal variation. We therefore exclude slip variations with depth, and restrict heterogeneity solely to variations along fault strike (Figure 6), where we impose a 1D von Karman distribution for fault slip, with correlation length of 20 km and Hurst exponent of one, which corresponds to the theoretical value for a 2D planar slip distribution. Predictions for 1,000 trial models are compared to actual ASTFs by normalized cross correlation. Similar correlation coefficients are achieved for a plane fault model (cc = 0.80) and a bended fault model (cc = 0.81), which cannot be discriminated from regional data. Fits are poor at several stations (Figure 6), confirming the anticipated inconsistencies in the set of ASTFs. In particular, the oscillatory nature of several ASTFs cannot be reproduced, suggesting that oscillations are an artifact from limited bandwidth and an incomplete fulfillment of the underlying assumptions in EGF analysis. Nevertheless, slip distributions are consistent with bulk rupture parameters inverted from apparent durations (Table 2) and show that: 1) Slip propagates longer into the NE segment, consistent with the asymmetry obtained from apparent durations, 2) Slip in the NE segment is larger, consistent with the overall faulting mechanism (N39°E), and 3) maximum slip (∼0.6 m) occurs close the hypocenter, explaining the impulsive onsets of all ASTFs.

FIGURE 6. (A) Fault model with 40 km length and 12 km width and slip initiation in the center (star), including two planar, N10°E and N40°E striking fault segments of the same length (perspective view from SE). (B) Comparison between best possible predictions (red) and deconvolved apparent source time functions (ASTFs) (black) from finite fault inversion at constant rupture speed (3 km/s) and slip varying only in horizontal direction. Inversion yields a partial fit, suggesting that oscillations in ASTFs may be attributed to deconvolution artifacts rather than actual partitioning of the rupture into subevents. Amplitudes of ASTFs are given for moment rate in Nm/s. (C) Best fitting slip distributions from random model search (slip amplitudes in m).

In comparison with previous studies, our total rupture duration of 5.2 s from regional recordings is significantly shorter than rupture duration inferred from the teleseismic wavefield. In the later case, the estimates range from a rupture duration of 11 s from automated processing of teleseismic source time functions (Vallée et al., 2011) or 12 s from inversion of teleseismic body waves (Buforn et al., 2017) up to significantly longer source time functions with about 20 s duration for the main slip pulse (Gràcia et al., 2019). A possible starting point to understand this inconsistency is the fact that teleseismic waves are modeled with synthetic Green functions, while on regional scale the empirical waveforms are used for this purpose. Regional source durations may suffer from deconvolution artifacts, while teleseismic source durations may be influenced by secondary body wave phases generated in slow sediments and the water layer, not reproduced completely by the models. Our estimate of a high rupture speed of 3 km/s is in agreement with rupture speed of 3.0 km/s from finite fault inversion and rupture speed of 2.9–3.0 km/s from modeling teleseismic Rayleigh waves by Buforn et al. (2017), but is in disagreement with the slow rupture speed of 2 km/s proposed by Gràcia et al. (2019). For offshore earthquakes like the 2016 sequence, geodetic information, that would help to reduce the non-uniqueness of finite fault inversion, is not available, and finite source models have to be inferred from seismic recordings alone. In this context, regional source analysis from small to moderate aftershocks plays an important role. In our case, the rupture length of 21 km estimated from regional data is close to the expectation for M_W 6.3, strike-slip earthquakes (22 km according to Thingbaijam et al., 2017), while teleseismic analyses yield longer ruptures.

Discussion

In regions of slow and distributed tectonic deformation, essential tectonic information has to be inferred from small-to-moderate earthquakes. In case of the Alboran region, even the analysis of the largest earthquakes of the instrumental period so far (M_W 6.3 for the 2004 Al Hoceima and 2016 Alboran Sea earthquakes) depend on the processing of aftershocks and requires the use of regional seismograms to understand its source process, context and implications. For seismic waves, the distance range places conditions on the geometry of observations. While the teleseismic wavefield is radiated steeply downward from the source, the part of the wavefield that appears in regional seismograms is injected laterally into the crustal waveguide. As a consequence, propagation effects are different. Teleseismic wavefields are characterized by well-understood multipathing and wave conversions, leading to well-defined packages of body waves. Regional wavefields for shallow earthquakes are characterized by complex reverberations within the crustal waveguide, leading to an emergent and ringing appearance of body waves and high sensitivity to changes in hypocenter position. These characteristics do not affect significantly the performance of moment tensor inversion in an intermediate period band, but turn out to complicate the extraction of broadband ASTFs in case of the 2016 mainshock. Yet, the analysis of regional seismograms is a key to understand the earthquake. On one hand, regional data suggest a shorter rupture than previously thought, limiting the section of the Al Idrisi fault that slipped during this earthquake. On the other hand, the regional analysis of aftershocks reveals the variability of focal mechanisms during this sequence.

The magnitude M_W 6.3, January 25th, 2016 Alboran Sea earthquake nucleated at a right-stepping bend of the Al-Idrisi fault (Figure 4). This has two interesting consequences. First, the 2016 mainshock is the first major earthquake recorded in the Alboran basin that can be assigned to a major structure. The left-lateral Al-Idrisi fault is an important tectonic accident in the Alboran Sea, composed of a ∼60 km long, ∼N20°E striking northern segment, a ∼20 km long, N40°E striking central restraining bend and a southern segment striking ∼N10°E, with length of at least 20 km (Martínez García et al., 2011; Martínez García et al., 2013; Gràcia et al., 2019). At the northern termination, the Al-Idrisi fault ends at a network of secondary right-lateral faults that separate this fault from the offshore Carboneras fault (Perea et al., 2018). The southern segment reaches into an area of distributed transtensional deformation at the Moroccan margin (d’Acremont et al., 2014; Lafosse et al., 2017; Lafosse et al., 2020). No other mapped fault in the epicentral area shows appropriate strike and length to host the 2016 earthquake with a ∼20 km long and NNE–SSW oriented rupture area, according to finite fault modeling from regional data. Our moment tensor, with near-vertical dip (88°), and non-horizontal rake estimate (−38°), agrees with high-resolution seismic profiles across the Al-Idrisi fault bend that show near-vertical dip and an exhumed crystalline basement with recent (Neogene) volcanic rocks in the eastern block of the fault, while subsidence and sedimentation occurs on the western block (Martínez García et al., 2011). Second, the location of the 2016 mainshock at the bend of the Al-Idrisi fault introduces local transpression around the restraining segment of the fault. According to moment tensor inversion, transpressional deformation is partitioned into strike-slip and dip-slip displacements in different fault structures.

The 2016 earthquake is the last out of three strong earthquakes that affected the Moroccan Mediterranean margin in recent times, following the 1994, M_W 5.9 and the 2004, M_W 6.3 Al Hoceima earthquakes. Unlike the last earthquake, the 1994 and 2004 earthquakes occurred on unmapped faults that have no recognized surface expression, (e.g., Stich et al., 2005; Biggs et al., 2006; van der Woerd et al., 2014). Both events show strike-slip mechanisms similar to the 2016 earthquake, (e.g., Biggs et al., 2006), but due to a lack of surface rupture, the discrimination between left-lateral and right-lateral faulting remains debatable. In SAR images suggest rupture along a NW–SE, right-lateral fault in case of the 2004 earthquake, and a NE–SW, left-lateral fault in 1994, thus attributing the two Al Hoceima earthquakes to conjugate faulting on previously unknown structures (Akoglu et al., 2006; Biggs et al., 2006; Cakir, 2006). The 1994 and 2004 Al Hoceima earthquakes reproduce a general trend of Alboran earthquakes, where small and moderate events correlate poorly with mapped faults, and the significance of these events and faults is not clear (Grevemeyer et al., 2015; Custódio et al., 2016). Among the larger events, also the 1910, M_W 6.1 Adra earthquake has been related to a secondary, oblique-right-lateral fault within the shear zone (Stich et al., 2003b; Gràcia et al., 2012). This leaves the 2016 earthquake as the first event associated to one of the principal left-lateral faults of the Trans-Alboran shear zone, and highlights the role of the Al-Idrisi fault in accommodating regional tectonic deformation in the Alboran region (Buforn et al., 2017; Medina and Cherkaoui, 2017; Kariche et al., 2018; Gràcia et al., 2019).

At regional scale, the Trans-Alboran shear zone can be considered a transform fault, accommodating differential deformation between the Algerian and Moroccan sectors of the plate boundary, the latter affected by stretching due to slab rollback and basal drag under the Rif (Rutter et al., 2012; Spakman et al., 2018). In particular, the Al-Idrisi fault aligns with the limit of independent southwestward motion of the Gibraltar Arc with respect to Nubia documented by onshore GNSS observations (Vernant et al., 2010; Koulali et al., 2011). The importance of long term motion along the Al-Idrisi fault is evidenced from subsurface imaging (Martínez García et al., 2011; Gràcia et al., 2019), as well as possibly by the interruption and left-lateral offset of about 10 km of the Alboran Ridge, a prominent basement high with overall trend of ∼N60°E and length of ∼130 km (Figure 4; Galindo-Zaldívar et al., 2018). This amount appears consistent with a Pliocene age of the southern Al-Idrisi fault (∼3.5 Myr; Martínez-García et al., 2011; Gràcia et al., 2019) and fault slip rates of about 3 mm/a. Such slip rates coincide with GNSS-based estimates for the total left-lateral strain rates along the Trans-Alboran Shear Zone (Stich et al., 2006), indicating that the Al-Idrisi is one of the principal faults in the Alboran basin, and may possible concentrate a large part of regional strain in a single structure. This seems to contrast with the situation further south at the Moroccan margin, where the mapped fault network is becoming more distributed (d’Acremont et al., 2014; Lafosse et al., 2017; Lafosse et al., 2020), including the occurrence of strong events like the Al Hoceima earthquakes on previously unrecognized structures, which suggests the formation of new faults, and possibly a contemporary reorganization of tectonic deformation at the margin (Galindo-Zaldivar et al., 2018).

From finite source modeling, we observe asymmetric rupture propagation out from the nucleation zone at the Al-Idrisi fault bend. The rupture was propagating mainly into the restraining fault segment northeast of the hypocenter. We can compare the extend of the rupture surface with the spatial separation between the hypocenter of the main event and the locations of the two groups of aftershock toward south and northeast. In the southern aftershock cluster (group B), strike slip solutions with ∼N10°E striking left-lateral planes emulate the orientation of the southern Al-Idrisi fault. From faulting style and locations, we may suppose that the southern cluster is located on the Al-Idrisi fault itself, and likely outlines the limit of the mainshock rupture in this direction. The distance to this aftershock cluster is ∼8 km, identical to the length estimate inferred for the southern leg of the mainshock rupture. The separation of the northeastern aftershock cluster (group C) is less, and the entire cluster falls within the estimated rupture length of ∼13 km for the northeastern leg of the mainshock. In this case, aftershock activity occurs off the Al-Idrisi fault, according to the fundamentally different fault dip and faulting style. This area aligns with the fault segment inside the double restraining bend of the Al-Idrisi fault, and shows a clear expression of slip partitioning: In a scenario of local transpressional deformation, we observe the concurrence of strike-slip motion on the vertical Al-Idrisi fault during the mainshock, with dip-slip motion during the aftershock sequence. The dip-slip aftershocks occur on subparallel, low-angle south-dipping reverse faults at the western end of the Alboran Ridge. The contrast between regional transtensional deformation and local transpressional deformation can be demonstrated by a comparison of seismicity in the southern Alboran basin before and during 2016 (Figure 7). Before 2016, only three out of 23 available moment tensor solutions (13%) indicate reverse faulting earthquakes, (i.e., the T-axes are the most vertical principal axes of the moment tensor), while during the 2016 earthquake sequence, 21 out of 50 solutions (42%) were reverse faulting, and a pronounced peak appears in the rake angle histogram for directions between 60° and 120° (Figure 7).

FIGURE 7. Comparison of seismicity and focal mechanisms for the south-central Alboran Sea before (A) and during (B) the 2016 earthquake sequence. Bathymetry from EMODnet Bathymetry Consortium (2018); seismicity [red for events during 2016 in (B) yellow for other shallow earthquakes between 1990 and 2020 in (A)] according to IGN catalog (www.ign.es). For both sets of moment tensor solutions, for 23 and 50 events, respectively, we show orientations of principal compression (P-axes) and tension directions (T-axes) in polar plot, pie charts with distribution of faulting style (NF for normal faulting, i.e., P-axes is most vertical principal axes; RF for reverse faulting, i.e., T-axes is most vertical principal axes; and SS for strike-slip faulting, i.e., N-axes is most vertical principal axes) and histograms of rake angles measured indiscriminately on both nodal planes.

Slip partitioning at the restraining bend appears consistent with the shallow centroid depths from moment tensor inversion and an influence of the stress-free boundary condition at the Earth’s surface. However, an ideal realization of an Andersonian model of faulting disagrees with the outcome of a non-horizontal rake vector in the 2016 mainshock focal mechanism and the visible offset of the seafloor produced by accumulated motion on the Al-Idrisi fault (Galindo-Zaldívar et al., 2018; Gràcia et al., 2019). Otherwise, partitioning of fault motions in the brittle crust is favored by the contact with ductile rheologies at depth (Bowman et al., 2003), which is the case in the Alboran basin with shallow brittle-to-ductile transition (Fernández-Ibáñez and Soto, 2008). A more trivial way to understanding partitioning in the 2016 Alboran sequence is through the nearly vertical dip of the Al-Idrisi fault. The fault is not capable of taking up a component of horizontal shortening, and contraction produced around the restraining segment may localize on secondary faults, in this case reverse faults with low dip angles (mostly below 30°; Table 1) and the ability to take up large amounts of horizontal shortening. The 2016 sequence may have reactivated Alboran Ridge reverse faults that accommodated NW–SE shortening since the Pliocene (Martínez-García et al., 2011; Martínez-García et al., 2013) and puts on the agenda a potential role of the Alboran Ridge as source of tsunamis (Álvarez-Gómez et al., 2011). Slip partitioning during the 2016 Al Idrisi rupture provides a framework to understand the long and intense aftershock sequence of this earthquake. We finally recall other examples of local transpression along the Trans-Alboran shear zone. Slip partitioning into a reverse dip-slip component and left-lateral strike-slip has been proposed for a restraining bend in the Alhama de Murcia fault, based on geological evidence (Ferrater et al., 2017). Deep-rooted positive flower structures, cutting through the entire crust, were imaged in receiver function studies in southeastern Spain (Mancilla et al., 2018).

Finite source modeling reveals that the 2016 rupture was able to propagate at high rupture speed into both directions from the fault bend, suggesting that the northern and southern segments of the Al Idrisi fault link across a continuous fault bend at seismogenic depth, thus confirming the geological observations of a continuously curved surface trace of the fault (Martínez-García et al., 2011; Galindo-Zaldívar et al., 2018; Gràcia et al., 2019). Compared to the N10°E striking southern segment or the N20°E striking northern segment, the N40°E striking segment inside the bend of the Al-Idrisi fault is less favorably oriented with respect to the regional stress field: it forms a high angle with the maximum horizontal stress direction in the southern Alboran Sea (σ₁ at ∼N345°E; Stich et al., 2006; Fernández-Ibáñez et al., 2007; De Vicente et al., 2008; Soumaya et al., 2018). Observed slip partitioning confirms the misalignment of the restraining segment and the regional shear stress. Fault bends may be efficient obstacles for rupture propagation, and frequently form the starting or termination points for earthquakes (e.g., King and Nábělek, 1985; Elliott et al., 2015), in agreement with numerical simulations of dynamic ruptures (e.g., Nielsen and Knopoff, 1998; Aochi et al., 2002; Kase and Day, 2006). In particular, during the rupture propagation into a restraining bend, the strength excess becomes larger due to an increase of compressional normal stress. The 2016 rupture nucleated in a high-angle bend of the Al Idrisi fault and died on apparently planar domains of the fault. For uniform long-term deformation along the fault, slip along the less favorably oriented fault segments is predicted to be smaller due to slip partitioning. The opposite pattern of fault slip in the 2016 earthquake is likely the result of an inherited slip deficit. In such a model, previous large ruptures that were approaching the restraining segment could not propagate past the 30° fault bend, but have contributed to stress buildup at the nucleation zone of the 2016 rupture.

Conclusions

• We may associate the 2016 Alboran Sea earthquake with the Al-Idrisi strike-slip fault, making it the first major earthquake in the Alboran basin that can be assigned to a principal fault within the Trans-Alboran shear zone.

• In contrast to adjacent areas with distributed deformation, the Al-Idrisi fault appears capable of concentrating in a single structure large part of the ∼3 mm/yr, left-lateral deformation due to slab rollback and tearing beneath the Alboran basin.

• The M_W 6.3, 2016 Alboran Sea earthquake involved ∼21 km of left-lateral faulting, and was nucleated in the bend between the southern, N10°E and central, N40°E segments of the Al Idrisi fault.

• ASTFs suggest a fast (∼3 km/s) and asymmetric rupture, propagating primarily toward north, and showing main slip (>0.5 m) in the N40°E fault segments, consistent with the N39°E strike of the moment tensor mechanisms.

• Hence, main slip is located inside a double restraining bend, where fault strike forms a high angle with the maximum horizontal stress. Deformation along the restraining segment is partitioned between left-lateral strike-slip and reverse dip-slip faulting.

• Local transpression around the Al-Idrisi fault bend activates reverse faults at the Alboran Ridge, that make an important contribution to the intense and long-lasting aftershock activity in the 2016 Alboran Sea earthquake sequence.

• Strike-slip faulting aftershocks are observed mainly along the N10°E segments of the Al Idrisi fault, beyond the end of the 2016 mainshock rupture.

• The fore- and main shocks show nearby locations and similar mechanisms, supporting activity of the same fault segment involved in the M_W 5.1 rupture of the foreshock, and the rupture nucleation for the mainshock 4 days later.

• The 2016 earthquake appears to have settled a slip deficit in the restraining segment, suggesting that this segment may have remained unbroken in previous ruptures that could not propagate across the fault bend.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

DS and RM have realized the point source modeling, DS and JA have realized the finite source modeling. All authors have contributed to the interpretation and discussion. DS prepared the initial draft, with contributions and editing from all authors.

Funding

This study was supported by FEDER/MINECO projects CGL2015–67130-C2-2-R and PID2019–109608GB-I00, FEDER/Junta de Andalucía project A-RNM-421-UGR18, and is part of the research group RNM104 of the Junta de Andalucía. JA has also received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 754446 and UGR Research and Knowledge Transfer Found–Athenea3i; and by project 407141557 of the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

For this study we use seismic broadband data from stations operated by the Instituto Andaluz de Geofísica (IAG) in Granada, Spain, Instituto Geográfico Nacional (IGN) in Madrid, Spain, Center National pour la Recherche Scientifique et Technique (CNRST) in Rabat, Morocco, Instituto Português do Mar e da Atmosfera (IPMA) in Lisbon, Portugal, as well as the Western Mediterranean Network (WestMed) operated by the Real Observatorio de la Armada (Cádiz), Universidad Complutense de Madrid and the GeoForschungsZentrum in Potsdam, (ROA/UCM/GFZ).

References