Geodetic Deformation versus Seismic Crustal Moment-Rates: Insights from the Ibero-Maghrebian Region

: Seismic and geodetic moment-rate comparisons can reveal regions with unexpected potential seismic hazards. We performed such a comparison for the Southeastern Iberia—Maghreb region. Located at the western Mediterranean border along the Eurasia–Nubia plate convergence, the region has been subject to a number of large earthquakes (M ≥ 6.5) in the last millennium. To this end, on the basis of available geological, tectonic, and seismological data, we divided the study area into twenty-ﬁve seismogenic source zones. Many of these seismogenic source zones, comprising the Western Betics, the Western Rif mountains, and the High, Middle, and Saharan Atlas, are characterized by seismic / geodetic ratio values lower than 23%, evidencing their prevailing aseismic behavior. Intermediate seismic / geodetic ratio values (between 35% and 60%) have been observed for some zones belonging to the Eastern Betics, the central Rif, and the Middle Atlas, indicating how crustal seismicity accounts only for a moderate fraction of the total deformation-rate budget. High seismic / geodetic ratio values ( > 95%) have been observed along the Tell Atlas, highlighting a fully seismic deformation.


Introduction
Seismic and geodetic moment-rate comparisons provide relevant insights into the seismic hazard of regions subjected to significant crustal deformation. Valuable examples come from numerous tectonic regions worldwide, such as Iran [1,2], western Canada [3], western USA [4,5], Greece [6,7], southern Italy [8][9][10], Himalayas [11][12][13], and the East African Rift [14]. Geodetic moment-rate represents a measure of crustal deformation and might include both anelastic and elastic components, whereas seismic moment-rates measure the deformation accommodated by faulting. Although seismic and geodetic moment-rate estimations are affected by several physical uncertainties, their comparison has enabled identifying regions where the total deformation-rate budget is entirely released by crustal seismicity [3,4,10], as well as regions where the excess deformation-rate can be released either in aseismic slip across faults or through large future earthquakes [1,2,5,13,14]. Regarding this last point, a common practice is to estimate the number of large earthquakes (usually called "missing earthquakes") needed so that the moment released by seismicity balances the deficit of moment derived from geodetic measurements [3,4,13].   [24,25]. The map is plotted in an oblique Mercator projection.
A number of on-and off-shore active fault systems dissect the study area and are responsible for a large amount of the current seismic release of the region. Among them, the Trans Alboran shear zone (a major NE-SW trending left-lateral strike-slip system) cuts the Alboran basin, and branches into the Eastern Betics shear zone and into the Nekor/Al-Hoceima seismic zone, northward and southward, respectively (e.g., [26][27][28]). In addition, the Trans Alboran shear zone is westward connected to southern Spanish off-shore by the right-lateral Yusuf fault [29]. In western Algeria, active faulting occurs within the Tell Atlas, along NE-SW-trending, right-stepping en-echelon reverse faults [30]. In eastern Algeria, active reverse faulting occurs in a broader area [31][32][33] and is coupled with right-lateral strike slip faulting on EW-trending faults [34]. Active faulting is recognized along the Aurès Mountains and the southern Atlas region [35].
Southern Iberia and Maghreb regions have been subject by the occurrence of large earthquakes (with estimated magnitude M ≥ 6.5) in the last centuries [36][37][38]. Global Navigation Satellite System (GNSS)-based measurements show that active deformation involves both Southern Iberia and Maghreb regions with rates up to 6 mm/yr (e.g., [23,39]). Such a GNSS-based active deformation has been interpreted in terms of elastic block modeling (e.g., [39,40]. However, the detection of a significant aseismic deformation component on the Betic-Rif system (∼75%; see [41] for details) highlights that the comparison of seismic and geodetic deformation-rates would not be balanced across the region. Results achieved in this study support the inference on the aseismic deformation behavior for most of the Betics and Rif regions, revealing a fully seismic deformation along most of the Tell Atlas.

Historical Earthquakes
Available historical seismic catalogs for southern Iberia and Maghreb regions highlight the occurrence of large (M ≥ 6.5) earthquakes in Tunisia since AD 856 [36][37][38]42], including not well-documented damaging earthquakes since several centuries before Christ (BC). However, in spite of efforts made in recent updates [25] focused on seismic hazard studies, the accuracy of these catalogs is not uniform due both to the off-shore occurrence of numerous earthquakes and to the sparsely Remote Sens. 2020, 12, 952 4 of 27 populated area. A number of studies focusing on this region have evidenced that many small and moderate earthquakes have most likely been overlooked, and that available seismic catalogs can be considered complete for M ≥ 6 earthquakes since approximately 1500, 1700, and 1850, for southern Iberia, northern Morocco and northern Algeria areas, respectively (e.g., [25,37,38] and references therein). Historical earthquakes are largely concentrated along the Azores-Gibraltar fault zone, along the Betic-Rif system, and along the Tell Atlas fold-thrust belt (Figure 2a). In detail, the strongest historical earthquakes (M ≥ 7.0) are mainly concentrated on the Azores-Gibraltar fault system (AD 1309, 1356, 1761, and 1755 events) and along the Tell Atlas fold-thrust belt (AD 856, 1365(AD 856, , 1716(AD 856, , 1722(AD 856, , 1790(AD 856, , 1832(AD 856, , 1867, and 1891 events).

Instrumental Seismicity
The distribution of instrumental seismicity (earthquakes with M ≥ 2.5) as well as the main seismotectonic features (focal plane solutions for crustal events with M ≥ 3.0) of the study area are reported in Figure 2b,c, respectively. The bulk of instrumental seismicity, corresponding to low and moderate magnitude earthquakes, is spread across the Nubia-Eurasia plate boundary. It shows clear epicenter alignments between Azores and Gibraltar, and from Algeria to Tunisia, providing evidence of active faulting along the plate boundary belt. Along the Azores-Gibraltar fault zone and the Gulf of Cádiz, earthquakes have depths up to 40 km, being characterized by fault plane solutions with both reverse and strike-slip features (Figure 2b,c). Along North Algeria and Tunisia regions, earthquakes usually have hypocentral depths shallower than 30 km and are characterized by fault plane solutions with prevailing reverse features, coupled with subordinate strike-slip faulting features (Figure 2b,c). Crustal seismicity spreads over a~300 km wide band between Northern Morocco and Southern Spain, where a number of epicenter alignments suggests the presence of individual active seismogenic faults; therefore, depicting a complex network of NE-SW and WNW-ESE faults (with lengths from tens to some hundreds of kilometers) that accommodate the Nubia-Eurasia convergence (Figure 2b). Seismicity at intermediate depths concentrates in the 50-100 km depth interval, in a 150-km-long band that extends from Southern Spain to the Alboran Sea and abruptly disappears toward east. In addition, the occurrence of deep earthquakes down to more than 600 km depth below the Granada Basin has been documented since 1954 (see [43,44] and references therein for additional details). These deep earthquakes are likely related to a sunk lithospheric slab whose relation with the current surface deformation is still debated (e.g., [45]).
Eurasia convergence (Figure 2b). Seismicity at intermediate depths concentrates in the 50-100 km depth interval, in a 150-km-long band that extends from Southern Spain to the Alboran Sea and abruptly disappears toward east. In addition, the occurrence of deep earthquakes down to more than 600 km depth below the Granada Basin has been documented since 1954 (see [43,44] and references therein for additional details). These deep earthquakes are likely related to a sunk lithospheric slab whose relation with the current surface deformation is still debated (e.g., [45]).   [46]. Symbols are colored according to focal depth. (c) Lower hemisphere, equal area projection for fault plane solutions (with M ≥ 3) compiled from the investigated area [23,[47][48][49]; fault plane solutions are colored according to rake: red indicates pure thrust faulting, blue is pure normal faulting, and yellow is strike-slip faulting.

Seismogenic Source Zones
The study area has been divided into twenty-five crustal seismogenic source zones ( Figure 1) by taking into account the seismic zonation available in recent literature [24,25,50], which have been defined on the basis of the latest tectonic, geological, and seismological considerations. Some of these zones are well-defined on the basis of constraints coming from the analysis of the earthquake catalog (stationarity of the completeness periods, evaluation of the mean activity rate) and from a set of geological and seismotectonic considerations, such as style, geometry, and distribution of fault systems (with direct evidence of Holocene activity) and their relation to the local stress and deformation regimes. Other zones, although defined on similar seismotectonic or geological information, are lacking a clear correspondence between the contemporary seismic activity and the Holocene tectonic activity.
The seismic sources adopted for the Moroccan-Algeria-Tunisian region come from [24] and slightly updated in [50]. For this region, the general criterion adopted in [24,50] was to delimit areas with similar tectonic features and a more or less homogeneous stress pattern and seismicity parameters. Individual faults were not specifically considered because of the very scarce availability of kinematic, seismic and paleoseismic observations in this region. They were only considered ensuring that known faults or fault segments were not divided and partially included in different adjacent sources.
The seismic sources adopted for southeastern Spain belong to the last official upgrade of the Spanish seismic hazard maps [25], known as the 'ByA12' zonation (Zonificación de Bernal y asociados). This zonation is prevailing based on the available seismotectonic and geological information and was aimed at seismic hazard computations [25]. Such a zonation is characterized by some small zones sampled by zero or just one GNSS site; therefore, we joined some of them into large zones following more generalized seismic and tectonic criteria ( Figure 1). We performed such an additional step because the computed strain-rates cannot be taken as an accurate estimate for small zones, due to the lack of spatial constraints from geodetic data (i.e., the strain-rates information would be derived from the velocity interpolation only).
A short description of the twenty-five seismogenic source zones is given below (see also Table A1 in Appendix A for additional details).

•
Tell seismogenic source zones (T1, T2, T3, T4, T5, and T6). Tell region is the one with the largest density of earthquakes, including some of the most destructive events in the area, such as the 10 October, 1980 Chlef (M W 7.3) and the more recent 21 May, 2003 Boumerdes (M W 6.9) earthquakes. The differentiation among these six zones is mainly related to the different seismic features, clearly revealed when computing seismicity parameters. These zones are based on the previous works by [51][52][53]. T1 to T4 sources display a dominant compression striking NNW-SSE, whereas the T5 source exhibits a N-S-striking compressional strike-slip regime (e.g., [24,49,54]. Few focal mechanisms are available for the T6 source and therefore the mean stress regime is poorly constrained (Figure 2c).

•
Rif seismogenic sources (R1a, R1b, and R2). The Rif region has been divided into three source zones on the basis of well-known active faults (e.g., the Jebah and Nekor faults) and seismic activity distribution in the Rif and in the southwestern part of the Alboran Sea. Large earthquakes, such as the recent 25 January, 2016 (M w 6.3) Al Hoceima event or the historical 1624 May, 11 (M S 6.7) Fez event, occurred along some of these active faults. R1b and R2 sources show a horizontal compression, striking NW-SE to NNW-SSE, as well as a perpendicular horizontal extension [50]. There are few mechanisms for R1a source, but they show a pure WSW-ENE horizontal extension [54]. Jerada events. The MA-HP source is based on a previous one already described in [52,53]. The MM source is characterized by the presence of NNE-SSW and NE-SW oriented faults, as well as faults parallel to the Atlantic coast [57]. The stress field for both seismogenic source zones is poorly constrained because of the limited number of available focal mechanisms (Figure 2c).

•
Betic seismic source zones (BET1, BET2, BET3, BET4, BET5, and BET6). These seismogenic sources correspond to the seismogenic region called Betic Cordilleras, also embracing the Subbetic, western part of the Prebetic units, and the Guadalquivir Basin, in its northern part. With the exception of the 29 March, 1954 (M w 7.8) deep seismic event (depth of~625 km), both instrumental and historical earthquakes did not exceed M 6.3 [24]. As mentioned above, initial sources delimited in [25] have been grouped into larger seismogenic sources. BET1 source zone embraces a complex network of faults striking WSW-ENE on the Guadalquivir Basin and NE-SW to NNE-SSW on the Sierra de Segura area, which hosted some historical destructive earthquakes, for instance the 1357 (M w 6.1) Andalucia and the 1504 (M w 6.1) Carmona events. BET2 source zone is characterized by faults striking N140/160E and N70E to E-W, as well as faults striking NW-SE which hosted destructive earthquakes, such as the 1431 (M w 6.11) Granada event and the 1884 (M w 6.30) Arenas del Rey event. BET3 source zone is characterized by low to moderate instrumental seismicity with a M w 6.17 destructive event occurring in 1531 close to the town of Baza. BET4 source zone is characterized by low to moderate instrumental seismicity with a Mw 6.10 destructive event occurring in 1910 close to the small coastal town of Adra. BET5 source zone includes, from west to east, the Sierra Nevada, the Sierra de los Filabres, and the Almanzora Corridor areas. This zone is characterized by NW-SE to NNW-SSE oriented faults which are considered the seismogenic sources for seismicity striking the region. BET6 source zone includes the Adra-Sierra Alhamilla and Almería areas, incorporating to some extent the off-shore platform. In this source, existing fault systems mostly strike NW-SE and NNE-SSW, hosting historical destructive earthquakes, such as the 1522 (M w 6.07) Almería event. Regarding the stress regime, BET2 and BET4 source zones are characterized by a primary NNW-SSE sub-horizontal compression, while BET1, BET3, BET5, and BET6 zones are characterized by a complex pattern resulting from the interaction of different sources of stress (see [54,55] for additional details). • LEV1 and LEV2 seismogenic source zones. This area, named as Levante seismogenic region [25], covers the easternmost part of the Prebetic and Internal Betics and is dominated by a prevailing strike slip stress field with a near E-W extension and a N-S compression [54]. Both seismogenic source zones are characterized by the occurrence of several low to moderate magnitude earthquakes and some damaging historical events, rarely exceeding M > 6 such as the 1396 (M W 6.00) Tavernes de Valldigna and the 1645 (M S 6.25) Muro de Alcoy for LEV1 and the 1829 (M W 6.25) Torrevieja event for LEV2 [36]. Among the main tectonic features included in this source, the Cádiz-Alicante and the Caudete-Elda-Elche fault systems are worth noting. Regarding the LEV2 seismogenic source zone, different fault systems host the seismicity, such as the striking N140-160E San Miguel de Salinas and Torrevieja faults and the NW-SE Bajo Segura fault.

Seismic Moment-Rates
For each seismogenic source zone, the seismic moment-rate was calculated by adopting a truncated cumulative Gutenberg-Richter distribution [16]: where M max is the magnitude of the largest earthquake that could occur within the investigated region; ϕ is a correction for the stochastic magnitude-moment relation; taking into account an average error of 0.2 on magnitudes, we assumed ϕ = 1.27 [16]. c and d are the coefficients of the magnitude (M)-scalar moment (M seis ) relation: Remote Sens. 2020, 12, 952 8 of 27 According to [58] we set c = 1.5 and d = 9.05. a and b are coefficients (seismicity rate and slope, respectively) of the Gutenberg-Richter recurrence relation: being N M the cumulative number of earthquakes with magnitude M and larger. Table 1. Earthquake catalog statistic and estimated seismic moment-rates ( . M seis ) with associated errors for each seismic source zone (SSZ). Moment-rates have been estimated by using both a (#) truncated cumulative Gutenberg-Richter distribution [16] and the (*) the moment summation [6] approach. Seismicity rate (a) and slope (b) of the Gutenberg-Richter recurrence relation with its associated error are reported, as well as the values of the maximum magnitude observed (M obs ) and estimated (M max ). Abbreviations are as follows: P, [37]; H, [38], S, [36]. Because of the large uncertainties of MA-HP and SA1, we also reported the recomputed a and b values (see the main text for details); in curved brackets the original values. Adopted parameters are detailed in Table 1; as a first step, we used the values already reported in literature [24,25]. As described in [24,25], a and b parameters have been directly computed from available catalogs by adopting the maximum likelihood method [59]. The M max values and associated uncertainties have been estimated by adopting different approaches. Regarding the Spanish seismogenic source zones, M max values have been estimated by adopting empirical laws [60] linking magnitude values to the dimensions of mapped active faults. For the zones characterized by a lack tectonic data, M max values have been estimated directly from the seismic catalog by adopting a cumulative exponential-truncated Gutenberg-Richter recurrence relation. Regarding the seismic source zones located along the African margin, M max values have been estimated from the seismic catalog by adopting the Kijko-Sellevoll approach [61], a method usually applied on regions where only a limited number of large earthquakes is available. Because the seismogenic source zones MA-HP and SA1 are characterized by b-values affected by large uncertainties (0.41 and 0.26, respectively), we performed new estimations by using as input event records reported in the International Seismological Centre (ISC) bulletin (see Figure A2 in Appendix B for details). Finally, estimated seismic moment-rates with an associated 67% confidence interval are reported in Table 1.

SSZ
The scalar seismic moment-rate can be estimated also by adopting the moment summation approach [15]: where N is the number of events occurring during a given time interval ∆T in the volume A·H s (with A, the surface area and H s , the seismogenic thickness), M seis is the scalar seismic moment of the n-th earthquake from the N total earthquakes. In the following, we focused mainly on the results obtained from Equation (1); however, in order to test sensitivity, we considered results achieved by adopting Equation (4) (see Table 1).

GNSS Data Processing
Here we analyse an extensive GNSS dataset which, covering 20 years of observations (from 1999.00 up to 2019.00), includes more than 300 continuous GNSS sites available at EUREF Permanent GNSS Network [62], at Crustal Dynamics Data Information System [63], at UNAVCO [64] and from local and regional networks [23]. Time series of this GNSS dataset cover different time spans, ranging from 3.5 to 20 years with an average duration of 8.7 years. We also included 25 episodic GNSS sites located in Morocco with measurements carried out during the 1999.80-2006.71-time interval [40]. The GNNS phase observations were processed by using the GAMIT/GLOBK 10.7 software [65] following the strategy described in [66]. As a final processing step, a consistent set of positions and velocities in the International Terrestrial Reference Frame ITRF2014 reference frame [67] has been computed.
In order to improve the spatial density of the geodetic velocity field over the studied area, we integrated our solutions with those reported in [39,[68][69][70]. In particular, we aligned published solutions to our ITRF2014 one by solving for Helmert transformation parameters that minimize the root mean square (RMS) of differences between velocities at common sites. The resulting velocity field, aligned to the ITRF2014, is reported in Table S1. Such a velocity field is characterized by an average site spatial density of 3.5 × 10 −4 site km 2 (i.e., ca. 11 sites within a radius of 100 km), with differences between the regions (Figure 3). The Betics have a high density of sites (8.3 × 10 −4 site km 2 ), the Rif region has a density of sites of 3.0 × 10 −4 site km 2 , while the Tell and the Saharan regions have a density of sites of 0.9 × 10 −5 site km 2 . In order to highlight regional and crustal deformation patterns over the investigated area, we align our ITRF2014 velocities to both a Eurasian and Nubian reference frames (see Table S1). The estimated velocity fields are reported in Figure 3.  Table S1 in the Supplementary Material section). Continuous and episodic GNSS stations are reported as yellow and red arrows, respectively. c) Geodetic strain-rate field: red and blue arrows represent the greatest extensional (εHmax) and contractional (εhmin) horizontal strain-rates, respectively; the colors in the background show the rotation strain-rate (which describes the rotations with respect to a downward positive vertical axis with clockwise (cw) positive and counterclockwise (ccw) negative, as in normal geological conventions). The map is plotted in an oblique Mercator projection. Table 2 lists the estimated geodetic moment-rates for all the identified source zones along with the associated uncertainties given from the ones related to the strain-rates, the seismogenic thicknesses and the shear modulus (for these last two parameters we assumed an uncertainty of 10% and 5% with respect to the chosen value, respectively).  Table S1 in the Supplementary Material section). Continuous and episodic GNSS stations are reported as yellow and red arrows, respectively. (c) Geodetic strain-rate field: red and blue arrows represent the greatest extensional (ε Hmax ) and contractional (ε hmin ) horizontal strain-rates, respectively; the colors in the background show the rotation strain-rate (which describes the rotations with respect to a downward positive vertical axis with clockwise (cw) positive and counterclockwise (ccw) negative, as in normal geological conventions). The map is plotted in an oblique Mercator projection.
Looking at the velocity field referred to a Eurasian reference frame (Figure 3a), the primary pattern is given by a crustal shortening of~5 mm/yr related to the Nubia-Eurasia oblique convergence. Such a crustal shortening is largely adsorbed on the Alboran and the Balearic-Algerian basins along off-shore structures. Eastern Betics are characterized by geodetic velocities showing a NW-to-NE fan-shaped pattern with rates ranging from~3 mm/yr near the coast to~0.8 mm/yr inland (Figure 3). To the east, the Baleares promontory is characterized by a motion of~0.5 mm/yr toward East. Stations located on the southern sector of the Betics move toward SW while stations located in the eastern Rif move toward NW; therefore, highlighting a NNW-SSE to N-S contraction of the Alboran Basin. Stations located on the central sector of the Gibraltar Arc are characterized by a westward motion, leading to an E-W stretching of the western side of the Alboran Basin. Moreover, the Gibraltar Arc is affected by an E-W contraction of~0.3 mm/yr, as deduced by differentiating the motion between the stations located externally and internally to the arc (Figure 3a). On northern Algeria along the coastal area, geodetic velocities decrease progressively from west (~3.0 mm/yr) to east (~0.5 mm/yr). Southward, the region encompassing the Middle, High, and Saharan Atlas as well as the western High Plateaus, is characterized by velocities similar to those predicted by the rigid Eurasia-Nubia plate motion [23,55], suggesting that the region belongs to stable Nubia. This last feature is also readily visible in the Nubian reference frame (Figure 3b), which shows small geodetic velocities (~0.3 mm/yr) in the Middle, High, Saharan Atlas and in the western High Plateaus, and velocities with values up to 5 mm/yr along the coastal area of the Tell Atlas and the Aurès Mountains.

Strain-Rate Field Estimation
The horizontal strain-rates have been estimated on a regular 0.5 • × 0.5 • grid over the investigated area by adopting the method reported in [71]. This method allows introducing different spatial weighting functions of data (e.g., uniform Gaussian or quadratic spatial weighting function), enabling to obtain a finer resolution, especially on regions characterized by sparsely distributed data. Based on some preliminary tests (see the Appendix C for additional details), the horizontal strain-rate field has been estimated by adopting a weighting threshold of 24 and by using a Gaussian function for distance weighting and Voronoi cell for areal weighting, respectively.
The estimated horizontal strain-rates and the rotation-rates are reported in Figure 3b. A large region comprising the eastern part of the Saharan Atlas, the Aurès Mountains and the Balearic-Algerian basin, is characterized by a prevailing shortening of~10 nanostrain/yr along the NW-SE orientation. On the Algerian-Tunisian side, this shortening is coupled with clockwise rotation-rates up to 15 nanoradian/yr, suggesting the occurrence of transpressive right-lateral faulting between the Aurès Mountains and the contiguous eastern Saharan Atlas/High Plateau block. The High Atlas and the western part of the Saharan Atlas are characterized by an elongation of~5 nanostrain/yr along a prevailing E-W direction. Westward, the regions comprising the High Atlas and the Middle Atlas are characterized by a shortening of~5 nanostrain/yr along the N-S orientation, passing to a NW-SE orientation along the westernmost sector of the Rif. All these regions are characterized by small rotation-rates (~2 nanoradian/yr), highlighting their low tectonic activity. The region comprising the eastern Rif, the Alboran Basin and the Betics are characterized by a shortening of~12 nanostrain/yr along the NW-SE orientation coupled with an elongation up to~16 nanostrain/yr along the NE-SW orientation. This region is characterized by a complex rotation-rate pattern: clockwise rotation-rates up to 20 nanoradian/yr can be observed along the north-western Rif and the western Betics, while counterclockwise rotation-rates of~15 nanoradian/yr and~10 nanoradian/yr can be observed along the eastern Betics and the eastern Rif, respectively.

Geodetic Moment-Rates Computation
The geodetic moment-rate (and its associated uncertainty) was estimated for each seismic source zone with volume A·H s as [72]: where µ is the shear modulus of the crust, ε Hmax and ε hmin are the principal horizontal strain-rates, and Max is a function returning the highest value of the arguments. Since the geodetic moment-rate  [73] and references therein). Our estimated H s value (18 km) for the Tell Atlas region falls within the range of 15-20 km proposed in recent literature (e.g., [31,74,75]). Therefore, in the following computations we adopted the values of 15 and 18 km as estimated from the depth distributions of earthquakes for the Betics and Tell regions, respectively (see Table 2). Regarding the Rif-western Atlas regions, [41] assumed a fixed H s value of 15 km for a wide area including the Rif region, while [76] observed that~85% of seismicity located in Central High Atlas has focal depths lesser than 20 km. Based on these considerations, and since the completeness of the catalog will vary considerably throughout the period of time studied and throughout the study region, we cautiously have chosen an H s value of 20 km for the Rif-western Atlas region. M geod ) with associated error and adopted parameter values for each seismic source zone (SSZ). A, area of the SSZ; µ, shear modulus for crustal rocks; H s , seismogenic thickness with a 10% variation for the lower and upper ranges (see Appendix B for additional details). Regarding the shear modulus µ, a value of 30 × 10 10 N/m 2 has been used, typical of average crustal rocks [77], and allowed for a 5% variation for the lower and upper ranges.  Table 2 lists the estimated geodetic moment-rates for all the identified source zones along with the associated uncertainties given from the ones related to the strain-rates, the seismogenic thicknesses and the shear modulus (for these last two parameters we assumed an uncertainty of 10% and 5% with respect to the chosen value, respectively).

Results and Discussions
Here we provide a detailed statistical evaluation of the seismic/geodetic moment-rates ratio for the Ibero-Maghrebian region by comparing GNSS-based moment-rates with those derived from seismic catalogs.

Moment-Rates
The estimated seismic moment-rates are reported in Table 1 and Figure 4. As above mentioned, the estimations have been performed by adopting two different approaches: the truncated Gutenberg-Richter distribution ( [16]; Figure 4a) and the moment summation ( [15]; Figure 4b).
Remote Sens. 2020, 12, 952 13 of 27 suffers from the possible lack of large earthquakes with high recurrence rate compared to the catalog duration. Based on these observations, in the following we refer to the results coming from the truncated Gutenberg-Richter distribution approach.  Table 1 and Table 2 for details). SSZ characterized by moment-rates larger than 7.  Estimations coming from the former approach, with the exception of very few zones (e.g., HA, SA2, and R1a; see Figure A4 in Appendix D), are generally larger (~40% on average) than the ones estimated from the latter approach; however, both estimations overlap within their associated uncertainties. The truncated Gutenberg-Richter distribution estimates the magnitude distribution of the earthquakes through the computation of the b value (Equation (3)), assuming that earthquake distribution is ergodic over the investigated region and represents an adequate sample of long-term seismicity, therefore allowing to take into account the probable incompleteness of the existing catalog (e.g., [78]). The moment summation method is strongly influenced by the catalog completeness and suffers from the possible lack of large earthquakes with high recurrence rate compared to the catalog duration. Based on these observations, in the following we refer to the results coming from the truncated Gutenberg-Richter distribution approach. The estimated geodetic moment-rates and associated uncertainties are reported in Table 2; computed values range in the interval 4. zone. The geodetic moment-rate computed for the TA seismogenic source zone is based on strain-rates information derived from the velocity interpolation only, therefore this estimation cannot be taken as an accurate estimate, due to the lack of spatial constraints from geodetic data (TA is not covered by GNSS stations). In addition, the geodetic moment-rate inferred for zone HA is statistically negligible because of the large uncertainties associated with the estimated value. Therefore, based on these considerations, the TA and HA seismogenic source zones have been excluded from the comparison reported below.

Seismic/Geodetic Moment-Rates Ratio
In the following, we define the seismic/geodetic moment-rates ratio (expressed as a percentage) as the "Seismic coupling coefficient" (hereinafter SCC; [2]), the geodetic moment-rate being a measure of both elastic and anelastic loading rates and the seismic moment-rate a measure of the elastic unloading rate. The theoretical range of SCC is between 0% and 100%. A SCC value close to 100% indicates that a large amount of the deformation budget is released through earthquakes. Conversely, a low ratio indicates an apparent seismic moment deficit, suggesting either a proportion of aseismic deformation (i.e., ongoing unloading by creep and other plastic process) or accumulating strain not released by seismicity (i.e., elastic storage).
Indeed, SCC can exceed the value of 100%, i.e., the seismic moment-rate is larger than the geodetic ones. This is commonly observed when seismic catalogs cover a very short time interval and contain one or more large earthquakes [2]. The temporal length of a seismic catalog as well as its completeness governs the adequacy of catalogs to estimate seismic moment-rates over a given region (e.g., [78]). Roughly speaking, a seismic catalog with a short duration (100-300 years) could be little adequate to capture the seismic cycle of a given region; indeed, seismic moment-rates estimated from seismic catalogs require that the average earthquake recurrence interval should be shorter than the catalog duration (see [2,4] for additional details). This means that the used seismic catalog would be long enough to capture a complete earthquake cycle for an individual seismogenic fault or alternatively, to provide a statistical sample of all phases of the seismic cycle, including of course earthquakes, for a region containing multiple seismogenic faults. Moreover, 50-100 year-long instrumental catalogs are the most common source of data used to derive earthquake statistics worldwide, under the assumption that such a time span is adequate to derive earthquake return periods over timescales of 500-5000 years, typically used in probabilistic seismic hazard analysis [3]. As mentioned above, available historical seismic catalogs for southern Iberia and Maghreb regions highlight the occurrence of large (M ≥ 6.5) earthquakes since AD 856 and can be considered complete for M ≥ 6 earthquakes since approximately 1500, 1700 and 1850, for the Southern Iberia, northern Morocco and northern Algeria areas, respectively [25,37,38]. Being aware that on the one hand, the earthquake occurrences and statistics may not be steady state over the whole timescales of our catalog and that on the other hand, the truncated Gutenberg-Richter distribution allows taking into account the probable incompleteness of the catalog [78], in the following we carefully discuss and frame our results in comparison with available and independent geological and geophysical observations.
Regarding geodetic moment-rate estimations, the GNSS measurements should sample a large spatial scale, so that they are not affected by local strain accumulation on individual faults, and a long enough time interval, so that measurement uncertainties have a minimal effect on velocity estimations and adequately sample both seismic, aseismic and long-term deformation transients. As noted above, our GNSS dataset shows different spatial density over the investigated area, with high density of sites on the Iberian region and low density along the Maghrebian region. Time series of this GNSS dataset cover different sample periods ranging from 3.5 to 20 years (with an average duration of 8.7 years) and, with the exception of a few sites located in the R1b source zone, generally do not include large earthquakes that could significantly contribute with co-seismic and post-seismic displacements to the estimated geodetic velocities. Therefore, our estimated regional velocity and the related strain rate field are statistically significant in most of the study area.
Estimated SCC values are reported in Table 3 and in Figure 5; based on the achieved SCC values, we considered three principal ranges of values: less than 23%, between 35% and 60%, and more than 95%. Table 3. Comparison of the seismic coupling coefficient (expressed as a percentage), which takes into account the moment-rates estimated from the (#) truncated cumulative Gutenberg-Richter distribution [16] and the (*) the moment summation [15] approach. Confidence intervals (67%) are also reported. SCC = seismic coupling coefficient. A large sector of the study area, comprising Betics (BET1, BET3, BET4, BET5, BET6), Rif (R1a, R2), High, Middle, and Saharan Atlas (MM, MA-HP, HA-AA, SA1, SA2), and one region in the Tell Atlas (T5), is characterized by SCC values lower than 23%. This result is highly consistent with previous estimations computed for the Betics and Rif areas (see [41] for details) which highlighted how about 24% of crustal deformation is released seismically. A number of indicators such as surface heat flow, seismic tomography and local rheological models provided evidence of a weak crustal rheology; therefore, supporting the inference on the aseismic behavior of some sectors of the region. In detail, surface heat flow measurements available for the area show high values (100-120 mW m −2 ; [79,80]) along the eastern sector of the Alboran basin and moderate values (60-75 mW m −2 ) along the Betics [81] and Rif [82] regions. In the central and southern Gulf of Cadiz, the surface heat flow is characterized by values close to 40-50 mW m −2 [74], assuming typical values for stable continental/oceanic lithosphere [83]. Local seismic tomographies have inferred a low-velocity zone at 18-km depth beneath the central Betics, indicating variations in lithology and/or in the rigidity of the lower crust rocks (e.g., [84]). Local rheological models (e.g., [85]) highlighted how the crustal yield strength as well as the inferred depth of the brittle-ductile transition follow the curved shape of the Gibraltar arc with maximum depths of 12-9 km, whereas in the Rif and the Betics, such a transition became shallower eastward (~6-5 km depth; Figure 5). All these geological and geophysical evidences support the inference on the aseismic behavior, at least of the Betics and Rif regions.

SSZ
Intermediate SCC values (between 35% and 60%) have been observed for HA-MA, R1b, LEV1, LEV2 and BET2 seismogenic source zones; therefore, suggesting how these regions account only for a moderate seismic fraction of the total deformation-rate budget. Some of these regions are characterized by well-known active faults which frequently generate moderate earthquakes such as LEV2 and R1b. On these zones, a significant contribution to the measured crustal deformation can be attributed to aseismic post-seismic mechanisms such as afterslip and/or viscoelastic relaxation, as for instance documented for earthquakes striking the R1b source zone in the last decades [86].  Table 3). As discussed in the main text, three principal ranges of SCC values have been identified: low (less than 23%), intermediate  Table 3). As discussed in the main text, three principal ranges of SCC values have been identified: low (less than 23%), intermediate (between 35% and 60%) and high (more than 95%). In grey are reported the SSZ excluded from the comparison, i.e., HA and TA (see the text for more details). The white dashed line represents the depth (in km) of the brittle-ductile transition as reported in [85]. (b) Comparison of the seismic coupling coefficient (expressed as a percentage) which takes into account the moment-rates estimated from both the approaches discussed in the text. Confidence intervals (67%) are also reported.
Intermediate SCC values (between 35% and 60%) have been observed for HA-MA, R1b, LEV1, LEV2 and BET2 seismogenic source zones; therefore, suggesting how these regions account only for a moderate seismic fraction of the total deformation-rate budget. Some of these regions are characterized by well-known active faults which frequently generate moderate earthquakes such as LEV2 and R1b. On these zones, a significant contribution to the measured crustal deformation can be attributed to aseismic post-seismic mechanisms such as afterslip and/or viscoelastic relaxation, as for instance documented for earthquakes striking the R1b source zone in the last decades [86].
The higher SCC values (>95%) have been observed along the Tell Atlas zones (T1, T2, T3, T4, T6), suggesting how the measured crustal deformation over these zones is mostly released through earthquakes. All these source zones are characterized by SCC values larger than 100%. These values are not surprising, since, along the Tell Atlas zones, the crustal shortening of~5 mm/yr related to the Nubia-Eurasia oblique convergence is largely adsorbed on off-shore reverse structures bordering the North African margin [23,39], therefore a relevant portion of the long-term elastic strain accumulation is not captured by the on-shore GNSS stations. In addition, a number of moderate to large earthquakes have epicenters concentrated off-shore (not sampled by GNSS stations), such as the 1790 Oran (M~7) [87], the 1856 Djijelli Mw 7.2 [88], and the 2003 Boumerdes Mw 6.8 [89] earthquakes, therefore, their released seismic energy corresponds to strain with no full counterpart on the on-shore surface.
Both these features lead to an underestimation of geodetic strain-rate. As above mentioned, the geodetic moment-rate estimations are affected by the assumed computational parameters; even if increasing the seismogenic thickness value up to 25 km and/or computational parameters (grid size and weighting threshold) for strain-rates estimation, the SCC values for Tell Atlas zones remain confined to values larger than 100%. All these considerations, coupled with the occurrence of large destructive events (M > 7) both in historical and instrumental periods adds realistic constraints on the seismic behavior of the region. In such a context the low SCC value inferred for T5 source zone (~20%) appears puzzling. The seismic moment-rate is based on a maximum magnitude M max of 5.8, about 1 order of magnitude unit lesser than the values estimated for the other Tell Atlas source zones. It must be noted that an increase of 1 magnitude unit leads to an increase by a factor of~5.4 of the moment-rate (Equation (1)), therefore leading to a SCC value close to 100%. These considerations, while on one hand suggest that the low SCC value inferred for T5 source zone could be poorly representative of its seismic behavior, on the other suggest that the possibility of forthcoming earthquakes in the region may increase. Since in this region the occurrence of a moderate earthquake (1985 Constantine Mw 5.8 event) is documented in the available seismic catalogs, it appears to be a potential seismic gap like other well studied regions of the world (e.g., [90,91].
Although the seismic and geodetic moment-rates comparison implies numerous assumptions and simplifications, it provides significant insights into the potential seismic hazard of tectonically active regions, essentially for time-dependent seismic hazard assessments. As noted above, valuable examples come from numerous tectonic regions worldwide such as Iran [1,2], western Canada [3], western USA [4,5], Greece [6,7], southern Italy [8][9][10], Himalaya [11][12][13], and the East African Rift [14]. Despite the different approaches adopted by these authors, all their estimations have enabled achieving valuable results that document cases with a good agreement between seismic and geodetic moment-rates and also cases where geodetic moment-rates are significantly larger than the seismic ones. However, uncertainties related to the physical significance of the deformation-rates mismatch over varying spatial and temporal scales are currently poorly understood. The continuous growth in continuous seismic and GNNS networks is allowing the acquisition of spatially extensive datasets at an increasing number of tectonic areas worldwide, therefore leading to an improved comprehension of such a physical significance.

Conclusive Remarks
Based on the SCC estimations previously presented and discussed, we identified three main regions with different seismic behavior along the western Mediterranean border.

•
A large sector of the study area, comprising the western Betics, the western Rif, and the High, Middle, and Saharan Atlas is characterized by SCC values lower than 23%. Such a result coupled with geological and geophysical evidence supports the inference on the aseismic deformation (aseismic release processes) behavior at least for most of the Betics and Rif regions.

•
Intermediate SCC values (between 35% and 60%) have been observed for some regions belonging to the eastern Betics, the central Rif, and the Middle Atlas. On these regions, crustal seismicity accounts only for a moderate fraction of the total deformation-rate budget.

•
Higher SCC values (>95%) have been observed along the Tell Atlas, highlighting a fully seismic deformation. By considering the low SCC value inferred for T5 source zone, we speculated on the possibility of this being a seismic gap and hosting impending earthquakes.  share a similar pattern, we therefore merged all the seismic source zones into three main regions ( Figure S1a) in order to obtain more robust results. Hence we re-analyzed, for each merged region, the depth distribution-time plots (Figure A1b-d) and the cumulative frequency plots ( Figure A1e-g). Cumulative frequency plots of depth distributions of earthquakes show that 90% of events occur at depths less than 15, 18 and 30 km for Betics (BET1, BET2, BET3, BET4, BET5, BET6, LEV1, and LEV2), Tell Atlas (T2, T3, T4, T5, T6,  The ISC routinely produces catalogs of earthquake hypocenter locations relying on data contributed by seismological agencies from around the world. A large number of earthquakes have their focal depth assigned to 0, 5, 10, 15, 20, 25, and 30 km, because of difficulties in computing reliable depth estimations. All events with M ≥ 3.5 are manually reviewed and relocated by an ISC analyst. The location provided by the ISC is based entirely on P-wave travel-time tables derived from a global 1D Earth velocity model (e.g., the ak135 model [95]). Therefore, many ISC hypocenters are poorly constrained in focal depth, especially in regions not extensively covered by seismic stations, as are the Rif and Atlas ones, and must be interpreted with caution. The ISC routinely produces catalogs of earthquake hypocenter locations relying on data contributed by seismological agencies from around the world. A large number of earthquakes have their focal depth assigned to 0, 5, 10, 15, 20, 25, and 30 km, because of difficulties in computing reliable depth estimations.
All events with M ≥ 3.5 are manually reviewed and relocated by an ISC analyst. The location provided by the ISC is based entirely on P-wave travel-time tables derived from a global 1D Earth velocity model (e.g., the ak135 model [95]). Therefore, many ISC hypocenters are poorly constrained in focal depth, especially in regions not extensively covered by seismic stations, as are the Rif and Atlas ones, and must be interpreted with caution.

Appendix C
A variety of methods, based on different approaches, have been developed to compute strains (see [71] for an overview). Since the spatial distribution of our velocity field data is heterogeneous, in order to verify that the resulting strain-rate field provides a good representation of the regional deformation field of the investigated area, we performed some additional computations by adopting the method described in [71].
Four solutions are computed using a different combination of data weighting functions (12, 18, 24, and 32). The computation was performed on a 0.5° × 0.5° regular grid by adopting a Voronoi cell approach for areal weighting. Achieved results are reported in Figure A3: arrows show the greatest extensional (εHmax) and contractional (εhmin) horizontal strain-rates, while the color in the background shows the second invariant of the strain-rate tensor (defined as = + + 2 , where , and are three strain-rate components in an east and north Cartesian coordinate system). Considering results reported in Figure A3, the strain-rate pattern suggests that as the weighting threshold decreases from 32 ( Figure A3d) to 18 ( Figure A3b), the strain-rate model picks up more tectonic strain signals along main tectonic features within the investigated area. However, as the weighting threshold decreases from 18 ( Figure A3b) to 12 ( Figure A3a), the differential strain-rate pattern starts to highlight very local signals. Some of these local signals are located in regions with dense data population (e.g., south-eastern Spain) while others are located in regions not covered by GNSS data, therefore suggesting a significant deteriorating of the strain-rate pattern at some places.
Balancing the trade-off between the resolution and robustness, we chose the strain-rate pattern resulting from a weighting threshold of 24 as the optimal model for the characterization of the strainrate field in the Ibero-Maghrebian region.

Appendix C
A variety of methods, based on different approaches, have been developed to compute strains (see [71] for an overview). Since the spatial distribution of our velocity field data is heterogeneous, in order to verify that the resulting strain-rate field provides a good representation of the regional deformation field of the investigated area, we performed some additional computations by adopting the method described in [71].
Four solutions are computed using a different combination of data weighting functions (12, 18, 24, and 32). The computation was performed on a 0.5 • × 0.5 • regular grid by adopting a Voronoi cell approach for areal weighting. Achieved results are reported in Figure A3: arrows show the greatest extensional (ε Hmax ) and contractional (ε hmin ) horizontal strain-rates, while the color in the background shows the second invariant of the strain-rate tensor (defined as ε 2inv = ε 2 E + ε 2 N + 2ε 2 EN , where ε E , ε N and ε EN are three strain-rate components in an east and north Cartesian coordinate system).
Considering results reported in Figure A3, the strain-rate pattern suggests that as the weighting threshold decreases from 32 ( Figure A3d) to 18 ( Figure A3b), the strain-rate model picks up more tectonic strain signals along main tectonic features within the investigated area. However, as the weighting threshold decreases from 18 ( Figure A3b) to 12 ( Figure A3a), the differential strain-rate pattern starts to highlight very local signals. Some of these local signals are located in regions with dense data population (e.g., south-eastern Spain) while others are located in regions not covered by GNSS data, therefore suggesting a significant deteriorating of the strain-rate pattern at some places.
Balancing the trade-off between the resolution and robustness, we chose the strain-rate pattern resulting from a weighting threshold of 24 as the optimal model for the characterization of the strain-rate field in the Ibero-Maghrebian region.