# The Al Hoceima earthquake sequence of 1994, 2004 and 2016: Stress transfer and poroelasticity in the Rif and Alboran Sea region

The Al Hoceima earthquake sequence of 1994, 2004 and 2016: Stress transfer and poroelasticity in... Abstract The 2016 January 25 earthquake (Mw 6.3) follows in sequence from the1994 May 26 earthquake (Mw 6.0) and the 2004 February 24 earthquake (Mw 6.4) in the Rif Mountains and Alboran Sea. The earlier two seismic events which were destructive took place on inland conjugate faults, and the third event occurred on an offshore fault. These earthquake sequences occurred within a period of 22 yr at ∼25 km distance and 11–16-km depth. The three events have similar strike-slip focal mechanism solutions with NNE-SSW trending left-lateral faulting for the 1994 and 2016 events and NW-SE trending right-lateral faulting for the 2004 event. This shallow seismic sequence offers the possibility (i) to model the change in Coulomb Failure Function (ΔCFF with low μ΄ including the pore pressure change) and understand fault-rupture interaction, and (ii) to analyse the effect of pore fluid on the rupture mechanism, and infer the clock-time advance. The variation of static stress change has a direct impact on the main shock, aftershocks and related positive lobes of the 2004 earthquake rupture with a stress change increase of 0.7–1.1 bar. Similarly, the 2004 main shock and aftershocks indicate loading zones with a stress change (>0.25 bar) that includes the 2016 earthquake rupture. The tectonic loading of 19–24 nanostrain yr−1 obtained from the seismicity catalogue of Morocco is comparable to the 5.0 × 1017 N·m yr−1 seismic strain release in the Rif Mountains. The seismic sequence is apparently controlled by the poroelastic properties of the seismogenic layer that depend on the undrained and drained fluid conditions. The short interseismic period between main shocks and higher rate of aftershocks with relatively large magnitudes (4 <  Mw < 5.5) implies the pore-fluid physical effect in undrained and drained conditions. The stress-rate ranges between 461 and 582 Pa yr−1 with a ΔCFF of 0.2–1.1 bar. The computed clock-time advance reaches 239 ± 22 yr in agreement with the ∼10 yr delay between main shocks. The calculated static stress change of 0.9–1.3 bar, under pore-fluid stimulus added with well-constrained geodetic and seismic strain rates are critical for any seismic hazard assessment. Seismicity and tectonics, Dynamics: seismotectonics, Fractures, faults, and high strain deformation zones, Neotectonics INTRODUCTION Three significant earthquakes have occurred on 1994 May 26 (Mw 6.0), 2004 February 24 (Mw 6.4) and 2016 January 25 (Mw 6.5) in the Rif Mountains of Morocco and southern Alboran Sea area within a period of 22 yr (Fig. 1). The earlier two earthquakes caused severe damage due to their location inland, but the third offshore event was only felt on the nearby Moroccan coastline. Aftershock distribution (El Alami et al.1998; Bezzeghoud & Buforn 1999) and surface deformation as deduced from Synthetic Aperture Radar Interferogram (InSAR) (Akoglu et al.2006; Cakir et al.2006; Tahayt et al.2009), indicate that the 1994 and 2004 events occurred on NNE–SSW and NW–SE trending conjugate strike-slip faults, respectively. The 2016 event located about 20 km offshore is associated with a NNE–SSW trending rupture with a similar mechanism to the 1994 event. This seismic sequence is unusual in the North Africa active zones because earthquake ruptures are within ∼25 km area and the time interval between main shocks is about 10–12 yr. Figure 1. View largeDownload slide Seismicity of the Al Hoceima region showing the seismic sequence of 1994 May 26 (El Alami et al.1998), 2004 February 24 (Tahayt et al.2009) and 2016 January 25 (CSEM, http://www.emsc-csem.org/#2). The 1994 aftershocks are in green, the 2004 aftershocks are in blue and 2016 aftershocks are in grey. Focal mechanisms are Harvard–CMT (Table 1). Inset represents the plate boundary in the Alboran Sea with convergence rate in mm yr−1 (Koulali et al.2011). Figure 1. View largeDownload slide Seismicity of the Al Hoceima region showing the seismic sequence of 1994 May 26 (El Alami et al.1998), 2004 February 24 (Tahayt et al.2009) and 2016 January 25 (CSEM, http://www.emsc-csem.org/#2). The 1994 aftershocks are in green, the 2004 aftershocks are in blue and 2016 aftershocks are in grey. Focal mechanisms are Harvard–CMT (Table 1). Inset represents the plate boundary in the Alboran Sea with convergence rate in mm yr−1 (Koulali et al.2011). In this paper, we first present the seismic sequence and suggest a fault-rupture interaction using Coulomb modeling on fixed planes. Secondly, the computation using optimally oriented planes is added to constrain aftershocks distribution. The three main shocks and related aftershocks appear to be closely related and their location implies a stress transfer with triggering. We show that the modeled stress distribution and seismicity rate change suggest a pore-fluid effect correlated with elastic dislocation in undrained and drained conditions. The Coulomb Failure Function change (ΔCFF) and pore-fluid flow seem to control the 10–12 yr recurrence of main seismic events with a clock advance. The plate boundary tectonic condition, the seismicity rate change and poroelastic properties of the seismogenic crust seem to play a significant role in the triggering of earthquakes in the Rif Mountains and Alboran Sea. SEISMOTECTONIC SETTING The seismicity of the Rif Mountains and Alboran Sea is due to the convergence between Africa and Eurasia (Iberia) in the western Mediterranean. The E–W trending Rif Mountains run along the northern coast of Morocco forming the southern branch of the Betic-Rif arc that includes the Alboran Sea and belongs to the transpression plate boundary system in the Western Mediterranean region (Morel & Meghraoui 1996; Meghraoui & Pondrelli 2012). Tahayt et al. (2009) interpret the region as a transrotational regime applied to the Oriental Rif block with a clockwise rotation. This complex tectonic domain also results from a Neogene subdued subduction zone with lithospheric delamination where the Alboran Sea appears as an oceanic microplate (Calvert et al.1997). As indicated by the three main shocks (Fig. 1 and Table 1), the present-day tectonic framework of the Al Hoceima region is dominated by a strike-slip fault regime where moment magnitudes do not exceed 6.5. Limited fault-rupture dimensions are likely due to the local structural geology made of overthrusting nappes on highly deformed continental crustal rocks (Chalouan et al.2008; Timoulali et al. 2014). Because of the limited number of local seismic stations in the Rif, the location of aftershocks of the 1994 and offshore 2016 earthquakes are poorly resolved, which is not the case for the 2004 earthquake (Tahayt et al.2009). Table 1. Physical characteristics of the earthquakes used in this study, HRV designed Harvard solution and CMT designed Centroid Moment Tensor solution (see also Fig. 1). Earthquake  Long.  Lat.  Mo (1018 nm)  Mw  U (m)  L (km)  W (km)  Strike  Dip  Rake  1994 May 26  −3.99  35.28  1.01  6.1  0.8  16  10  17  85  −7  2004 February 24  −3.99  35.142  3.0  6.4  1.0  19  16.5  340  87  −161  2016 January 25  −3.70  35.67  4.69  6.5  0.8  25  13.5  214  78  19  Earthquake  Long.  Lat.  Mo (1018 nm)  Mw  U (m)  L (km)  W (km)  Strike  Dip  Rake  1994 May 26  −3.99  35.28  1.01  6.1  0.8  16  10  17  85  −7  2004 February 24  −3.99  35.142  3.0  6.4  1.0  19  16.5  340  87  −161  2016 January 25  −3.70  35.67  4.69  6.5  0.8  25  13.5  214  78  19  View Large Several authors have studied the seismicity and suggested moment tensor solutions of major earthquakes showing a NNW–SSE contraction stress regime in a predominantly strike-slip faulting domain associated with normal and thrust mechanisms (Hatzfeld et al.1977; Cherkaoui et al. 1990; Medina 1995; Stich et al.2006, 2010; Palano et al.2013). Although the seismicity may appear diffuse in the Rif-Alboran Sea, the 1994, 2004 and 2016 seismic events reveal a clear migration of earthquake ruptures (Fig. 1). The seismicity and tectonics of the region indicate a clear correlation between continental and offshore faults from both the Betics and Rif Mountains towards the Alboran Sea (Grevemeyer et al.2015). From detailed bathymetry and seismicity distribution, Gevemeyer et al. (2015) identify a fault zone crossing the Alboran ridge at the location of the 2016 January 25 earthquake (Mw 6.5). Using InSAR, Cakir et al. (2006) and Akoglu et al. (2006) constrain the coseismic earthquake surface deformation and provide the rupture parameters of the N23° E trending 1994 and N45° W trending 2004 earthquakes using elastic dislocations. They suggest that the two main shocks occurred on blind and conjugate strike-slip faults, with left- and right-lateral slip, respectively. The time-series analysis of SAR data (Envisat) for 7 yr following the 2004 earthquake shows that post-seismic deformation reaches up to 4 cm at the surface and infers 0.3 m displacement at shallow depth (<7 km), mainly above the high coseismic slip patches which can be explained by 5.0 × 1017 N·m (Mw 5.73) cumulative moment release (Cetin 2015). The 2016 earthquake rupture, located about 20 km further north, is modeled for the source time function to obtain a strike-slip faulting mechanism with a NNE–SSW main rupture in agreement with aftershock distribution (Vallée 2016). The three focal mechanism solutions and the 1994 and 2004 rupture geometries inferred from surface deformation are used as an input for the ΔCFF modeling (Table 1). MODELING AL HOCEIMA SEQUENCE BY COULOMB FAILURE FUNCTION The Rif region of northern Morocco experienced a seismic sequence with three moderate to large earthquakes within 22 yr. The sequence suggests earthquake triggering, fault interaction and stress transfer as observed in other earthquake areas (Hudnut et al.1989; Stein et al.1997). The case study uses the applied stress change calculated as ΔCFF (Reasenberg & Simpson 1992; King et al.1994) expressed by:   $$\Delta {\rm{CFF}}\ = \ \Delta \tau \ -\ \mu \ \left( {\ \Delta {\sigma _n} - \Delta P} \right)$$ (1)  $$\Delta {\rm{CFF}}\ = \Delta \tau \ - \ \mu '\Delta {\sigma _n}$$ (2)where τ is the shear stress, σn is the normal stress (compression positive), P is the pore-fluid pressure, μ and μ΄ are the coefficient of friction and effective coefficient of friction, respectively and Δ refers to changes during the earthquake. The apparent friction is given by (Reasenberg & Simpson 1992)   $$\mu ' = \ \mu \ \left( {1\ - B} \right)$$ (3)where B is the Skempton coefficient which defines the relation between the stress change and pore pressure change (Beeler et al.2000)   $$B\ = \ \ \frac{{\Delta P}}{{\Delta {\sigma _m}}} = \ 3\ \frac{{\Delta P}}{{\Delta {\sigma _{kk}}}}$$ (4)where Δσm is the mean stress change and Δσkk is the sum over the diagonal elements of the stress tensor. It is important to note that for an isotropic model, the apparent friction coefficient used in triggered seismicity is defined by the combination of pore pressure and friction coefficient:   $$\mu ' = \ \mu \ \left( {1 - \frac{{{\rm{\Delta }}P}}{{{\rm{\Delta }}{\sigma _n}}}} \right)$$ (5)Substituting eqs (5) and (4) in eq. (2), we obtain   $${\rm{\Delta CFF}}\ = \ {\rm{\Delta }}\tau + \ \mu \ \left( {{\rm{\Delta }}{\sigma _n} - B{\rm{\Delta }}{\sigma _m}} \right)$$ (6a) Beeler et al. (2000) suggest that due to the pore-fluid effect and for an isotropic poroelastic model, the two expressions defined in eqs (2) and (6) yield different results in some modeling configurations. The variable effective friction coefficient related to the variation of the Skempton coefficient B and the pore pressure change along the fault zone gives more realistic solutions especially at high pore pressure change than the imposed constant effective friction commonly used in the Coulomb stress modeling (see Appendices A and B for more details). The variation of B must be considered when different porosity and different diffusive processes are present in the fault zone (Scholz 1990). From the structural point of view, fluid migration in the host rock may occur on nearby subsidiary fractures linked to bounding faults directly related to earthquake and aftershock behaviour (Kirkpatrick et al.2008). The cumulative moment for Al Hoceima region (inland and offshore) reaches 1.1 × 1019 N·m, since 1994, and the seismic strain release tested on a fault network reveals a low effective coefficient of friction on the fault (≤0.2) to model the slip rates (Negredo et al.2002). In order to model the active tectonics of the Ibero-Maghrebian region, a strain rate of 15–40 nanostrain yr−1 is obtained using slip rates on faults and Global Positioning System (GPS) data (Negredo et al.2002; Koulali et al.2011; Palano et al.2013). In their calculation of frictional strength, Negredo et al. (2002) assume no cohesion in the media. In our work, the modeling is performed using Coulomb 3.4 software (Toda et al.2011) based on the conversion of DC3D subroutines (Okada 1992) to calculate the ΔCFF. The static stress change is computed on both fixed receiver earthquake ruptures and optimally oriented faults using an effective coefficient of friction μ΄ = 0.4, and fault parameters summarized in Table 2. The computed stress changes are presented in Table 3, where the value of static stress change according to the receiver fault represents the ΔCFFmax and also expresses the Coulomb stress drop on the source faults. Table 2. Rupture parameters of significant earthquakes from InSAR results (Akoglu et al.2006) and source time function (http://geoscope.ipgp.fr/index.php/en/catalog/earthquake-description?seis = us10004gy9) used for the Coulomb stress transfer modeling. Long_c and Lat_c indicate the centre of each dislocation. EQ  Long_c (°)  Lat_c (°)  Top-depth (km)  Bot-depth (km)  L (km)  W (km)  Strike (°)  Dip (°)  Rev.slip (m)  Right_lat. slip (m)  M0 (dyne·cm) ×1025  1994  −4.01  35.17  1.00  11.00  16.00  10.04  17  85  −0.04  −0.32  1.7  2004  −3.97  35.13  1.00  16.00  20.63  15.02  340  87  −0.22  0.63  6.6  2016  −3.84  35.5  1.00  14.47  28.70  13.50  205  86  0.2  −0.57  7.5  EQ  Long_c (°)  Lat_c (°)  Top-depth (km)  Bot-depth (km)  L (km)  W (km)  Strike (°)  Dip (°)  Rev.slip (m)  Right_lat. slip (m)  M0 (dyne·cm) ×1025  1994  −4.01  35.17  1.00  11.00  16.00  10.04  17  85  −0.04  −0.32  1.7  2004  −3.97  35.13  1.00  16.00  20.63  15.02  340  87  −0.22  0.63  6.6  2016  −3.84  35.5  1.00  14.47  28.70  13.50  205  86  0.2  −0.57  7.5  View Large Table 3. Shear, normal and Coulomb stress change for the major earthquakes. Fault geometries used in the Coulomb stress modeling are defined in Table 2. SF and RF are the source and receiver faults, respectively. EQ (SF-RF)  Calc_Location  Receiver faults (°)  Stress computation (bar)    Long (°)  Lat (°)  Z (km)  Strike  Dip  Rake  Shear  Normal  Coulomb  1994(SF)  −4.010  35.166  7.00  340  87  −161  1.997  −8.547  −1.422  2004(RF)  −3.931  35.088  7.00        0.656  1.104  1.062  1994(SF)  −3.987  35.131  7.00  195  78  19  −8.146  −1.227  −8.637  2016(RF)  −3.870  35.429  7.00        0.093  0.048  0.112  2004(SF)  −4.008  35.166  7.00  195  78  19  −4.594  −11.509  −9.198  2016(RF)  −3.870  35.429  7.00        0.228  0.033  0.241  EQ (SF-RF)  Calc_Location  Receiver faults (°)  Stress computation (bar)    Long (°)  Lat (°)  Z (km)  Strike  Dip  Rake  Shear  Normal  Coulomb  1994(SF)  −4.010  35.166  7.00  340  87  −161  1.997  −8.547  −1.422  2004(RF)  −3.931  35.088  7.00        0.656  1.104  1.062  1994(SF)  −3.987  35.131  7.00  195  78  19  −8.146  −1.227  −8.637  2016(RF)  −3.870  35.429  7.00        0.093  0.048  0.112  2004(SF)  −4.008  35.166  7.00  195  78  19  −4.594  −11.509  −9.198  2016(RF)  −3.870  35.429  7.00        0.228  0.033  0.241  View Large In the ΔCFF modeling computed with fixed strike-slip receiver fault plane, the correlation between loading lobes and aftershock distribution suggests a close interaction between the 1994, 2004 and 2016 earthquake ruptures (Figs 2a–c, 3a and Table 3, and Figs S1a and b, S2 and Table S1, Supporting Information). Unlike the 2016 earthquake epicentre which is clearly within positive CFF lobes, the 2004 earthquake location given by the InSAR analysis (Akoglu et al.2006) is at the transition from negative to positive ΔCFF lobes. Taking the 1994 as a source fault and the 2004 rupture as a receiver fault, and using a low friction coefficient (μ΄ < 0.4 and for isotropic models, Figs 3a and b), our CFF modeling show the 2004 fault rupture clearly located in a loading zone with positive stress change. Figure 2. View largeDownload slide (a) Calculated ΔCFF with the 1994 source fault, and related aftershock distribution (see the text for explanation) and the 2004 as a receiver fault (strike/dip/rake = 340°/87°/−161°), the blue, white and black stars are epicentres of the 1999, 2004 and 2016, respectively (same symbols in figures b–d). (b) Computed ΔCFF with the 2004 as a source fault and the 2016 as a receiver fault (strike/dip/rake = 195°/78°/19°), and related aftershock distribution (see the text for explanation). (c) Computed ΔCFF with the 2016 source fault for fixed planes (strike/dip/rake = 195°/78°/19°), and related aftershock distribution (see the text for explanation). (d) Computed cumulative ΔCFF with the three source faults on a fixed planes (strike/dip/rake = 195°/78°/19°), and related aftershock distribution (see the text for explanation). The 1994 and 2004 mapped fault ruptures are from Akoglu et al. (2006), the 2016 fault-rupture model is from M. Vallée (http://geoscope.ipgp.fr/index.php/en/catalog/earthquake-description?seis=us10004gy9). Figure 2. View largeDownload slide (a) Calculated ΔCFF with the 1994 source fault, and related aftershock distribution (see the text for explanation) and the 2004 as a receiver fault (strike/dip/rake = 340°/87°/−161°), the blue, white and black stars are epicentres of the 1999, 2004 and 2016, respectively (same symbols in figures b–d). (b) Computed ΔCFF with the 2004 as a source fault and the 2016 as a receiver fault (strike/dip/rake = 195°/78°/19°), and related aftershock distribution (see the text for explanation). (c) Computed ΔCFF with the 2016 source fault for fixed planes (strike/dip/rake = 195°/78°/19°), and related aftershock distribution (see the text for explanation). (d) Computed cumulative ΔCFF with the three source faults on a fixed planes (strike/dip/rake = 195°/78°/19°), and related aftershock distribution (see the text for explanation). The 1994 and 2004 mapped fault ruptures are from Akoglu et al. (2006), the 2016 fault-rupture model is from M. Vallée (http://geoscope.ipgp.fr/index.php/en/catalog/earthquake-description?seis=us10004gy9). Figure 3. View largeDownload slide (a) ΔCFF for various effective friction coefficients (μ΄) along the 2004 rupture (strike direction) as receiver fault with strike/dip/rake = 340°/87°/−161° at 7 km depth. The profiles start at the intersection of the two cross faults and terminate at the end of the aftershock sequence (SE of epicentre area). The increases of pore pressure could trigger events in regions where reduced Coulomb stress predict and absence of activity. We show that at the epicentre area, the rupture nucleation occurs when pore fluids are redistributed, the value of μ΄ ≤ 0.4 seems to be more adaptable for a chosen receiver fault geometry. The computation is based on the effective constant friction model. (b) Stress change caused by the 1994 earthquake (source fault) along the 2004 fault zone due to (1) coseismic stress change due to elastic dislocation (red line), and (2) stress change due to the coupled poroelastic effect (green line). The maximum stress load on the 2004 fault zone are given for an isotropic model when μ΄ < 0.4 that implies Skempton coefficient B = 0.9 near the 1994 rupture and for μ΄ = 0.4 with B = 0.47 far from the 1994 rupture. Figure 3. View largeDownload slide (a) ΔCFF for various effective friction coefficients (μ΄) along the 2004 rupture (strike direction) as receiver fault with strike/dip/rake = 340°/87°/−161° at 7 km depth. The profiles start at the intersection of the two cross faults and terminate at the end of the aftershock sequence (SE of epicentre area). The increases of pore pressure could trigger events in regions where reduced Coulomb stress predict and absence of activity. We show that at the epicentre area, the rupture nucleation occurs when pore fluids are redistributed, the value of μ΄ ≤ 0.4 seems to be more adaptable for a chosen receiver fault geometry. The computation is based on the effective constant friction model. (b) Stress change caused by the 1994 earthquake (source fault) along the 2004 fault zone due to (1) coseismic stress change due to elastic dislocation (red line), and (2) stress change due to the coupled poroelastic effect (green line). The maximum stress load on the 2004 fault zone are given for an isotropic model when μ΄ < 0.4 that implies Skempton coefficient B = 0.9 near the 1994 rupture and for μ΄ = 0.4 with B = 0.47 far from the 1994 rupture. The cumulative ΔCFF from 1994 to 2016 (Fig. 2d) explains the present day and most important part of the seismicity and aftershock distribution and Coulomb stress drop in the Al Hoceima region. In order to analyse, the relationship between the main shocks and aftershocks distribution in detail, we model a static stress change caused by the 2004 and the 2016 earthquakes on some major aftershocks (see Table S1, Supporting Information). The modeling on fixed receiver planes suggest a stress load range between 0.5 and 0.8 bar for 2004 aftershock sequences with an optimal strike and dip value ranging between 207° and 298° and 66° and 84°, respectively. For the 2016 earthquake, major aftershocks are reached by a positive ΔCFF ranging between 0.1 and 0.4 bar, with an optimal value of [250°–260°] for strike and [40°–45°] for dip. In order to take into account all aftershock sequences, we compute a static stress change on optimally oriented fault planes as this approach does not need to include the focal mechanism of each rupture. As obtained by previous works on the northern Morocco tectonics (Medina 1995; Akoglu et al.2006; Fernandez-Ibanez et al.2007), a regional stress field (with a NW–SE principal stress direction as σ1) is added as a pre-existing stress field on the stress modeling. Based on the seismic tensor inversion and GPS data, the inferred stress field and maximum horizontal stress for Al-Hoceima–Alboran region is in good agreement with convergence models along the plate boundary (Demets et al.2010; Meghraoui & Pondrelli 2012). Note that the magnitude of the principal pre-existing stresses does not change the static coulomb stress modeling, because the stress levels are largely dominated by the coseismic rupture process in the near field. At seismogenic depth and for optimal failure planes, the static stress change modeling due to the 2004 earthquake on optimally oriented fault planes suggest that ∼30 per cent of aftershocks hypocentres were pushed closer to failure for high effective friction coefficients (Figs S2c and d, Supporting Information), while this percentage rises to ∼90 per cent when pore fluid are redistributed (Figs S2a and b, Supporting Information). The elastic modeling in Fig. 3(a) represents the stress change caused by the 1994 earthquake along the 2004 fault zone (with strike-slip mechanism, see also Table 3). The ΔCFF profiles are constrained by the aftershock distributions and the computation is performed for receiver faults with strike/dip/rake = 340°/87°/−161°. For the modeling procedure, we assume a 0.25 typical value of Poisson ratio with 8 × 105 bar for the Young modulus and 3.3 × 105 bar for the shear modulus in the seismogenic layer (5–15 km thickness). Here, the ΔCFF modeling require a low effective friction coefficient (μ΄ ≤ 0.4) denoting a pore-fluid effect. At the 2004 epicentre area, the modeling suggest that the rupture nucleation occurs when pore fluids are redistributed, with the μ΄ ≤ 0.4 considered as an optimal value also explains the nucleation process for a chosen receiver fault geometry. The increase of pore pressure may trigger seismic events in regions where reduced Coulomb stress due to the high effective friction coefficient predict and absence of activity (see also Fig. S2, Supporting Information). In fact, the low friction coefficient implies an optimum value of stress loading where the pore-fluid component takes an important role in stress transfer and earthquake triggering. The earthquake triggering caused by the 1994 seismic event results from 0.1 to 1.0 bar ΔCFF at the 2016 and 2004 receiver faults, respectively (Fig. 3a and Table 3). The ΔCFF modeling shows that the stress transfer due to the 1994 earthquake promotes the 2004 earthquake failure (Figs 2a and 3 and Table 3, and Fig. S1a, Supporting Information), and both 1994 and 2004 earthquakes promote the failure of the 2016 earthquake (Figs 2a and b and Table 3, and Fig. S1b, Supporting Information). The cumulative post-seismic deformation due to the elastic dislocation increases the stress loading on the 2016 rupture from 0.24 to 0.3 bars. ROLE OF PORE FLUID IN THE EARTHQUAKE SEQUENCE A stress change may result from pore-fluid diffusion. If the stress field satisfies the strain compatibility equation $$\frac{{{\partial ^2}}}{{\partial x_j^2}}[ {\frac{{2( {{\nu _u} - \nu } )}}{{B( {1 - \nu } )( {1 + {\nu _u}} )}}\ P + \sigma } ]$$ = 0 (Rice & Cleary 1976; see also appendix A for more details) and if we consider the boundary condition (pore pressure is neglected far from the fault), a simple solution is given by Bosl & Nur (2002):   $$\frac{\sigma }{{{\sigma _{{\rm{init}}}}}} = \ \frac{{\left( {{\nu _u} - \nu } \right)}}{{\left( {1 - \nu } \right)\left( {1 + {\nu _u}} \right)}}$$ (6b)where σ is the change in the stress field due to the pore-fluid diffusion, σinit is the initial stress induced by the coseismic dislocation and νu and ν are the undrained and drained Poisson ratios, respectively. This relation also shows the correspondence between the post-seismic mean stress change induced by pore pressure relaxation and the mean stress caused by the initial dislocation. Taking into account the Rice & Cleary (1976) and Bosl & Nur (2002) solutions, the short-term poroelastic deformation is defined as a diffusive process and can be interpreted as a linear combination of pore pressure and mean stress changes. The drained and undrained Poisson ratios used in the coupled poroelastic stress modeling are 0.25 and 0.31, respectively (Fig. 3b); these values are typical in unconsolidated sedimentary aquifers and water saturated rocks in the upper few kilometres of the seismogenic zone (h ≤ 15-km depth). Our coupled poroelastic modeling suggest a value of short-term post-seismic stress equal to 0.6 bar due to the coupled poroelastic effect for μ΄ = 0.4 related to a an internal friction coefficient of μ = 0.75 and Skempton coefficient B = 0.47. The μ΄ < 0.4 effective friction coefficient gives a value of μ = 0.75 and Skempton coefficient B = 0.9. The variation of B at the intersection of the two ruptures to the end of the 2004 rupture zone can be interpreted as a variation of porosity and the diffusive process along the fault zone (Scholz 1990). The same phenomenon is observed in reservoir induced seismicity where the elastic and the coupled poroelastic effects are considered instantaneous (Scholz 1990). The diffusive process from the reservoir is associated with the fluid migration at short time delay, and moving from a region with a high B to a region with a low B, indicating a short-time poroelastic rebound where the fluid flow transfer from the 1994 and 2004 rupture zones can be compared as reservoir induced seismic activity. Taking into account the time of occurrence of aftershocks and the short-term post-seismic stress change induced by the main shock using the Bosl & Nur hypothesis (2002), we note that the stress load occurring 10–100 d after the 2004 main shock is far more important than the stress load due to the several years of post-seismic deformation computed by InSAR time-series (PS and SBAS, Cetin 2015). The comparable computations take into account the 2004 rupture as source fault and 2016 rupture as receiver fault, and the post-seismic deformation added as a cumulative moment is incorporated into the elastic dislocation modeling. Due to the absence of pore pressure in-situ data and difficulty to perform a 3-D model of stress change and related poroelastic dislocation, we use eqs (2) and (6) to evaluate the stress and pore pressure changes related to the Al Hoceima earthquake sequence. The apparent constant friction model (eq. 2) and variable (or isotropic) friction model (eq. 6 and in appendix B) are considered as parameters able to improve our knowledge on the pore-fluid effect in the Al Hoceima-Alboran region. Nevertheless, the isotropic poroelastic model appears to be the most appropriate for modeling the short-term diffusive process in a complex fault zone. The low value of μ΄ (≤0.4) at the intersection between the 1994 and the 2004 fault (Fig. 4) implies that the first shock reduces the pore pressure along the 1994 rupture, while it increases along the 2004 rupture. The low value of μ΄ implies a high values of stress change due to the short-term coupled poroelastic effect (Fig. 3b). Fig. 4 shows how the pore pressure can affect the triggered seismicity when hydrological processes are coupled to the rupture process. The large pore pressure change associated with the volumetric strain corresponds to fluid migration close to and from the 1994 rupture to the 2004 rupture zone (Fig. 4). Das & Scholz (1981) suggest that the stress effect on the fault might be enhanced as P (see eq. 1) and restored on the main fault, causing fluids to migrate into the receiver fault. The same observation is also made by Jonsson et al. (2003) for the Mw 6.5 2000 June earthquake of SW Iceland showing an increasing water level change associated with high pore pressure change in unconsolidated sedimentary aquifers. The evaluation of pore pressure change at depth in Fig. 4 is resolved for a homogeneous elastic half-space using the theory of linear elasticity (Rice & Cleary 1976); eq. (5) (see also Appendix A5) is used to evaluate the pore pressure change by computing the mean stress change due to the 1994 earthquake for a constant Skempton coefficient. Figure 4. View largeDownload slide Calculated pore-pressure change based on the coseismic volumetric strain and theory of linear poroelasticity (Rice & Cleary 1976) following the 1994 earthquake. The 1994 and the 2004 rupture are represented by black lines (see also Fig. 2d). The 1994 coseismic slip create a high pore pressure zone in the rupture nucleation zone of the next 2004 earthquake, according to Terzaghi (1925) definition of effective stress, the increase in pore pressure diminishes the normal stress acting on the fault and promote the 2004 failure. Figure 4. View largeDownload slide Calculated pore-pressure change based on the coseismic volumetric strain and theory of linear poroelasticity (Rice & Cleary 1976) following the 1994 earthquake. The 1994 and the 2004 rupture are represented by black lines (see also Fig. 2d). The 1994 coseismic slip create a high pore pressure zone in the rupture nucleation zone of the next 2004 earthquake, according to Terzaghi (1925) definition of effective stress, the increase in pore pressure diminishes the normal stress acting on the fault and promote the 2004 failure. The successive earthquakes in the Al Hoceima region may be a response to stress loading in the Rif Mountains and related pore-fluid diffusion in the upper crust. The coupling between crustal deformation and pore-fluid effect imply that the pore pressure may decrease the rock strength by reducing the effective stress and creating a slip-instability that favours earthquake triggering. The seismicity rate change following the 2004 earthquake shows additional aftershocks (see relative rate fluctuations in Fig. 5a), which cannot be explained as aftershock rate decrease as predicted by the Omori law. As the afterslip along seismic ruptures and related post-seismic deformation (Cetin 2015) may have enhanced the pore-fluid flow and contributed to the earthquake triggering in the Rif-Alboran region, in our case, we include a coupled poroelastic component in the ΔCFF modeling and analyse the post-seismic stress transfer and related effective friction coefficient (μ΄ ≤ 0.4). The ΔCFF due to the full poroelastic relaxation is simply computed from a dislocation in a homogeneous elastic half-space using both undrained and drained Poisson's ratios and obtain the difference as seen in Fig. S3 in the Supporting Information. Figure 5. View largeDownload slide (a) The seismic frequency and relative seismicity rate following the 2004 main shock (blue line) and cumulative number of seismic events (green line). The fluctuation in the seismicity rate change shows additional aftershocks possibly due to pore-fluid diffusion in the upper crust. The relative rate change are obtained from changes in slop of the cumulative number curve using a Habermann function regardless to the time of greatest change and comparing the rate in the two parts of the period (before and after the division point) by appropriate time windows function (Wyss & Habermann 1988; Wyss & Viemer 2000), the time variation function defines the local time variation between the rate before and after. (b) The seismicity rate change versus time in the Al Hoceima region. We show the complexity of aftershocks sequences as a realistic representation of the temporal post-seismic effect. Figure 5. View largeDownload slide (a) The seismic frequency and relative seismicity rate following the 2004 main shock (blue line) and cumulative number of seismic events (green line). The fluctuation in the seismicity rate change shows additional aftershocks possibly due to pore-fluid diffusion in the upper crust. The relative rate change are obtained from changes in slop of the cumulative number curve using a Habermann function regardless to the time of greatest change and comparing the rate in the two parts of the period (before and after the division point) by appropriate time windows function (Wyss & Habermann 1988; Wyss & Viemer 2000), the time variation function defines the local time variation between the rate before and after. (b) The seismicity rate change versus time in the Al Hoceima region. We show the complexity of aftershocks sequences as a realistic representation of the temporal post-seismic effect. The pore-fluid flow in drained condition should reduce the normal stress σn (eq. 2) and favour significant additional aftershocks able to affect the relative seismicity rate change (Nur & Booker 1972; Cocco & Rice 2003). In our case, the deep-seated water flow in the flysch units of the fold-and-thrust Rif Belt and substratum metamorphic complex (Chalouan et al.2008) can be considered as another mechanism responsible for the seismicity increase. Piombo et al. (2005) suggest that for significant earthquakes with Mw  > 6, the stress transfer may occur under a fluid diffusion within a 10–20 km radius. The temporal evolution of seismicity within individual fractures takes 10–12 yr to travel up to 20 km (5–6 km yr−1) in our case as well as in other case studies (e.g. Pytharouli et al.2011). Taking into account the complexity of aftershock sequences, a realistic representation of the temporal post-seismic factor, the modified Omori Law (Utsu 1969) or Omori-Utsu Law (Narteau et al.2009) can be expressed as:   $$\lambda \left( t \right) = \ \frac{k}{{{{\left( {t + c} \right)}^{ - x}}}}$$ (7)where λ is the aftershock frequency within a given magnitude range, t is the time from the main shock triggered event, k is the productivity of aftershocks that depends on the total number of events, x is the power law exponent and c is the time delay before the onset of the power-law aftershock decay rate and is dependent on the rate of activity in the earlier part of the seismic sequence. The change of c values characterizes the aftershocks sequence and can be correlated with the stress field orientation (Narteau et al.2009), Guo & Ogata (1997) obtain a range of c value between 0.003 and 0.3 d for various earthquake datasets. In our case, the c value has to be the lowest possible and is fixed as 0.01 d, in order to obtain sufficient aftershock productivity. Note that the c value is often retained connected to the incompleteness of seismic catalogues soon after strong earthquakes. For aftershocks, the seismicity decay requires a time-dependent process that is much faster than the large-scale tectonic loading and much slower than the propagation of elastic waves (Nur & Booker 1972). In our application, we observe that immediately after the 2004 main shock, the aftershock rate decays by $$\frac{1}{{\sqrt t }}$$ while it becomes equal to $$\frac{1}{t}$$ in subsequent months. It appears that the decay rate is due to fluid flow in the crust by means of a diffusion process that contributes to the aftershock sequence (Fig. 5b). Similar results are obtained by Shapiro et al. (1997) for the events related to pressure changes in operation wells and recently by Turkaya et al. (2015) in laboratory experiments. THE CLOCK-TIME ADVANCE AND PERIODIC FREQUENCY The Al Hoceima earthquake sequence shows about 10 and 12 yr recurrence with an aftershock distribution and stress loading that correlate with the location of earthquake ruptures (Figs 2d and Figs S2a and S3, Supporting Information). The mechanism controlling the dependence time of earthquake ruptures, aftershocks and related stress change is complex. The well-resolved 2004 aftershocks distribution during a short time interval (100 d; Tahayt et al.2009) confirms our observation that a significant seismicity rate change is observed 50 d after the main shocks (see Fig. 5a); the positive change in seismicity is correlated with the positive aftershocks productivity due to pore-fluid diffusivity (Fig. 5b). To study the influence of the coseismic stress change on recurrence time interval, Chéry et al. (2001) point that a positive shear stress change on a fault plane should advance the time of the next earthquake on this fault. Here, we consider that the nucleation will occur where the ΔCFF has a maximum value (Console et al.2010). The clock-time advance and related recurrence time (in years) depends on the stressing rate in the positive ΔCFF expressed by the linear equation (Stein et al.1997):   $${\rm{Tr'}} = {\rm{Tr}}-\Delta {\rm{CFF}}/\dot{\tau },$$ (8)where Tr΄ is the calculated recurrence time, $$\Delta t = \Delta {\rm{CFF}}/\dot{\tau }$$ is the clock-time advance ($$\dot{\tau }$$ is stressing rate) and Tr is the mean recurrence time before the earthquake. The stressing rate is computed from the strain rate which is derived from the seismic moment following the equation: $$\dot{M} = 2\mu \Sigma W\dot{\varepsilon }/k$$, where μ is the shear modulus, Σ is the surface area of the region, W is the seismogenic thickness, $$\dot{\varepsilon}$$is the strain rate and parameter k is a dimensionless constant that adjusts for the inefficiency of randomly oriented faults to accommodate strain. Note that Kostrov (1974) chooses k = 1, while Anderson (1979) chooses k = 0.75, and using their relations, our calculated strain rate gives 19 and 24 nanostrain yr−1, respectively (see Table 4). These strain rate values are comparable to the Negredo et al. (2002) and Palano et al. (2013) results and stressing rate of 461 and 582 Pa yr−1, respectively. Hence, for the 1.1 bar maximum stress change obtained from the Coulomb modeling on the 2004 receiver rupture (Table 3), the clock-time advance Δt for the 2004 receiver source able to generate a significant earthquake with Mw  ≥ 6 is 239 ± 22 yr using the Kostrov (1974) relation and 189 ± 17 yr using the Anderson (1979) relation. Table 4. Clock-time advance (Δt) associated with the ΔCFF and corresponding strain rate in the vicinity of the 1994 source rupture, and for the 2004 and 2016 receiver faults in the Al Hoceima region. k = 0.75 and 1 are the values used by Kostrov (1974) and Anderson (1979), respectively, to compute the strain rate and the conditional probability Pc is calculated over 10 yr for earthquakes with Mw  > 6. Seismic  Strain rate  $$\dot{\tau }$$ stress rate  ΔCFFmax    Pc no stress  Pc with stress  ruptures  (nanostrain yr−1)  (Pa yr−1)  (MPa)  $$\Delta t = \frac{{\Delta {\rm{CFF}}}}{{\tau '}}({\rm{yr}})$$  change (pre-1994)  change (post-1994)  2004  19 (k = 0.75)  461  0.11  239 ± 22  12 per cent  55 per cent  2004  24 (k = 1)  582  0.11  189 ± 17      2016  19 (k = 0.75)  461  0.02  21 ± 02      2016  24 (k = 1)  582  0.02  17 ± 02      Seismic  Strain rate  $$\dot{\tau }$$ stress rate  ΔCFFmax    Pc no stress  Pc with stress  ruptures  (nanostrain yr−1)  (Pa yr−1)  (MPa)  $$\Delta t = \frac{{\Delta {\rm{CFF}}}}{{\tau '}}({\rm{yr}})$$  change (pre-1994)  change (post-1994)  2004  19 (k = 0.75)  461  0.11  239 ± 22  12 per cent  55 per cent  2004  24 (k = 1)  582  0.11  189 ± 17      2016  19 (k = 0.75)  461  0.02  21 ± 02      2016  24 (k = 1)  582  0.02  17 ± 02      View Large From the seismicity catalogue, the conditional probability for a specified time interval depends only on the time interval between large earthquakes ΔT, and the long-term regional seismicity rate. The conditional probability for earthquake triggering with Mw  > 6 for the Al Hoceima region is given by (Cornell 1968):   $${\rm {Pc}} = 1-{{\rm e}^{ - \lambda \Delta {\rm{T}}}}$$ (9)where λ is the seismicity rate with magnitude M > 6, and ΔT is the elapsed time since the most recent large earthquake (M > 6) obtained from the Moroccan seismicity catalogue (Jabour, personal communication, 2014). The value of λ is obtained from the Gutenberg–Richter (G-R) law constrained by the parameters a and b. Here, we also observe that λ is different before and after the 1994 stress perturbation (λ ranges between 0.012 and 0.09 following the 1994 earthquake). Based on Cornell & Winterstein (1988) hypothesis, the conditional probability related to seismic recurrence models appears to be the most appropriate (i) where the hazard is dominated by the nearest fault segment, and (ii) in the absence of slip rate and strain rate on each fault segment. We show in Table 4 that after the 1994 earthquake, the 2004 earthquake fault is under a high value of clock-time advance, for example, Δt = 239 ± 22 yr or Δt = 189 ± 17 yr according to Kostrov (1974) and Anderson (1979) formula's, respectively. However, the conditional probability for a specified time interval of 10 yr to have an earthquake with Mw  > 6 rises from 12 per cent to 55 per cent(from eq. 9, see also Table 4), in agreement with the seismic rate activity. The change in pore-fluid pressure along ruptures induced by successive earthquakes results in a cluster of large seismic events (Mw  > 6) during a short period of time (∼22 yr). The conditional probability Pc is similar in the 2004 and the 2016 seismogenic areas, since we consider the regional probability condition. In his study of the 1992 Landers California earthquake (Mw 7.3), Hardebeck (2004) shows no significant difference in the uncertainty between the conditional probability obtained from stress change and that based on the G-R distribution. In our case, Pc depends only on the G-R seismicity temporal distribution, while other studies use different approaches based mainly on stress drop. DISCUSSION AND CONCLUSIONS A sequence of three earthquakes occurred in the Rif Mountains and nearby Alboran Sea in 1994 (Mw 6.0), 2004 (Mw 6.4) and 2016 (Mw 6.5). The static stress change modeling (ΔCFF in undrained condition) suggest a fault-rupture interaction with stress loading located on the selected receiver faults. The poroelastic relaxation (drained condition) and the coupled short-term poroelastic stress transfer help us to understand the seismic migration induced by the pore-fluid diffusion. Aftershock sequences of the three earthquakes correlate well with the ΔCFF distribution which confirms the role of pore fluid in the triggering of post-1994 earthquakes (Figs S2 and S3, Supporting Information). Besides the fault-rupture complexity, the modeling parameters require two levels of friction coefficient with μ΄ < 0.4 correlating with significant pore-fluid diffusivity, and μ΄ = 0.4 in fault zones with limited diffusivity (Fig. 3b). The stress change and related pore-fluid diffusion may explain the ∼10–12 yr interval and the seismic sequence migration. The role of strain rate and its impact on stress change and pore-fluid diffusion combined with the permeability along fault-rupture zones in the Al Hoceima region is crucial in the comprehension of the time delay between the earthquake sequences. The 1994 and the 2004 earthquake ruptures illustrate the stress level change, related value of friction coefficient and role of pore-fluid diffusion in conjugate fault geometry. Although the accuracy in location of the 1994 and 2016 aftershocks is limited due to the azimuthal gap and absence of near-field seismic stations, the distribution of seismic sequences concurs with the ΔCFF and cumulative loading areas (Fig. S3, Supporting Information). However, the incomplete seismicity catalogue with precise earthquake locations in the Al Hoceima region prevents a suitable study on the role of fluid pore pressure before the 1994 earthquake sequence. The earthquake fault locations inland retrieved from the InSAR analysis of coseismic and after slip surface deformation agree with the aftershock distribution. The limited distance between earthquake ruptures (<25 km) and fault geometries with strike-slip mechanisms also promote the stress transfer and failure on fixed fault planes. Earthquake faults in the Al Hoceima region are blind with basically no geomorphological signature at the surface (Tahayt et al.2009). Therefore, fault parameters such as slip-per-event; long-term active deformation and slip rate are missing in our study. The static strain release by the 1994 event induced a high pore pressure change with fluid flow on the 2004 rupture area. This hydrological phenomenon affects the fault zone permeability and promotes the failure of the 2004 event. Wang (2000) uses the correlation between fluid migration and rock permeability to explain the link between the two phenomena; he points out that if the pore pressure becomes too large, earthquakes occur and will increase permeability with groundwater fluid flow. For a strain rate ranging between 19 and 24 nanostrain yr−1 (Table 4), we observe that the regional aftershock frequency following the 2004 earthquake (Fig. 5b) is in good agreement with the simulated aftershock frequency based on the pore-fluid diffusion hypothesis (Bosl & Nur 2002). Pore-fluid effects comparable to our case study at the intersection between the 1994 and 2004 ruptures (Fig. 4) are also observed for conjugate earthquake ruptures during the Superstition Hills earthquakes (Hudnut et al.1989; Scholz 1990). In addition, the decay in the 2004 aftershock activity includes variable seismicity rate probably due to the pore-fluid diffusion. A similar behaviour is observed during the 1966 Parkfield-Cholame earthquake where aftershock productivity and related fluid pore pressure have a direct effect on rock strength (Nur & Booker 1972). Depending on the geological structures, substratum permeability and seismicity rate, the 2004 and 2016 earthquakes could have been predictable by the Coulomb modeling taking into account the pore-fluid effect (undrained and drained conditions). In fact, the occurrence of the 2016 January 21 Mw 5.0 foreshock (4 d before the main shock) may have allowed the fluids to migrate across the epicentral area promoting the 2016 earthquake rupture. Comparable phenomenon with foreshocks and fluid migration across a fault zone is described for the L’Aquila earthquake sequence (Lucente et al.2010). The timescale of post-stress redistribution for the 1994, 2004 and 2016 Al Hoceima earthquakes is larger, for instance, than the ∼11 hr Superstitions Hills sequences (Mw 6.2 and 6.6), suggesting different diffusion processes probably controlled by the permeability along the fault zone. Following the 1994 earthquake sequence, the probability for triggering an Mw  > 6 earthquake within 10 yr interval increases to 55 per cent with respect to the 12 per cent pre-1994 period (Table 4). With the computed 239 ± 22 yr clock-time advance for large earthquakes (Mw  > 6) on the 2004 rupture, the seismic strain rate and ΔCFF explains the 10–12 yr delay and the 55 per cent probability of promoting failure in the Rif Mountains. Acknowledgements JK benefits from a scholarship from the Algerian Ministry of Higher Education and Scientific Research (MESRS). The authors wish to thank Nacer Jabour (CNRST Rabat) and Abdelilah Tahayt (University Mohamed V, Rabat) for the access to the seismicity database. We are thankful to Zoe Shipton and Mark Stillings for their comments on an early version of the manuscript, and to Ross Stein and Shinji Toda for sharing the Coulomb 3.4 software. We are grateful to two anonymous reviewers for their contributions in improving an early version of the manuscript. This research programme was funded by the Direction of Research at MESRS and the DERCI-CNRS with INSU-UMR 7516 IPG Strasbourg. Some figures were prepared using the public domain GMT software (Wessel & Smith 1998). REFERENCES Akoglu A.M., Cakir Z., Meghraoui M., Belabbes S., El Alami S.O., Ergintav S., Akyüz H.S., 2006. The 1994–2004 Al Hoceima (Morocco) earthquake sequence: conjugate fault ruptures deduced from InSAR, Earth planet. Sci. Lett. , 252, 467– 480. Google Scholar CrossRef Search ADS   Anderson J.G., 1979. Estimating the seismicity from geological structure for seismic-risk studies, Bull. seism. Soc. Am. , 69, 135– 158 Beeler N.M., Simpson R.W., Hickman S.H., Lockner D.A., 2000. Pore fluid pressure, apparent friction, and Coulomb failure, J. geophys. Res. , 105( B11), 25 533–25 542. Google Scholar CrossRef Search ADS   Bezzeghoud M., Buforn E., 1999. Source parameters of the 1992 Melilla (Spain, Mw = 4.8), 1994 Al Hoceima (Morocco Mw = 5. 8) and Mascara (Algeria, Mw 5.7) earthquakes and seismotectonic implications, Bull. seism. Soc. Am. , 89, 359– 372. Bosl, W.J., Nur A., 2002. Aftershocks and pore fluid diffusion following the 1992 Landers earthquake, J. geophys. Res. , 107( B12), 2366, doi:10.1029/2001JB000155. Google Scholar CrossRef Search ADS   Cakir Z., Meghraoui M., Akoglu A.M., Jabour N., Belabbes S., Ait-Brahim L., 2006. Surface deformation associated with the Mw 6.4, February 24, 2004 Al Hoceima (Morocco) earthquake deduced from InSAR: implications for the active tectonics along North Africa, Bull. seism. Soc. Am. , 96, 1– 10. Google Scholar CrossRef Search ADS   Calvert A., Gómez F., Seber D., Barazangi M., Jabour N., Ibenbrahim A., Demnati A., 1997. An integrated geophysical investigation of recent seismicity in the Al-Hoceima region of North Morocco, Bull. seism. Soc. Am. , 87, 637– 651. Cetin E., 2015. Analysis and modeling of crustal deformation using InSAR time series along selected active faults within the Africa-Eurasia convergence zone, PhD thesis , University of Strasbourg, EOST, 135 pp. Chalouan A., Michard A., El Kadiri Kh., Negro F., Frizon de Lamotte D., Soto J.I., Saddiqi O., 2008. The Rif Belt, in Continental Evolution: The Geology of Morocco , pp. 203– 302, eds Michard A.et al.  , Springer. Google Scholar CrossRef Search ADS   Cherkaoui T.-E., Hatzfeld D., Jebli H., Medina F., Caillot V., 1990. Etude microsismique de la région d'Al Hoceima, Bull. Inst. Sci. Rabat , 14, 25– 34. Chéry J., Merkel S., Bouissou S., 2001. A physical basis for time clustering of large earthquakes, Bull. seism. Soc. Am. , 91( 6), 1685– 1693. Google Scholar CrossRef Search ADS   Cocco M., Rice J.R., 2003. Pore pressure and poro-elasticity effects in Coulomb stress analysis of earthquake interactions, J. geophys. Res. , 107, doi:10.1029/2000JB000138. Console R., Murru M., Falcone G., 2010. Probability gains of an epidemic-type aftershock sequence model in retrospective forecasting of M ≥ 5 earth-quakes in Italy, J. Seismol. , 14, 9– 26. Google Scholar CrossRef Search ADS   Cornell C.A., 1968. Engineering seismic risk analysis, Bull. seism. Soc. Am. , 58, 1583– 1606. Cornell C.A., Winterstein S.R., 1988. Temporal and magnitude dependence in earthquake recurrence models, Bull. seism. Soc. Am. , 78( 4), 1522– 1573. Das S., Scholz C.H., 1981. Theory of time-dependent rupture in the Earth, J. geophys. Res. , 86( B7), 6039–6051. Google Scholar CrossRef Search ADS   Demets C., Gordon R.G., Argus D.F., 2010. Geologically current plate motions, Geophys. J. Int. , 181, 1– 80. Google Scholar CrossRef Search ADS   El Alami S.O., Tadili B.A., Cherkaoui T.E., Medina F., Ramdani M., Ait Brahim L., Harnaff M., 1998. The Al Hoceima earthquake of May 26, 1994 and its aftershocks: a seismotectonic study, Ann. Geofis. , 41, 519– 537. Fernandez-Ibanez F., Soto J.I., Zoback M.D., Morales J., 2007. Present-day stress field in the Gibraltar Arc (western Mediterranean), J. geophys. Res. , 112, B08404, doi:10.1029/2006JB004683. Google Scholar CrossRef Search ADS   Grevemeyer I., Gràcia E., Villaseñor A., Leuchters W., Watts A.B., 2015. Seismicity and active tectonics in the Alboran Sea, Western Mediterranean: constraints from an offshore-onshore seismological network and swath bathymetry data, J. geophys. Res. , 120, doi:10.1002/2015JB012073. Guo Z., Ogata Y., 1997. Statistical relations between the parameters of aftershocks in time, space, and magnitude, J. geophys. Res. , 102( B2), 2857–2873. Google Scholar CrossRef Search ADS   Hardebeck J.L., 2004. Stress triggering and earthquake probability estimates, J. geophys. Res. , 109, B04310, doi:10.1029/2003JB002437. Hatzfeld D., Frogneux M., Girardin D., 1977. Etude de la sismicite dans la region de I'arc de Gibraltar, Bull. Soc. geol. Fr. , S7-XIX( 4), 741– 747. Google Scholar CrossRef Search ADS   Hudnut W., Seeber L., Pacheco J., 1989. Cross-fault triggering in the November 1987 Superstition Hills earthquake sequence, Southern California, Geophys. Res. Lett. , 16( 2), 199– 202. Google Scholar CrossRef Search ADS   Jonsson S., Segall P., Pedersen R., Bjornsson G., 2003. Post-earthquake ground movements correlated to pore-pressure transients, Nature , 424, 179– 183. Google Scholar CrossRef Search ADS PubMed  King, G.C.P., Stein, R.S., Lin, J., 1994. Static stress changes and the triggering of earthquakes, Bull. seism. Soc. Am. , 84( 3), 935– 953. Kirkpatrick J.D., Shipton Z.K., Evans J.P., Micklethwaite S., Lim S.J., McKillop P., 2008. Strike-slip fault terminations at seismogenic depths: the structure and kinematics of the Glacier Lakes fault, Sierra Nevada United States, J. geophys. Res. , 113, B04304, doi:10.1029/2007JB005311. Google Scholar CrossRef Search ADS   Kostrov V.V., 1974. Seismic moment and energy of earthquakes, and seismic flow of rocks, Izv. Acad. Sci. USSR Phys. Solid Earth , 1, 23–44 (Eng. Transl.). Koulali A.et al.  , 2011. New GPS constraints on active deformation along the Africa-Iberia plate boundary, Earth planet. Sci. Lett. , 308, 211– 217. Google Scholar CrossRef Search ADS   Lucente F.P., Gori P.D., Margheriti L., Piccinini D., Bona M.D., Chiarabba C., Agostinetti N.P., 2010. Temporal variation of seismic velocity and anisotropy before the 2009 Mw 6.3 L’Aquila earthquake, Italy, Geology , 38, 1015– 1018. Google Scholar CrossRef Search ADS   Medina F., 1995. Present-day state of stress in northern Morocco from focal mechanism analysis, J. Struct. Geol. , 17( 7), 1035– 1046. Google Scholar CrossRef Search ADS   Meghraoui M., Pondrelli S., 2012. Active faulting and transpression tectonics along the plate boundary in North Africa, Ann. Geophys. , 55( 5), 955– 967. Morel J.-L., Meghraoui M., 1996. The Goringe–Alboran–Tell tectonic zone, a transpression system along the Africa–Eurasia plate boundary, Geology , 24, 755– 758. Google Scholar CrossRef Search ADS   Narteau C., Byrdina S., Shebalin P., Schorlemmer D., 2009. Common dependence on stress for the two fundamental laws of statistical seismology, Nature , 462( 7273), 642– 645. Google Scholar CrossRef Search ADS PubMed  Negredo A.M., Bird P., de Galdeano C.-S., Buforn E., 2002. Noetectonic modeling of the Ibero-Maghrebian region, J. geophys. Res. , 107( B11), ETG10-1– ETG10-15. Google Scholar CrossRef Search ADS   Nur A., Booker J.R., 1972. Aftershocks caused by pore fluid flow?, Science , 175( 4024), 885– 887. Google Scholar CrossRef Search ADS PubMed  Okada Y., 1992. Internal deformation due to shear and tensile faults in a half-space, Bull. seism. Soc. Am. , 82( 2), 1018– 1040. Palano M., González P.J., Fernández J.A., 2013. Strain and stress fields along the Gibraltar orogenic arc: constraints on active geodynamics, Gondwana Res. , 23, 1071– 1088. Google Scholar CrossRef Search ADS   Piombo A., Martinelli G., Dragoni M. 2005. Post-seismic fluid flow and Coulomb stress changes in a poro-elastic medium, Geophys. J. Int. , 162, 507– 515. Google Scholar CrossRef Search ADS   Pytharouli S.I., Lunn R.J., Shipton Z.K., Kirkpatrick J.D., do Nascimento A.F., 2011. Microseismicity illuminates open fractures in the shallow crust, Geophys. Res. Lett. , 38, L02402, doi:10.1029/2010GL045875. Google Scholar CrossRef Search ADS   Reasenberg P.A., Simpson R.W., 1992. Response of regional seismicity to the static stress change produced by the Loma Prieta earthquake, Science , 255, 1687– 1690. Google Scholar CrossRef Search ADS PubMed  Rice J.R., Cleary M.P., 1976. Some basic stress diffusion solutions for fluid-saturated elastic porous media with compressible constituents, Rev. Geophys. Space Phys. , 14( 2), 227– 241. Google Scholar CrossRef Search ADS   Shapiro S.A., Huenges E., Borm G., 1997. Estimating the crust permeability from fluid-injection-induced seismic emission at the KTB site, Geophys. J. Int. , 131, F15– F18. Google Scholar CrossRef Search ADS   Scholz C.H., 1990. The Mechanics of Earthquakes and Faulting , Cambridge Univ. Press, p. 439. Stein R.S., Barka A., Dieterich J.H., 1997. Progressive failure on the North Anatolian fault since 1939 by earthquake stress triggering, Geophys. J. Int. , 128, 594– 604. Google Scholar CrossRef Search ADS   Stich D., Serpelloni E., Mancilla F., Morales J., 2006. Kinematics of the Iberia-Maghreb plate contact from seismic moment tensors and GPS observations, Tectonophysics , 426, 295– 317. Google Scholar CrossRef Search ADS   Stich D., Martin R., Morales J., 2010. Moment tensor inversion for Iberia-Maghreb earthquakes 2005–2008, Tectonophysics , 483, 390– 398. Google Scholar CrossRef Search ADS   Tahayt A. et al., 2009. The Al Hoceima (Morocco) earthquake of 24 February (2004) analysis and interpretation of data from ENVISAT ASAR and SPOT5 validated by ground-based observations, Remote Sens. Environ. , 113, 306– 316. Google Scholar CrossRef Search ADS   Terzaghi K., 1925. Erdbaumechanik auf Bodenphysikalischer Grundlage. Franz Deuticke , Liepzig-Vienna. Timoulali Y., Hahou Y., Jabour N., Merrouch R., El Kharrim A., 2014. Main features of the deep structure by local earthquake tomography and active tectonics: case of Rif Mountain (Morocco) and Betic Cordillera (Spain), J. Seismol. , 18( 2), 221– 234. Google Scholar CrossRef Search ADS   Toda S., Stein R.S., Sevilgen V., Lin J., 2011. Coulomb 3.3 graphic-rich deformation and stress-change software for earthquake, tectonic, and volcano research and teaching, User guide, U.S. Geol. Surv . Open-File Rept. 2011-1060, 63 pp. Available at: http://pubs.usgs.gov/of/2011/1060/, last accessed October 2016. Turkaya S., Toussaint R., Ericksen F.K., Zecevic M., Daniel G., Flekkøy E.G., Måløy M., 2015. Bridging aero-fracture evolution with the characteristics of the acoustic emissions in a porous medium, Front. Phys. , 3, doi:10.3389/fphy.2015.00070. Utsu T., 1969. Aftershocks and earthquake statistics (I) - Some parameters which characterize an aftershock sequence and their interrelations, J. Fac. Sci. Hokkaido Univ., Ser. , VII( 3), 121– 195. Vallée M., 2016. Automatic determination of source parameters using the SCARDEC method. Available at: http://geoscope.ipgp.fr/index.php/en/catalog/earthquake-description?seis=us10004gy9, last accessed 22 September 2016. Wang H.F., 2000. Theory of Linear Poroelasticity with Applications to Geomechanics and Hydrogeology , pp. 304, Princeton Univ. Press. Wessel P., Smith W.H.F., 1998. New, improved version of generic mapping tools released, EOS, Trans. Am. geophys. Un. , 79( 47), 579, doi:10.1029/98EO00426. Google Scholar CrossRef Search ADS   Wyss M., Habermann R.E., 1988. Precursory seismic quiescence, Pure appl. Geophys. , 126, 319– 332. Google Scholar CrossRef Search ADS   Wyss M., Wiemer S., 2000. Change in the probabilities for earthquakes in Southern California due to the Landers M 7.3 earthquake, Science , 290, 1334– 1338. Google Scholar CrossRef Search ADS PubMed  SUPPORTING INFORMATION Supplementary data are available at GJI online. Figure S1. (a) ΔCFF modeling on the 2004 fault plane due to the 1994 fault rupture, the computation are made for μ΄ = 0.2, the black stars indicate the 2004 hypocentre, (b) ΔCFF modeling on the 2016 fault plane due to the cumulative 1994 and 2004 earthquake, the grey stars indicate the 2016 hypocentre, note that the receiver planes are tapered to 40 fault patches system. Figure S2. ΔCFF modeling on optimally oriented fault planes caused by the 2004 earthquake at depth range [5–15 km] for: (a) μ΄ = 0.2, (b) μ΄ = 0.4, (c) μ΄ = 0.6 and (d) μ΄ = 0.8, the aftershocks distributions are represented by black circles, the white star represents the 2004 epicentre and the grey star is the 2016 epicentre. The increase of pressure could trigger events in regions where the Coulomb stress predict and absence of activity. Figure S3. Full poroelastic relaxation of the 1994 and the 2016 earthquakes from the undrained state to the fully drained state, the computation is obtained by linear elasticity theory with appropriate values of undrained and drained Poisson ratio: (a) full poroelastic relaxation due to the 1994 earthquake on the 2004 fault rupture with a typical sedimentary undrained and drained Poisson ratios, respectively, (b) full poroelastic relaxation due to the 1994 earthquake on the 2004 fault rupture with an extreme undrained and drained Poisson ratios, respectively, according the Bosl and Nur hypothesis, this value might be typical for upper crust when fractures exists, (c) ΔCFF ( red + green) and ΔP/B profiles along the 2004 rupture related to the full poroelastic relaxation of the 1994 earthquake, a significant change of ΔP/B is obtained for (vu, v) = (0.31,0.15), the yellow stars indicate the 2004 epicentre location on the axe, (d) full poroelastic relaxation due to the 2016 earthquake on fixed planes strike/dip/rake = 195°/78°/19° for (vu, v) = (0.31,0.15) and (e)full poroelastic relaxation due to the 2016 earthquake on fixed planes strike/dip/rake = 195°/78°/19° for (vu, v) = (0.31,0.25). Note that the white stars indicate the 2016 epicentre location. Figure S4. Angular relationship between pre-seismic local stress field and the 2004 fault plane acting on the central Rif block, the maximum horizontal stress (green) and the minimum horizontal (yellow) stress directions acting of the 2004 rupture at depth are obtained from stress estimation based on inversion of focal mechanism and GPS data obtained by several authors (Akoglu et al.2006; Fernandez-Ibanez et al.2007; Tahayt et al.2009), σ and τ denote the normal and shear stress according to the state of stress, the red black and the red line indicates the SS surface fault of the 1994 and 2004 earthquakes, respectively, the tectonic models used for this study are from Tahayt et al. (2009), where arrows indicate the relative movements of these blocks with respect to Africa plate, the 1994 and the 2004 rupture are from Akoglu et al. (2006). The majors reverse faults trace are represented by an appropriate symbol and the major strike-slip fault traces represented by lines. Table S1. Coulomb stress change caused by the 2004 and the 2016 earthquakes on aftershocks planes for μ΄ = 0.4, the strike/dip/rake represent the best geometry related to the optimally oriented ΔCFF loading. SF represents the source fault. Appendix A: Constitutive equations of poro-elasticity. Appendix B: Relationship between the poro-elastic model and the 2004 earthquake triggering. Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper. © The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# The Al Hoceima earthquake sequence of 1994, 2004 and 2016: Stress transfer and poroelasticity in the Rif and Alboran Sea region

, Volume 212 (1) – Jan 1, 2018
12 pages

/lp/ou_press/the-al-hoceima-earthquake-sequence-of-1994-2004-and-2016-stress-20SnwUr3PX
Publisher
The Royal Astronomical Society
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx385
Publisher site
See Article on Publisher Site

### Abstract

Abstract The 2016 January 25 earthquake (Mw 6.3) follows in sequence from the1994 May 26 earthquake (Mw 6.0) and the 2004 February 24 earthquake (Mw 6.4) in the Rif Mountains and Alboran Sea. The earlier two seismic events which were destructive took place on inland conjugate faults, and the third event occurred on an offshore fault. These earthquake sequences occurred within a period of 22 yr at ∼25 km distance and 11–16-km depth. The three events have similar strike-slip focal mechanism solutions with NNE-SSW trending left-lateral faulting for the 1994 and 2016 events and NW-SE trending right-lateral faulting for the 2004 event. This shallow seismic sequence offers the possibility (i) to model the change in Coulomb Failure Function (ΔCFF with low μ΄ including the pore pressure change) and understand fault-rupture interaction, and (ii) to analyse the effect of pore fluid on the rupture mechanism, and infer the clock-time advance. The variation of static stress change has a direct impact on the main shock, aftershocks and related positive lobes of the 2004 earthquake rupture with a stress change increase of 0.7–1.1 bar. Similarly, the 2004 main shock and aftershocks indicate loading zones with a stress change (>0.25 bar) that includes the 2016 earthquake rupture. The tectonic loading of 19–24 nanostrain yr−1 obtained from the seismicity catalogue of Morocco is comparable to the 5.0 × 1017 N·m yr−1 seismic strain release in the Rif Mountains. The seismic sequence is apparently controlled by the poroelastic properties of the seismogenic layer that depend on the undrained and drained fluid conditions. The short interseismic period between main shocks and higher rate of aftershocks with relatively large magnitudes (4 <  Mw < 5.5) implies the pore-fluid physical effect in undrained and drained conditions. The stress-rate ranges between 461 and 582 Pa yr−1 with a ΔCFF of 0.2–1.1 bar. The computed clock-time advance reaches 239 ± 22 yr in agreement with the ∼10 yr delay between main shocks. The calculated static stress change of 0.9–1.3 bar, under pore-fluid stimulus added with well-constrained geodetic and seismic strain rates are critical for any seismic hazard assessment. Seismicity and tectonics, Dynamics: seismotectonics, Fractures, faults, and high strain deformation zones, Neotectonics INTRODUCTION Three significant earthquakes have occurred on 1994 May 26 (Mw 6.0), 2004 February 24 (Mw 6.4) and 2016 January 25 (Mw 6.5) in the Rif Mountains of Morocco and southern Alboran Sea area within a period of 22 yr (Fig. 1). The earlier two earthquakes caused severe damage due to their location inland, but the third offshore event was only felt on the nearby Moroccan coastline. Aftershock distribution (El Alami et al.1998; Bezzeghoud & Buforn 1999) and surface deformation as deduced from Synthetic Aperture Radar Interferogram (InSAR) (Akoglu et al.2006; Cakir et al.2006; Tahayt et al.2009), indicate that the 1994 and 2004 events occurred on NNE–SSW and NW–SE trending conjugate strike-slip faults, respectively. The 2016 event located about 20 km offshore is associated with a NNE–SSW trending rupture with a similar mechanism to the 1994 event. This seismic sequence is unusual in the North Africa active zones because earthquake ruptures are within ∼25 km area and the time interval between main shocks is about 10–12 yr. Figure 1. View largeDownload slide Seismicity of the Al Hoceima region showing the seismic sequence of 1994 May 26 (El Alami et al.1998), 2004 February 24 (Tahayt et al.2009) and 2016 January 25 (CSEM, http://www.emsc-csem.org/#2). The 1994 aftershocks are in green, the 2004 aftershocks are in blue and 2016 aftershocks are in grey. Focal mechanisms are Harvard–CMT (Table 1). Inset represents the plate boundary in the Alboran Sea with convergence rate in mm yr−1 (Koulali et al.2011). Figure 1. View largeDownload slide Seismicity of the Al Hoceima region showing the seismic sequence of 1994 May 26 (El Alami et al.1998), 2004 February 24 (Tahayt et al.2009) and 2016 January 25 (CSEM, http://www.emsc-csem.org/#2). The 1994 aftershocks are in green, the 2004 aftershocks are in blue and 2016 aftershocks are in grey. Focal mechanisms are Harvard–CMT (Table 1). Inset represents the plate boundary in the Alboran Sea with convergence rate in mm yr−1 (Koulali et al.2011). In this paper, we first present the seismic sequence and suggest a fault-rupture interaction using Coulomb modeling on fixed planes. Secondly, the computation using optimally oriented planes is added to constrain aftershocks distribution. The three main shocks and related aftershocks appear to be closely related and their location implies a stress transfer with triggering. We show that the modeled stress distribution and seismicity rate change suggest a pore-fluid effect correlated with elastic dislocation in undrained and drained conditions. The Coulomb Failure Function change (ΔCFF) and pore-fluid flow seem to control the 10–12 yr recurrence of main seismic events with a clock advance. The plate boundary tectonic condition, the seismicity rate change and poroelastic properties of the seismogenic crust seem to play a significant role in the triggering of earthquakes in the Rif Mountains and Alboran Sea. SEISMOTECTONIC SETTING The seismicity of the Rif Mountains and Alboran Sea is due to the convergence between Africa and Eurasia (Iberia) in the western Mediterranean. The E–W trending Rif Mountains run along the northern coast of Morocco forming the southern branch of the Betic-Rif arc that includes the Alboran Sea and belongs to the transpression plate boundary system in the Western Mediterranean region (Morel & Meghraoui 1996; Meghraoui & Pondrelli 2012). Tahayt et al. (2009) interpret the region as a transrotational regime applied to the Oriental Rif block with a clockwise rotation. This complex tectonic domain also results from a Neogene subdued subduction zone with lithospheric delamination where the Alboran Sea appears as an oceanic microplate (Calvert et al.1997). As indicated by the three main shocks (Fig. 1 and Table 1), the present-day tectonic framework of the Al Hoceima region is dominated by a strike-slip fault regime where moment magnitudes do not exceed 6.5. Limited fault-rupture dimensions are likely due to the local structural geology made of overthrusting nappes on highly deformed continental crustal rocks (Chalouan et al.2008; Timoulali et al. 2014). Because of the limited number of local seismic stations in the Rif, the location of aftershocks of the 1994 and offshore 2016 earthquakes are poorly resolved, which is not the case for the 2004 earthquake (Tahayt et al.2009). Table 1. Physical characteristics of the earthquakes used in this study, HRV designed Harvard solution and CMT designed Centroid Moment Tensor solution (see also Fig. 1). Earthquake  Long.  Lat.  Mo (1018 nm)  Mw  U (m)  L (km)  W (km)  Strike  Dip  Rake  1994 May 26  −3.99  35.28  1.01  6.1  0.8  16  10  17  85  −7  2004 February 24  −3.99  35.142  3.0  6.4  1.0  19  16.5  340  87  −161  2016 January 25  −3.70  35.67  4.69  6.5  0.8  25  13.5  214  78  19  Earthquake  Long.  Lat.  Mo (1018 nm)  Mw  U (m)  L (km)  W (km)  Strike  Dip  Rake  1994 May 26  −3.99  35.28  1.01  6.1  0.8  16  10  17  85  −7  2004 February 24  −3.99  35.142  3.0  6.4  1.0  19  16.5  340  87  −161  2016 January 25  −3.70  35.67  4.69  6.5  0.8  25  13.5  214  78  19  View Large Several authors have studied the seismicity and suggested moment tensor solutions of major earthquakes showing a NNW–SSE contraction stress regime in a predominantly strike-slip faulting domain associated with normal and thrust mechanisms (Hatzfeld et al.1977; Cherkaoui et al. 1990; Medina 1995; Stich et al.2006, 2010; Palano et al.2013). Although the seismicity may appear diffuse in the Rif-Alboran Sea, the 1994, 2004 and 2016 seismic events reveal a clear migration of earthquake ruptures (Fig. 1). The seismicity and tectonics of the region indicate a clear correlation between continental and offshore faults from both the Betics and Rif Mountains towards the Alboran Sea (Grevemeyer et al.2015). From detailed bathymetry and seismicity distribution, Gevemeyer et al. (2015) identify a fault zone crossing the Alboran ridge at the location of the 2016 January 25 earthquake (Mw 6.5). Using InSAR, Cakir et al. (2006) and Akoglu et al. (2006) constrain the coseismic earthquake surface deformation and provide the rupture parameters of the N23° E trending 1994 and N45° W trending 2004 earthquakes using elastic dislocations. They suggest that the two main shocks occurred on blind and conjugate strike-slip faults, with left- and right-lateral slip, respectively. The time-series analysis of SAR data (Envisat) for 7 yr following the 2004 earthquake shows that post-seismic deformation reaches up to 4 cm at the surface and infers 0.3 m displacement at shallow depth (<7 km), mainly above the high coseismic slip patches which can be explained by 5.0 × 1017 N·m (Mw 5.73) cumulative moment release (Cetin 2015). The 2016 earthquake rupture, located about 20 km further north, is modeled for the source time function to obtain a strike-slip faulting mechanism with a NNE–SSW main rupture in agreement with aftershock distribution (Vallée 2016). The three focal mechanism solutions and the 1994 and 2004 rupture geometries inferred from surface deformation are used as an input for the ΔCFF modeling (Table 1). MODELING AL HOCEIMA SEQUENCE BY COULOMB FAILURE FUNCTION The Rif region of northern Morocco experienced a seismic sequence with three moderate to large earthquakes within 22 yr. The sequence suggests earthquake triggering, fault interaction and stress transfer as observed in other earthquake areas (Hudnut et al.1989; Stein et al.1997). The case study uses the applied stress change calculated as ΔCFF (Reasenberg & Simpson 1992; King et al.1994) expressed by:   $$\Delta {\rm{CFF}}\ = \ \Delta \tau \ -\ \mu \ \left( {\ \Delta {\sigma _n} - \Delta P} \right)$$ (1)  $$\Delta {\rm{CFF}}\ = \Delta \tau \ - \ \mu '\Delta {\sigma _n}$$ (2)where τ is the shear stress, σn is the normal stress (compression positive), P is the pore-fluid pressure, μ and μ΄ are the coefficient of friction and effective coefficient of friction, respectively and Δ refers to changes during the earthquake. The apparent friction is given by (Reasenberg & Simpson 1992)   $$\mu ' = \ \mu \ \left( {1\ - B} \right)$$ (3)where B is the Skempton coefficient which defines the relation between the stress change and pore pressure change (Beeler et al.2000)   $$B\ = \ \ \frac{{\Delta P}}{{\Delta {\sigma _m}}} = \ 3\ \frac{{\Delta P}}{{\Delta {\sigma _{kk}}}}$$ (4)where Δσm is the mean stress change and Δσkk is the sum over the diagonal elements of the stress tensor. It is important to note that for an isotropic model, the apparent friction coefficient used in triggered seismicity is defined by the combination of pore pressure and friction coefficient:   $$\mu ' = \ \mu \ \left( {1 - \frac{{{\rm{\Delta }}P}}{{{\rm{\Delta }}{\sigma _n}}}} \right)$$ (5)Substituting eqs (5) and (4) in eq. (2), we obtain   $${\rm{\Delta CFF}}\ = \ {\rm{\Delta }}\tau + \ \mu \ \left( {{\rm{\Delta }}{\sigma _n} - B{\rm{\Delta }}{\sigma _m}} \right)$$ (6a) Beeler et al. (2000) suggest that due to the pore-fluid effect and for an isotropic poroelastic model, the two expressions defined in eqs (2) and (6) yield different results in some modeling configurations. The variable effective friction coefficient related to the variation of the Skempton coefficient B and the pore pressure change along the fault zone gives more realistic solutions especially at high pore pressure change than the imposed constant effective friction commonly used in the Coulomb stress modeling (see Appendices A and B for more details). The variation of B must be considered when different porosity and different diffusive processes are present in the fault zone (Scholz 1990). From the structural point of view, fluid migration in the host rock may occur on nearby subsidiary fractures linked to bounding faults directly related to earthquake and aftershock behaviour (Kirkpatrick et al.2008). The cumulative moment for Al Hoceima region (inland and offshore) reaches 1.1 × 1019 N·m, since 1994, and the seismic strain release tested on a fault network reveals a low effective coefficient of friction on the fault (≤0.2) to model the slip rates (Negredo et al.2002). In order to model the active tectonics of the Ibero-Maghrebian region, a strain rate of 15–40 nanostrain yr−1 is obtained using slip rates on faults and Global Positioning System (GPS) data (Negredo et al.2002; Koulali et al.2011; Palano et al.2013). In their calculation of frictional strength, Negredo et al. (2002) assume no cohesion in the media. In our work, the modeling is performed using Coulomb 3.4 software (Toda et al.2011) based on the conversion of DC3D subroutines (Okada 1992) to calculate the ΔCFF. The static stress change is computed on both fixed receiver earthquake ruptures and optimally oriented faults using an effective coefficient of friction μ΄ = 0.4, and fault parameters summarized in Table 2. The computed stress changes are presented in Table 3, where the value of static stress change according to the receiver fault represents the ΔCFFmax and also expresses the Coulomb stress drop on the source faults. Table 2. Rupture parameters of significant earthquakes from InSAR results (Akoglu et al.2006) and source time function (http://geoscope.ipgp.fr/index.php/en/catalog/earthquake-description?seis = us10004gy9) used for the Coulomb stress transfer modeling. Long_c and Lat_c indicate the centre of each dislocation. EQ  Long_c (°)  Lat_c (°)  Top-depth (km)  Bot-depth (km)  L (km)  W (km)  Strike (°)  Dip (°)  Rev.slip (m)  Right_lat. slip (m)  M0 (dyne·cm) ×1025  1994  −4.01  35.17  1.00  11.00  16.00  10.04  17  85  −0.04  −0.32  1.7  2004  −3.97  35.13  1.00  16.00  20.63  15.02  340  87  −0.22  0.63  6.6  2016  −3.84  35.5  1.00  14.47  28.70  13.50  205  86  0.2  −0.57  7.5  EQ  Long_c (°)  Lat_c (°)  Top-depth (km)  Bot-depth (km)  L (km)  W (km)  Strike (°)  Dip (°)  Rev.slip (m)  Right_lat. slip (m)  M0 (dyne·cm) ×1025  1994  −4.01  35.17  1.00  11.00  16.00  10.04  17  85  −0.04  −0.32  1.7  2004  −3.97  35.13  1.00  16.00  20.63  15.02  340  87  −0.22  0.63  6.6  2016  −3.84  35.5  1.00  14.47  28.70  13.50  205  86  0.2  −0.57  7.5  View Large Table 3. Shear, normal and Coulomb stress change for the major earthquakes. Fault geometries used in the Coulomb stress modeling are defined in Table 2. SF and RF are the source and receiver faults, respectively. EQ (SF-RF)  Calc_Location  Receiver faults (°)  Stress computation (bar)    Long (°)  Lat (°)  Z (km)  Strike  Dip  Rake  Shear  Normal  Coulomb  1994(SF)  −4.010  35.166  7.00  340  87  −161  1.997  −8.547  −1.422  2004(RF)  −3.931  35.088  7.00        0.656  1.104  1.062  1994(SF)  −3.987  35.131  7.00  195  78  19  −8.146  −1.227  −8.637  2016(RF)  −3.870  35.429  7.00        0.093  0.048  0.112  2004(SF)  −4.008  35.166  7.00  195  78  19  −4.594  −11.509  −9.198  2016(RF)  −3.870  35.429  7.00        0.228  0.033  0.241  EQ (SF-RF)  Calc_Location  Receiver faults (°)  Stress computation (bar)    Long (°)  Lat (°)  Z (km)  Strike  Dip  Rake  Shear  Normal  Coulomb  1994(SF)  −4.010  35.166  7.00  340  87  −161  1.997  −8.547  −1.422  2004(RF)  −3.931  35.088  7.00        0.656  1.104  1.062  1994(SF)  −3.987  35.131  7.00  195  78  19  −8.146  −1.227  −8.637  2016(RF)  −3.870  35.429  7.00        0.093  0.048  0.112  2004(SF)  −4.008  35.166  7.00  195  78  19  −4.594  −11.509  −9.198  2016(RF)  −3.870  35.429  7.00        0.228  0.033  0.241  View Large In the ΔCFF modeling computed with fixed strike-slip receiver fault plane, the correlation between loading lobes and aftershock distribution suggests a close interaction between the 1994, 2004 and 2016 earthquake ruptures (Figs 2a–c, 3a and Table 3, and Figs S1a and b, S2 and Table S1, Supporting Information). Unlike the 2016 earthquake epicentre which is clearly within positive CFF lobes, the 2004 earthquake location given by the InSAR analysis (Akoglu et al.2006) is at the transition from negative to positive ΔCFF lobes. Taking the 1994 as a source fault and the 2004 rupture as a receiver fault, and using a low friction coefficient (μ΄ < 0.4 and for isotropic models, Figs 3a and b), our CFF modeling show the 2004 fault rupture clearly located in a loading zone with positive stress change. Figure 2. View largeDownload slide (a) Calculated ΔCFF with the 1994 source fault, and related aftershock distribution (see the text for explanation) and the 2004 as a receiver fault (strike/dip/rake = 340°/87°/−161°), the blue, white and black stars are epicentres of the 1999, 2004 and 2016, respectively (same symbols in figures b–d). (b) Computed ΔCFF with the 2004 as a source fault and the 2016 as a receiver fault (strike/dip/rake = 195°/78°/19°), and related aftershock distribution (see the text for explanation). (c) Computed ΔCFF with the 2016 source fault for fixed planes (strike/dip/rake = 195°/78°/19°), and related aftershock distribution (see the text for explanation). (d) Computed cumulative ΔCFF with the three source faults on a fixed planes (strike/dip/rake = 195°/78°/19°), and related aftershock distribution (see the text for explanation). The 1994 and 2004 mapped fault ruptures are from Akoglu et al. (2006), the 2016 fault-rupture model is from M. Vallée (http://geoscope.ipgp.fr/index.php/en/catalog/earthquake-description?seis=us10004gy9). Figure 2. View largeDownload slide (a) Calculated ΔCFF with the 1994 source fault, and related aftershock distribution (see the text for explanation) and the 2004 as a receiver fault (strike/dip/rake = 340°/87°/−161°), the blue, white and black stars are epicentres of the 1999, 2004 and 2016, respectively (same symbols in figures b–d). (b) Computed ΔCFF with the 2004 as a source fault and the 2016 as a receiver fault (strike/dip/rake = 195°/78°/19°), and related aftershock distribution (see the text for explanation). (c) Computed ΔCFF with the 2016 source fault for fixed planes (strike/dip/rake = 195°/78°/19°), and related aftershock distribution (see the text for explanation). (d) Computed cumulative ΔCFF with the three source faults on a fixed planes (strike/dip/rake = 195°/78°/19°), and related aftershock distribution (see the text for explanation). The 1994 and 2004 mapped fault ruptures are from Akoglu et al. (2006), the 2016 fault-rupture model is from M. Vallée (http://geoscope.ipgp.fr/index.php/en/catalog/earthquake-description?seis=us10004gy9). Figure 3. View largeDownload slide (a) ΔCFF for various effective friction coefficients (μ΄) along the 2004 rupture (strike direction) as receiver fault with strike/dip/rake = 340°/87°/−161° at 7 km depth. The profiles start at the intersection of the two cross faults and terminate at the end of the aftershock sequence (SE of epicentre area). The increases of pore pressure could trigger events in regions where reduced Coulomb stress predict and absence of activity. We show that at the epicentre area, the rupture nucleation occurs when pore fluids are redistributed, the value of μ΄ ≤ 0.4 seems to be more adaptable for a chosen receiver fault geometry. The computation is based on the effective constant friction model. (b) Stress change caused by the 1994 earthquake (source fault) along the 2004 fault zone due to (1) coseismic stress change due to elastic dislocation (red line), and (2) stress change due to the coupled poroelastic effect (green line). The maximum stress load on the 2004 fault zone are given for an isotropic model when μ΄ < 0.4 that implies Skempton coefficient B = 0.9 near the 1994 rupture and for μ΄ = 0.4 with B = 0.47 far from the 1994 rupture. Figure 3. View largeDownload slide (a) ΔCFF for various effective friction coefficients (μ΄) along the 2004 rupture (strike direction) as receiver fault with strike/dip/rake = 340°/87°/−161° at 7 km depth. The profiles start at the intersection of the two cross faults and terminate at the end of the aftershock sequence (SE of epicentre area). The increases of pore pressure could trigger events in regions where reduced Coulomb stress predict and absence of activity. We show that at the epicentre area, the rupture nucleation occurs when pore fluids are redistributed, the value of μ΄ ≤ 0.4 seems to be more adaptable for a chosen receiver fault geometry. The computation is based on the effective constant friction model. (b) Stress change caused by the 1994 earthquake (source fault) along the 2004 fault zone due to (1) coseismic stress change due to elastic dislocation (red line), and (2) stress change due to the coupled poroelastic effect (green line). The maximum stress load on the 2004 fault zone are given for an isotropic model when μ΄ < 0.4 that implies Skempton coefficient B = 0.9 near the 1994 rupture and for μ΄ = 0.4 with B = 0.47 far from the 1994 rupture. The cumulative ΔCFF from 1994 to 2016 (Fig. 2d) explains the present day and most important part of the seismicity and aftershock distribution and Coulomb stress drop in the Al Hoceima region. In order to analyse, the relationship between the main shocks and aftershocks distribution in detail, we model a static stress change caused by the 2004 and the 2016 earthquakes on some major aftershocks (see Table S1, Supporting Information). The modeling on fixed receiver planes suggest a stress load range between 0.5 and 0.8 bar for 2004 aftershock sequences with an optimal strike and dip value ranging between 207° and 298° and 66° and 84°, respectively. For the 2016 earthquake, major aftershocks are reached by a positive ΔCFF ranging between 0.1 and 0.4 bar, with an optimal value of [250°–260°] for strike and [40°–45°] for dip. In order to take into account all aftershock sequences, we compute a static stress change on optimally oriented fault planes as this approach does not need to include the focal mechanism of each rupture. As obtained by previous works on the northern Morocco tectonics (Medina 1995; Akoglu et al.2006; Fernandez-Ibanez et al.2007), a regional stress field (with a NW–SE principal stress direction as σ1) is added as a pre-existing stress field on the stress modeling. Based on the seismic tensor inversion and GPS data, the inferred stress field and maximum horizontal stress for Al-Hoceima–Alboran region is in good agreement with convergence models along the plate boundary (Demets et al.2010; Meghraoui & Pondrelli 2012). Note that the magnitude of the principal pre-existing stresses does not change the static coulomb stress modeling, because the stress levels are largely dominated by the coseismic rupture process in the near field. At seismogenic depth and for optimal failure planes, the static stress change modeling due to the 2004 earthquake on optimally oriented fault planes suggest that ∼30 per cent of aftershocks hypocentres were pushed closer to failure for high effective friction coefficients (Figs S2c and d, Supporting Information), while this percentage rises to ∼90 per cent when pore fluid are redistributed (Figs S2a and b, Supporting Information). The elastic modeling in Fig. 3(a) represents the stress change caused by the 1994 earthquake along the 2004 fault zone (with strike-slip mechanism, see also Table 3). The ΔCFF profiles are constrained by the aftershock distributions and the computation is performed for receiver faults with strike/dip/rake = 340°/87°/−161°. For the modeling procedure, we assume a 0.25 typical value of Poisson ratio with 8 × 105 bar for the Young modulus and 3.3 × 105 bar for the shear modulus in the seismogenic layer (5–15 km thickness). Here, the ΔCFF modeling require a low effective friction coefficient (μ΄ ≤ 0.4) denoting a pore-fluid effect. At the 2004 epicentre area, the modeling suggest that the rupture nucleation occurs when pore fluids are redistributed, with the μ΄ ≤ 0.4 considered as an optimal value also explains the nucleation process for a chosen receiver fault geometry. The increase of pore pressure may trigger seismic events in regions where reduced Coulomb stress due to the high effective friction coefficient predict and absence of activity (see also Fig. S2, Supporting Information). In fact, the low friction coefficient implies an optimum value of stress loading where the pore-fluid component takes an important role in stress transfer and earthquake triggering. The earthquake triggering caused by the 1994 seismic event results from 0.1 to 1.0 bar ΔCFF at the 2016 and 2004 receiver faults, respectively (Fig. 3a and Table 3). The ΔCFF modeling shows that the stress transfer due to the 1994 earthquake promotes the 2004 earthquake failure (Figs 2a and 3 and Table 3, and Fig. S1a, Supporting Information), and both 1994 and 2004 earthquakes promote the failure of the 2016 earthquake (Figs 2a and b and Table 3, and Fig. S1b, Supporting Information). The cumulative post-seismic deformation due to the elastic dislocation increases the stress loading on the 2016 rupture from 0.24 to 0.3 bars. ROLE OF PORE FLUID IN THE EARTHQUAKE SEQUENCE A stress change may result from pore-fluid diffusion. If the stress field satisfies the strain compatibility equation $$\frac{{{\partial ^2}}}{{\partial x_j^2}}[ {\frac{{2( {{\nu _u} - \nu } )}}{{B( {1 - \nu } )( {1 + {\nu _u}} )}}\ P + \sigma } ]$$ = 0 (Rice & Cleary 1976; see also appendix A for more details) and if we consider the boundary condition (pore pressure is neglected far from the fault), a simple solution is given by Bosl & Nur (2002):   $$\frac{\sigma }{{{\sigma _{{\rm{init}}}}}} = \ \frac{{\left( {{\nu _u} - \nu } \right)}}{{\left( {1 - \nu } \right)\left( {1 + {\nu _u}} \right)}}$$ (6b)where σ is the change in the stress field due to the pore-fluid diffusion, σinit is the initial stress induced by the coseismic dislocation and νu and ν are the undrained and drained Poisson ratios, respectively. This relation also shows the correspondence between the post-seismic mean stress change induced by pore pressure relaxation and the mean stress caused by the initial dislocation. Taking into account the Rice & Cleary (1976) and Bosl & Nur (2002) solutions, the short-term poroelastic deformation is defined as a diffusive process and can be interpreted as a linear combination of pore pressure and mean stress changes. The drained and undrained Poisson ratios used in the coupled poroelastic stress modeling are 0.25 and 0.31, respectively (Fig. 3b); these values are typical in unconsolidated sedimentary aquifers and water saturated rocks in the upper few kilometres of the seismogenic zone (h ≤ 15-km depth). Our coupled poroelastic modeling suggest a value of short-term post-seismic stress equal to 0.6 bar due to the coupled poroelastic effect for μ΄ = 0.4 related to a an internal friction coefficient of μ = 0.75 and Skempton coefficient B = 0.47. The μ΄ < 0.4 effective friction coefficient gives a value of μ = 0.75 and Skempton coefficient B = 0.9. The variation of B at the intersection of the two ruptures to the end of the 2004 rupture zone can be interpreted as a variation of porosity and the diffusive process along the fault zone (Scholz 1990). The same phenomenon is observed in reservoir induced seismicity where the elastic and the coupled poroelastic effects are considered instantaneous (Scholz 1990). The diffusive process from the reservoir is associated with the fluid migration at short time delay, and moving from a region with a high B to a region with a low B, indicating a short-time poroelastic rebound where the fluid flow transfer from the 1994 and 2004 rupture zones can be compared as reservoir induced seismic activity. Taking into account the time of occurrence of aftershocks and the short-term post-seismic stress change induced by the main shock using the Bosl & Nur hypothesis (2002), we note that the stress load occurring 10–100 d after the 2004 main shock is far more important than the stress load due to the several years of post-seismic deformation computed by InSAR time-series (PS and SBAS, Cetin 2015). The comparable computations take into account the 2004 rupture as source fault and 2016 rupture as receiver fault, and the post-seismic deformation added as a cumulative moment is incorporated into the elastic dislocation modeling. Due to the absence of pore pressure in-situ data and difficulty to perform a 3-D model of stress change and related poroelastic dislocation, we use eqs (2) and (6) to evaluate the stress and pore pressure changes related to the Al Hoceima earthquake sequence. The apparent constant friction model (eq. 2) and variable (or isotropic) friction model (eq. 6 and in appendix B) are considered as parameters able to improve our knowledge on the pore-fluid effect in the Al Hoceima-Alboran region. Nevertheless, the isotropic poroelastic model appears to be the most appropriate for modeling the short-term diffusive process in a complex fault zone. The low value of μ΄ (≤0.4) at the intersection between the 1994 and the 2004 fault (Fig. 4) implies that the first shock reduces the pore pressure along the 1994 rupture, while it increases along the 2004 rupture. The low value of μ΄ implies a high values of stress change due to the short-term coupled poroelastic effect (Fig. 3b). Fig. 4 shows how the pore pressure can affect the triggered seismicity when hydrological processes are coupled to the rupture process. The large pore pressure change associated with the volumetric strain corresponds to fluid migration close to and from the 1994 rupture to the 2004 rupture zone (Fig. 4). Das & Scholz (1981) suggest that the stress effect on the fault might be enhanced as P (see eq. 1) and restored on the main fault, causing fluids to migrate into the receiver fault. The same observation is also made by Jonsson et al. (2003) for the Mw 6.5 2000 June earthquake of SW Iceland showing an increasing water level change associated with high pore pressure change in unconsolidated sedimentary aquifers. The evaluation of pore pressure change at depth in Fig. 4 is resolved for a homogeneous elastic half-space using the theory of linear elasticity (Rice & Cleary 1976); eq. (5) (see also Appendix A5) is used to evaluate the pore pressure change by computing the mean stress change due to the 1994 earthquake for a constant Skempton coefficient. Figure 4. View largeDownload slide Calculated pore-pressure change based on the coseismic volumetric strain and theory of linear poroelasticity (Rice & Cleary 1976) following the 1994 earthquake. The 1994 and the 2004 rupture are represented by black lines (see also Fig. 2d). The 1994 coseismic slip create a high pore pressure zone in the rupture nucleation zone of the next 2004 earthquake, according to Terzaghi (1925) definition of effective stress, the increase in pore pressure diminishes the normal stress acting on the fault and promote the 2004 failure. Figure 4. View largeDownload slide Calculated pore-pressure change based on the coseismic volumetric strain and theory of linear poroelasticity (Rice & Cleary 1976) following the 1994 earthquake. The 1994 and the 2004 rupture are represented by black lines (see also Fig. 2d). The 1994 coseismic slip create a high pore pressure zone in the rupture nucleation zone of the next 2004 earthquake, according to Terzaghi (1925) definition of effective stress, the increase in pore pressure diminishes the normal stress acting on the fault and promote the 2004 failure. The successive earthquakes in the Al Hoceima region may be a response to stress loading in the Rif Mountains and related pore-fluid diffusion in the upper crust. The coupling between crustal deformation and pore-fluid effect imply that the pore pressure may decrease the rock strength by reducing the effective stress and creating a slip-instability that favours earthquake triggering. The seismicity rate change following the 2004 earthquake shows additional aftershocks (see relative rate fluctuations in Fig. 5a), which cannot be explained as aftershock rate decrease as predicted by the Omori law. As the afterslip along seismic ruptures and related post-seismic deformation (Cetin 2015) may have enhanced the pore-fluid flow and contributed to the earthquake triggering in the Rif-Alboran region, in our case, we include a coupled poroelastic component in the ΔCFF modeling and analyse the post-seismic stress transfer and related effective friction coefficient (μ΄ ≤ 0.4). The ΔCFF due to the full poroelastic relaxation is simply computed from a dislocation in a homogeneous elastic half-space using both undrained and drained Poisson's ratios and obtain the difference as seen in Fig. S3 in the Supporting Information. Figure 5. View largeDownload slide (a) The seismic frequency and relative seismicity rate following the 2004 main shock (blue line) and cumulative number of seismic events (green line). The fluctuation in the seismicity rate change shows additional aftershocks possibly due to pore-fluid diffusion in the upper crust. The relative rate change are obtained from changes in slop of the cumulative number curve using a Habermann function regardless to the time of greatest change and comparing the rate in the two parts of the period (before and after the division point) by appropriate time windows function (Wyss & Habermann 1988; Wyss & Viemer 2000), the time variation function defines the local time variation between the rate before and after. (b) The seismicity rate change versus time in the Al Hoceima region. We show the complexity of aftershocks sequences as a realistic representation of the temporal post-seismic effect. Figure 5. View largeDownload slide (a) The seismic frequency and relative seismicity rate following the 2004 main shock (blue line) and cumulative number of seismic events (green line). The fluctuation in the seismicity rate change shows additional aftershocks possibly due to pore-fluid diffusion in the upper crust. The relative rate change are obtained from changes in slop of the cumulative number curve using a Habermann function regardless to the time of greatest change and comparing the rate in the two parts of the period (before and after the division point) by appropriate time windows function (Wyss & Habermann 1988; Wyss & Viemer 2000), the time variation function defines the local time variation between the rate before and after. (b) The seismicity rate change versus time in the Al Hoceima region. We show the complexity of aftershocks sequences as a realistic representation of the temporal post-seismic effect. The pore-fluid flow in drained condition should reduce the normal stress σn (eq. 2) and favour significant additional aftershocks able to affect the relative seismicity rate change (Nur & Booker 1972; Cocco & Rice 2003). In our case, the deep-seated water flow in the flysch units of the fold-and-thrust Rif Belt and substratum metamorphic complex (Chalouan et al.2008) can be considered as another mechanism responsible for the seismicity increase. Piombo et al. (2005) suggest that for significant earthquakes with Mw  > 6, the stress transfer may occur under a fluid diffusion within a 10–20 km radius. The temporal evolution of seismicity within individual fractures takes 10–12 yr to travel up to 20 km (5–6 km yr−1) in our case as well as in other case studies (e.g. Pytharouli et al.2011). Taking into account the complexity of aftershock sequences, a realistic representation of the temporal post-seismic factor, the modified Omori Law (Utsu 1969) or Omori-Utsu Law (Narteau et al.2009) can be expressed as:   $$\lambda \left( t \right) = \ \frac{k}{{{{\left( {t + c} \right)}^{ - x}}}}$$ (7)where λ is the aftershock frequency within a given magnitude range, t is the time from the main shock triggered event, k is the productivity of aftershocks that depends on the total number of events, x is the power law exponent and c is the time delay before the onset of the power-law aftershock decay rate and is dependent on the rate of activity in the earlier part of the seismic sequence. The change of c values characterizes the aftershocks sequence and can be correlated with the stress field orientation (Narteau et al.2009), Guo & Ogata (1997) obtain a range of c value between 0.003 and 0.3 d for various earthquake datasets. In our case, the c value has to be the lowest possible and is fixed as 0.01 d, in order to obtain sufficient aftershock productivity. Note that the c value is often retained connected to the incompleteness of seismic catalogues soon after strong earthquakes. For aftershocks, the seismicity decay requires a time-dependent process that is much faster than the large-scale tectonic loading and much slower than the propagation of elastic waves (Nur & Booker 1972). In our application, we observe that immediately after the 2004 main shock, the aftershock rate decays by $$\frac{1}{{\sqrt t }}$$ while it becomes equal to $$\frac{1}{t}$$ in subsequent months. It appears that the decay rate is due to fluid flow in the crust by means of a diffusion process that contributes to the aftershock sequence (Fig. 5b). Similar results are obtained by Shapiro et al. (1997) for the events related to pressure changes in operation wells and recently by Turkaya et al. (2015) in laboratory experiments. THE CLOCK-TIME ADVANCE AND PERIODIC FREQUENCY The Al Hoceima earthquake sequence shows about 10 and 12 yr recurrence with an aftershock distribution and stress loading that correlate with the location of earthquake ruptures (Figs 2d and Figs S2a and S3, Supporting Information). The mechanism controlling the dependence time of earthquake ruptures, aftershocks and related stress change is complex. The well-resolved 2004 aftershocks distribution during a short time interval (100 d; Tahayt et al.2009) confirms our observation that a significant seismicity rate change is observed 50 d after the main shocks (see Fig. 5a); the positive change in seismicity is correlated with the positive aftershocks productivity due to pore-fluid diffusivity (Fig. 5b). To study the influence of the coseismic stress change on recurrence time interval, Chéry et al. (2001) point that a positive shear stress change on a fault plane should advance the time of the next earthquake on this fault. Here, we consider that the nucleation will occur where the ΔCFF has a maximum value (Console et al.2010). The clock-time advance and related recurrence time (in years) depends on the stressing rate in the positive ΔCFF expressed by the linear equation (Stein et al.1997):   $${\rm{Tr'}} = {\rm{Tr}}-\Delta {\rm{CFF}}/\dot{\tau },$$ (8)where Tr΄ is the calculated recurrence time, $$\Delta t = \Delta {\rm{CFF}}/\dot{\tau }$$ is the clock-time advance ($$\dot{\tau }$$ is stressing rate) and Tr is the mean recurrence time before the earthquake. The stressing rate is computed from the strain rate which is derived from the seismic moment following the equation: $$\dot{M} = 2\mu \Sigma W\dot{\varepsilon }/k$$, where μ is the shear modulus, Σ is the surface area of the region, W is the seismogenic thickness, $$\dot{\varepsilon}$$is the strain rate and parameter k is a dimensionless constant that adjusts for the inefficiency of randomly oriented faults to accommodate strain. Note that Kostrov (1974) chooses k = 1, while Anderson (1979) chooses k = 0.75, and using their relations, our calculated strain rate gives 19 and 24 nanostrain yr−1, respectively (see Table 4). These strain rate values are comparable to the Negredo et al. (2002) and Palano et al. (2013) results and stressing rate of 461 and 582 Pa yr−1, respectively. Hence, for the 1.1 bar maximum stress change obtained from the Coulomb modeling on the 2004 receiver rupture (Table 3), the clock-time advance Δt for the 2004 receiver source able to generate a significant earthquake with Mw  ≥ 6 is 239 ± 22 yr using the Kostrov (1974) relation and 189 ± 17 yr using the Anderson (1979) relation. Table 4. Clock-time advance (Δt) associated with the ΔCFF and corresponding strain rate in the vicinity of the 1994 source rupture, and for the 2004 and 2016 receiver faults in the Al Hoceima region. k = 0.75 and 1 are the values used by Kostrov (1974) and Anderson (1979), respectively, to compute the strain rate and the conditional probability Pc is calculated over 10 yr for earthquakes with Mw  > 6. Seismic  Strain rate  $$\dot{\tau }$$ stress rate  ΔCFFmax    Pc no stress  Pc with stress  ruptures  (nanostrain yr−1)  (Pa yr−1)  (MPa)  $$\Delta t = \frac{{\Delta {\rm{CFF}}}}{{\tau '}}({\rm{yr}})$$  change (pre-1994)  change (post-1994)  2004  19 (k = 0.75)  461  0.11  239 ± 22  12 per cent  55 per cent  2004  24 (k = 1)  582  0.11  189 ± 17      2016  19 (k = 0.75)  461  0.02  21 ± 02      2016  24 (k = 1)  582  0.02  17 ± 02      Seismic  Strain rate  $$\dot{\tau }$$ stress rate  ΔCFFmax    Pc no stress  Pc with stress  ruptures  (nanostrain yr−1)  (Pa yr−1)  (MPa)  $$\Delta t = \frac{{\Delta {\rm{CFF}}}}{{\tau '}}({\rm{yr}})$$  change (pre-1994)  change (post-1994)  2004  19 (k = 0.75)  461  0.11  239 ± 22  12 per cent  55 per cent  2004  24 (k = 1)  582  0.11  189 ± 17      2016  19 (k = 0.75)  461  0.02  21 ± 02      2016  24 (k = 1)  582  0.02  17 ± 02      View Large From the seismicity catalogue, the conditional probability for a specified time interval depends only on the time interval between large earthquakes ΔT, and the long-term regional seismicity rate. The conditional probability for earthquake triggering with Mw  > 6 for the Al Hoceima region is given by (Cornell 1968):   $${\rm {Pc}} = 1-{{\rm e}^{ - \lambda \Delta {\rm{T}}}}$$ (9)where λ is the seismicity rate with magnitude M > 6, and ΔT is the elapsed time since the most recent large earthquake (M > 6) obtained from the Moroccan seismicity catalogue (Jabour, personal communication, 2014). The value of λ is obtained from the Gutenberg–Richter (G-R) law constrained by the parameters a and b. Here, we also observe that λ is different before and after the 1994 stress perturbation (λ ranges between 0.012 and 0.09 following the 1994 earthquake). Based on Cornell & Winterstein (1988) hypothesis, the conditional probability related to seismic recurrence models appears to be the most appropriate (i) where the hazard is dominated by the nearest fault segment, and (ii) in the absence of slip rate and strain rate on each fault segment. We show in Table 4 that after the 1994 earthquake, the 2004 earthquake fault is under a high value of clock-time advance, for example, Δt = 239 ± 22 yr or Δt = 189 ± 17 yr according to Kostrov (1974) and Anderson (1979) formula's, respectively. However, the conditional probability for a specified time interval of 10 yr to have an earthquake with Mw  > 6 rises from 12 per cent to 55 per cent(from eq. 9, see also Table 4), in agreement with the seismic rate activity. The change in pore-fluid pressure along ruptures induced by successive earthquakes results in a cluster of large seismic events (Mw  > 6) during a short period of time (∼22 yr). The conditional probability Pc is similar in the 2004 and the 2016 seismogenic areas, since we consider the regional probability condition. In his study of the 1992 Landers California earthquake (Mw 7.3), Hardebeck (2004) shows no significant difference in the uncertainty between the conditional probability obtained from stress change and that based on the G-R distribution. In our case, Pc depends only on the G-R seismicity temporal distribution, while other studies use different approaches based mainly on stress drop. DISCUSSION AND CONCLUSIONS A sequence of three earthquakes occurred in the Rif Mountains and nearby Alboran Sea in 1994 (Mw 6.0), 2004 (Mw 6.4) and 2016 (Mw 6.5). The static stress change modeling (ΔCFF in undrained condition) suggest a fault-rupture interaction with stress loading located on the selected receiver faults. The poroelastic relaxation (drained condition) and the coupled short-term poroelastic stress transfer help us to understand the seismic migration induced by the pore-fluid diffusion. Aftershock sequences of the three earthquakes correlate well with the ΔCFF distribution which confirms the role of pore fluid in the triggering of post-1994 earthquakes (Figs S2 and S3, Supporting Information). Besides the fault-rupture complexity, the modeling parameters require two levels of friction coefficient with μ΄ < 0.4 correlating with significant pore-fluid diffusivity, and μ΄ = 0.4 in fault zones with limited diffusivity (Fig. 3b). The stress change and related pore-fluid diffusion may explain the ∼10–12 yr interval and the seismic sequence migration. The role of strain rate and its impact on stress change and pore-fluid diffusion combined with the permeability along fault-rupture zones in the Al Hoceima region is crucial in the comprehension of the time delay between the earthquake sequences. The 1994 and the 2004 earthquake ruptures illustrate the stress level change, related value of friction coefficient and role of pore-fluid diffusion in conjugate fault geometry. Although the accuracy in location of the 1994 and 2016 aftershocks is limited due to the azimuthal gap and absence of near-field seismic stations, the distribution of seismic sequences concurs with the ΔCFF and cumulative loading areas (Fig. S3, Supporting Information). However, the incomplete seismicity catalogue with precise earthquake locations in the Al Hoceima region prevents a suitable study on the role of fluid pore pressure before the 1994 earthquake sequence. The earthquake fault locations inland retrieved from the InSAR analysis of coseismic and after slip surface deformation agree with the aftershock distribution. The limited distance between earthquake ruptures (<25 km) and fault geometries with strike-slip mechanisms also promote the stress transfer and failure on fixed fault planes. Earthquake faults in the Al Hoceima region are blind with basically no geomorphological signature at the surface (Tahayt et al.2009). Therefore, fault parameters such as slip-per-event; long-term active deformation and slip rate are missing in our study. The static strain release by the 1994 event induced a high pore pressure change with fluid flow on the 2004 rupture area. This hydrological phenomenon affects the fault zone permeability and promotes the failure of the 2004 event. Wang (2000) uses the correlation between fluid migration and rock permeability to explain the link between the two phenomena; he points out that if the pore pressure becomes too large, earthquakes occur and will increase permeability with groundwater fluid flow. For a strain rate ranging between 19 and 24 nanostrain yr−1 (Table 4), we observe that the regional aftershock frequency following the 2004 earthquake (Fig. 5b) is in good agreement with the simulated aftershock frequency based on the pore-fluid diffusion hypothesis (Bosl & Nur 2002). Pore-fluid effects comparable to our case study at the intersection between the 1994 and 2004 ruptures (Fig. 4) are also observed for conjugate earthquake ruptures during the Superstition Hills earthquakes (Hudnut et al.1989; Scholz 1990). In addition, the decay in the 2004 aftershock activity includes variable seismicity rate probably due to the pore-fluid diffusion. A similar behaviour is observed during the 1966 Parkfield-Cholame earthquake where aftershock productivity and related fluid pore pressure have a direct effect on rock strength (Nur & Booker 1972). Depending on the geological structures, substratum permeability and seismicity rate, the 2004 and 2016 earthquakes could have been predictable by the Coulomb modeling taking into account the pore-fluid effect (undrained and drained conditions). In fact, the occurrence of the 2016 January 21 Mw 5.0 foreshock (4 d before the main shock) may have allowed the fluids to migrate across the epicentral area promoting the 2016 earthquake rupture. Comparable phenomenon with foreshocks and fluid migration across a fault zone is described for the L’Aquila earthquake sequence (Lucente et al.2010). The timescale of post-stress redistribution for the 1994, 2004 and 2016 Al Hoceima earthquakes is larger, for instance, than the ∼11 hr Superstitions Hills sequences (Mw 6.2 and 6.6), suggesting different diffusion processes probably controlled by the permeability along the fault zone. Following the 1994 earthquake sequence, the probability for triggering an Mw  > 6 earthquake within 10 yr interval increases to 55 per cent with respect to the 12 per cent pre-1994 period (Table 4). With the computed 239 ± 22 yr clock-time advance for large earthquakes (Mw  > 6) on the 2004 rupture, the seismic strain rate and ΔCFF explains the 10–12 yr delay and the 55 per cent probability of promoting failure in the Rif Mountains. Acknowledgements JK benefits from a scholarship from the Algerian Ministry of Higher Education and Scientific Research (MESRS). The authors wish to thank Nacer Jabour (CNRST Rabat) and Abdelilah Tahayt (University Mohamed V, Rabat) for the access to the seismicity database. We are thankful to Zoe Shipton and Mark Stillings for their comments on an early version of the manuscript, and to Ross Stein and Shinji Toda for sharing the Coulomb 3.4 software. We are grateful to two anonymous reviewers for their contributions in improving an early version of the manuscript. This research programme was funded by the Direction of Research at MESRS and the DERCI-CNRS with INSU-UMR 7516 IPG Strasbourg. Some figures were prepared using the public domain GMT software (Wessel & Smith 1998). REFERENCES Akoglu A.M., Cakir Z., Meghraoui M., Belabbes S., El Alami S.O., Ergintav S., Akyüz H.S., 2006. The 1994–2004 Al Hoceima (Morocco) earthquake sequence: conjugate fault ruptures deduced from InSAR, Earth planet. Sci. Lett. , 252, 467– 480. Google Scholar CrossRef Search ADS   Anderson J.G., 1979. Estimating the seismicity from geological structure for seismic-risk studies, Bull. seism. Soc. Am. , 69, 135– 158 Beeler N.M., Simpson R.W., Hickman S.H., Lockner D.A., 2000. Pore fluid pressure, apparent friction, and Coulomb failure, J. geophys. Res. , 105( B11), 25 533–25 542. Google Scholar CrossRef Search ADS   Bezzeghoud M., Buforn E., 1999. Source parameters of the 1992 Melilla (Spain, Mw = 4.8), 1994 Al Hoceima (Morocco Mw = 5. 8) and Mascara (Algeria, Mw 5.7) earthquakes and seismotectonic implications, Bull. seism. Soc. Am. , 89, 359– 372. Bosl, W.J., Nur A., 2002. Aftershocks and pore fluid diffusion following the 1992 Landers earthquake, J. geophys. Res. , 107( B12), 2366, doi:10.1029/2001JB000155. Google Scholar CrossRef Search ADS   Cakir Z., Meghraoui M., Akoglu A.M., Jabour N., Belabbes S., Ait-Brahim L., 2006. Surface deformation associated with the Mw 6.4, February 24, 2004 Al Hoceima (Morocco) earthquake deduced from InSAR: implications for the active tectonics along North Africa, Bull. seism. Soc. Am. , 96, 1– 10. Google Scholar CrossRef Search ADS   Calvert A., Gómez F., Seber D., Barazangi M., Jabour N., Ibenbrahim A., Demnati A., 1997. An integrated geophysical investigation of recent seismicity in the Al-Hoceima region of North Morocco, Bull. seism. Soc. Am. , 87, 637– 651. Cetin E., 2015. Analysis and modeling of crustal deformation using InSAR time series along selected active faults within the Africa-Eurasia convergence zone, PhD thesis , University of Strasbourg, EOST, 135 pp. Chalouan A., Michard A., El Kadiri Kh., Negro F., Frizon de Lamotte D., Soto J.I., Saddiqi O., 2008. The Rif Belt, in Continental Evolution: The Geology of Morocco , pp. 203– 302, eds Michard A.et al.  , Springer. Google Scholar CrossRef Search ADS   Cherkaoui T.-E., Hatzfeld D., Jebli H., Medina F., Caillot V., 1990. Etude microsismique de la région d'Al Hoceima, Bull. Inst. Sci. Rabat , 14, 25– 34. Chéry J., Merkel S., Bouissou S., 2001. A physical basis for time clustering of large earthquakes, Bull. seism. Soc. Am. , 91( 6), 1685– 1693. Google Scholar CrossRef Search ADS   Cocco M., Rice J.R., 2003. Pore pressure and poro-elasticity effects in Coulomb stress analysis of earthquake interactions, J. geophys. Res. , 107, doi:10.1029/2000JB000138. Console R., Murru M., Falcone G., 2010. Probability gains of an epidemic-type aftershock sequence model in retrospective forecasting of M ≥ 5 earth-quakes in Italy, J. Seismol. , 14, 9– 26. Google Scholar CrossRef Search ADS   Cornell C.A., 1968. Engineering seismic risk analysis, Bull. seism. Soc. Am. , 58, 1583– 1606. Cornell C.A., Winterstein S.R., 1988. Temporal and magnitude dependence in earthquake recurrence models, Bull. seism. Soc. Am. , 78( 4), 1522– 1573. Das S., Scholz C.H., 1981. Theory of time-dependent rupture in the Earth, J. geophys. Res. , 86( B7), 6039–6051. Google Scholar CrossRef Search ADS   Demets C., Gordon R.G., Argus D.F., 2010. Geologically current plate motions, Geophys. J. Int. , 181, 1– 80. Google Scholar CrossRef Search ADS   El Alami S.O., Tadili B.A., Cherkaoui T.E., Medina F., Ramdani M., Ait Brahim L., Harnaff M., 1998. The Al Hoceima earthquake of May 26, 1994 and its aftershocks: a seismotectonic study, Ann. Geofis. , 41, 519– 537. Fernandez-Ibanez F., Soto J.I., Zoback M.D., Morales J., 2007. Present-day stress field in the Gibraltar Arc (western Mediterranean), J. geophys. Res. , 112, B08404, doi:10.1029/2006JB004683. Google Scholar CrossRef Search ADS   Grevemeyer I., Gràcia E., Villaseñor A., Leuchters W., Watts A.B., 2015. Seismicity and active tectonics in the Alboran Sea, Western Mediterranean: constraints from an offshore-onshore seismological network and swath bathymetry data, J. geophys. Res. , 120, doi:10.1002/2015JB012073. Guo Z., Ogata Y., 1997. Statistical relations between the parameters of aftershocks in time, space, and magnitude, J. geophys. Res. , 102( B2), 2857–2873. Google Scholar CrossRef Search ADS   Hardebeck J.L., 2004. Stress triggering and earthquake probability estimates, J. geophys. Res. , 109, B04310, doi:10.1029/2003JB002437. Hatzfeld D., Frogneux M., Girardin D., 1977. Etude de la sismicite dans la region de I'arc de Gibraltar, Bull. Soc. geol. Fr. , S7-XIX( 4), 741– 747. Google Scholar CrossRef Search ADS   Hudnut W., Seeber L., Pacheco J., 1989. Cross-fault triggering in the November 1987 Superstition Hills earthquake sequence, Southern California, Geophys. Res. Lett. , 16( 2), 199– 202. Google Scholar CrossRef Search ADS   Jonsson S., Segall P., Pedersen R., Bjornsson G., 2003. Post-earthquake ground movements correlated to pore-pressure transients, Nature , 424, 179– 183. Google Scholar CrossRef Search ADS PubMed  King, G.C.P., Stein, R.S., Lin, J., 1994. Static stress changes and the triggering of earthquakes, Bull. seism. Soc. Am. , 84( 3), 935– 953. Kirkpatrick J.D., Shipton Z.K., Evans J.P., Micklethwaite S., Lim S.J., McKillop P., 2008. Strike-slip fault terminations at seismogenic depths: the structure and kinematics of the Glacier Lakes fault, Sierra Nevada United States, J. geophys. Res. , 113, B04304, doi:10.1029/2007JB005311. Google Scholar CrossRef Search ADS   Kostrov V.V., 1974. Seismic moment and energy of earthquakes, and seismic flow of rocks, Izv. Acad. Sci. USSR Phys. Solid Earth , 1, 23–44 (Eng. Transl.). Koulali A.et al.  , 2011. New GPS constraints on active deformation along the Africa-Iberia plate boundary, Earth planet. Sci. Lett. , 308, 211– 217. Google Scholar CrossRef Search ADS   Lucente F.P., Gori P.D., Margheriti L., Piccinini D., Bona M.D., Chiarabba C., Agostinetti N.P., 2010. Temporal variation of seismic velocity and anisotropy before the 2009 Mw 6.3 L’Aquila earthquake, Italy, Geology , 38, 1015– 1018. Google Scholar CrossRef Search ADS   Medina F., 1995. Present-day state of stress in northern Morocco from focal mechanism analysis, J. Struct. Geol. , 17( 7), 1035– 1046. Google Scholar CrossRef Search ADS   Meghraoui M., Pondrelli S., 2012. Active faulting and transpression tectonics along the plate boundary in North Africa, Ann. Geophys. , 55( 5), 955– 967. Morel J.-L., Meghraoui M., 1996. The Goringe–Alboran–Tell tectonic zone, a transpression system along the Africa–Eurasia plate boundary, Geology , 24, 755– 758. Google Scholar CrossRef Search ADS   Narteau C., Byrdina S., Shebalin P., Schorlemmer D., 2009. Common dependence on stress for the two fundamental laws of statistical seismology, Nature , 462( 7273), 642– 645. Google Scholar CrossRef Search ADS PubMed  Negredo A.M., Bird P., de Galdeano C.-S., Buforn E., 2002. Noetectonic modeling of the Ibero-Maghrebian region, J. geophys. Res. , 107( B11), ETG10-1– ETG10-15. Google Scholar CrossRef Search ADS   Nur A., Booker J.R., 1972. Aftershocks caused by pore fluid flow?, Science , 175( 4024), 885– 887. Google Scholar CrossRef Search ADS PubMed  Okada Y., 1992. Internal deformation due to shear and tensile faults in a half-space, Bull. seism. Soc. Am. , 82( 2), 1018– 1040. Palano M., González P.J., Fernández J.A., 2013. Strain and stress fields along the Gibraltar orogenic arc: constraints on active geodynamics, Gondwana Res. , 23, 1071– 1088. Google Scholar CrossRef Search ADS   Piombo A., Martinelli G., Dragoni M. 2005. Post-seismic fluid flow and Coulomb stress changes in a poro-elastic medium, Geophys. J. Int. , 162, 507– 515. Google Scholar CrossRef Search ADS   Pytharouli S.I., Lunn R.J., Shipton Z.K., Kirkpatrick J.D., do Nascimento A.F., 2011. Microseismicity illuminates open fractures in the shallow crust, Geophys. Res. Lett. , 38, L02402, doi:10.1029/2010GL045875. Google Scholar CrossRef Search ADS   Reasenberg P.A., Simpson R.W., 1992. Response of regional seismicity to the static stress change produced by the Loma Prieta earthquake, Science , 255, 1687– 1690. Google Scholar CrossRef Search ADS PubMed  Rice J.R., Cleary M.P., 1976. Some basic stress diffusion solutions for fluid-saturated elastic porous media with compressible constituents, Rev. Geophys. Space Phys. , 14( 2), 227– 241. Google Scholar CrossRef Search ADS   Shapiro S.A., Huenges E., Borm G., 1997. Estimating the crust permeability from fluid-injection-induced seismic emission at the KTB site, Geophys. J. Int. , 131, F15– F18. Google Scholar CrossRef Search ADS   Scholz C.H., 1990. The Mechanics of Earthquakes and Faulting , Cambridge Univ. Press, p. 439. Stein R.S., Barka A., Dieterich J.H., 1997. Progressive failure on the North Anatolian fault since 1939 by earthquake stress triggering, Geophys. J. Int. , 128, 594– 604. Google Scholar CrossRef Search ADS   Stich D., Serpelloni E., Mancilla F., Morales J., 2006. Kinematics of the Iberia-Maghreb plate contact from seismic moment tensors and GPS observations, Tectonophysics , 426, 295– 317. Google Scholar CrossRef Search ADS   Stich D., Martin R., Morales J., 2010. Moment tensor inversion for Iberia-Maghreb earthquakes 2005–2008, Tectonophysics , 483, 390– 398. Google Scholar CrossRef Search ADS   Tahayt A. et al., 2009. The Al Hoceima (Morocco) earthquake of 24 February (2004) analysis and interpretation of data from ENVISAT ASAR and SPOT5 validated by ground-based observations, Remote Sens. Environ. , 113, 306– 316. Google Scholar CrossRef Search ADS   Terzaghi K., 1925. Erdbaumechanik auf Bodenphysikalischer Grundlage. Franz Deuticke , Liepzig-Vienna. Timoulali Y., Hahou Y., Jabour N., Merrouch R., El Kharrim A., 2014. Main features of the deep structure by local earthquake tomography and active tectonics: case of Rif Mountain (Morocco) and Betic Cordillera (Spain), J. Seismol. , 18( 2), 221– 234. Google Scholar CrossRef Search ADS   Toda S., Stein R.S., Sevilgen V., Lin J., 2011. Coulomb 3.3 graphic-rich deformation and stress-change software for earthquake, tectonic, and volcano research and teaching, User guide, U.S. Geol. Surv . Open-File Rept. 2011-1060, 63 pp. Available at: http://pubs.usgs.gov/of/2011/1060/, last accessed October 2016. Turkaya S., Toussaint R., Ericksen F.K., Zecevic M., Daniel G., Flekkøy E.G., Måløy M., 2015. Bridging aero-fracture evolution with the characteristics of the acoustic emissions in a porous medium, Front. Phys. , 3, doi:10.3389/fphy.2015.00070. Utsu T., 1969. Aftershocks and earthquake statistics (I) - Some parameters which characterize an aftershock sequence and their interrelations, J. Fac. Sci. Hokkaido Univ., Ser. , VII( 3), 121– 195. Vallée M., 2016. Automatic determination of source parameters using the SCARDEC method. Available at: http://geoscope.ipgp.fr/index.php/en/catalog/earthquake-description?seis=us10004gy9, last accessed 22 September 2016. Wang H.F., 2000. Theory of Linear Poroelasticity with Applications to Geomechanics and Hydrogeology , pp. 304, Princeton Univ. Press. Wessel P., Smith W.H.F., 1998. New, improved version of generic mapping tools released, EOS, Trans. Am. geophys. Un. , 79( 47), 579, doi:10.1029/98EO00426. Google Scholar CrossRef Search ADS   Wyss M., Habermann R.E., 1988. Precursory seismic quiescence, Pure appl. Geophys. , 126, 319– 332. Google Scholar CrossRef Search ADS   Wyss M., Wiemer S., 2000. Change in the probabilities for earthquakes in Southern California due to the Landers M 7.3 earthquake, Science , 290, 1334– 1338. Google Scholar CrossRef Search ADS PubMed  SUPPORTING INFORMATION Supplementary data are available at GJI online. Figure S1. (a) ΔCFF modeling on the 2004 fault plane due to the 1994 fault rupture, the computation are made for μ΄ = 0.2, the black stars indicate the 2004 hypocentre, (b) ΔCFF modeling on the 2016 fault plane due to the cumulative 1994 and 2004 earthquake, the grey stars indicate the 2016 hypocentre, note that the receiver planes are tapered to 40 fault patches system. Figure S2. ΔCFF modeling on optimally oriented fault planes caused by the 2004 earthquake at depth range [5–15 km] for: (a) μ΄ = 0.2, (b) μ΄ = 0.4, (c) μ΄ = 0.6 and (d) μ΄ = 0.8, the aftershocks distributions are represented by black circles, the white star represents the 2004 epicentre and the grey star is the 2016 epicentre. The increase of pressure could trigger events in regions where the Coulomb stress predict and absence of activity. Figure S3. Full poroelastic relaxation of the 1994 and the 2016 earthquakes from the undrained state to the fully drained state, the computation is obtained by linear elasticity theory with appropriate values of undrained and drained Poisson ratio: (a) full poroelastic relaxation due to the 1994 earthquake on the 2004 fault rupture with a typical sedimentary undrained and drained Poisson ratios, respectively, (b) full poroelastic relaxation due to the 1994 earthquake on the 2004 fault rupture with an extreme undrained and drained Poisson ratios, respectively, according the Bosl and Nur hypothesis, this value might be typical for upper crust when fractures exists, (c) ΔCFF ( red + green) and ΔP/B profiles along the 2004 rupture related to the full poroelastic relaxation of the 1994 earthquake, a significant change of ΔP/B is obtained for (vu, v) = (0.31,0.15), the yellow stars indicate the 2004 epicentre location on the axe, (d) full poroelastic relaxation due to the 2016 earthquake on fixed planes strike/dip/rake = 195°/78°/19° for (vu, v) = (0.31,0.15) and (e)full poroelastic relaxation due to the 2016 earthquake on fixed planes strike/dip/rake = 195°/78°/19° for (vu, v) = (0.31,0.25). Note that the white stars indicate the 2016 epicentre location. Figure S4. Angular relationship between pre-seismic local stress field and the 2004 fault plane acting on the central Rif block, the maximum horizontal stress (green) and the minimum horizontal (yellow) stress directions acting of the 2004 rupture at depth are obtained from stress estimation based on inversion of focal mechanism and GPS data obtained by several authors (Akoglu et al.2006; Fernandez-Ibanez et al.2007; Tahayt et al.2009), σ and τ denote the normal and shear stress according to the state of stress, the red black and the red line indicates the SS surface fault of the 1994 and 2004 earthquakes, respectively, the tectonic models used for this study are from Tahayt et al. (2009), where arrows indicate the relative movements of these blocks with respect to Africa plate, the 1994 and the 2004 rupture are from Akoglu et al. (2006). The majors reverse faults trace are represented by an appropriate symbol and the major strike-slip fault traces represented by lines. Table S1. Coulomb stress change caused by the 2004 and the 2016 earthquakes on aftershocks planes for μ΄ = 0.4, the strike/dip/rake represent the best geometry related to the optimally oriented ΔCFF loading. SF represents the source fault. Appendix A: Constitutive equations of poro-elasticity. Appendix B: Relationship between the poro-elastic model and the 2004 earthquake triggering. Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper. © The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Jan 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations