Rapid Source Characterization of the 2023 $Mw$ 6.8 Al Haouz, Morocco, Earthquake Open Access

On 8 September 2023, at 11:11 p.m. local time (22:11 UTC), an $Mw$ 6.8 earthquake occurred in the High Atlas Mountains of Morocco, roughly 75 km southwest of Marrakesh, locally referred to as the Al Haouz earthquake after the province most impacted by shaking (Fig. 1). The Al Haouz earthquake resulted in nearly 3000 fatalities and 5500 injuries reported in the weeks following the event (International Blue Crescent Relief and Development Foundation [IBC], 2023). “Did You Feel It?” (Quitoriano and Wald, 2020) reported modified Mercalli intensity (MMI) 9 (violent shaking) near the epicenter, with MMI 7 (very strong shaking) in Marrakesh. The high reported intensities throughout the region correlate with widespread structural damage. The earthquake’s impact was aggravated by the fragility of the building stock and its nighttime occurrence when more people were home.

Although the background seismicity in the region is low, with only nine M 5 and larger earthquakes having occurred within 500 km of this event since 1900, Morocco has experienced damaging earthquakes in the past. Particularly, the 1960 M 5.9 Agadir earthquake sequence (Fig. 1) destroyed the city of Agadir, resulting in 12,000–15,000 fatalities, accounting for roughly 40% of Agadir’s population (Tillotson, 1960; see Data and Resources). The Agadir sequence was devastating primarily due to its proximity to a population center, shallow depth, and the susceptibility to ground shaking of the regions building stock of stone and brick structures (Paradise, 2005). The larger 2023 M 6.8 Al Houz earthquake was less devastating due to its greater source depth and because it was more distant from a large population center (e.g., Marrakesh). Details on the source characteristics of this event (e.g., location, magnitude, mechanism, and rupture kinematics) presented here are critically important to begin developing a better understanding of the seismotectonic environment of the High Atlas Mountains and the potential for the future large earthquakes in the region.

The Al Haouz earthquake ruptured within the Moroccan High Atlas (MHA) Mountains. The MHA is dominated by thick‐skinned thrusting and folding, and the mountains are considered a double‐vergent mountain belt (Lanari et al., 2020; Skikra et al., 2021). Although recognized as an intraplate mountain belt, shortening of the MHA has accommodated a substantial fraction of the total Africa–Eurasia plate convergence since the Miocene, indicating that it is part of a wide plate boundary zone in the western Mediterranean (Gomez et al., 2000) that includes the Betic‐Rif mountain belt to the north. Shortening in the MHA decreases along strike from east to west (Teixell et al., 2003, 2005). The 2023 Al Haouz earthquake occurred in the western segment of the MHA [west of 7° W, which is substantially narrower (30–50 km) than in the central High Atlas (100 km)]. The western MHA is characterized as a transpressive regime for which the MHA is bounded by the North Atlas and South Atlas fault zones (Styron and Pagani, 2020; Fig. 1). Seismic and geodetic modeling shows that the Al Haouz earthquake nucleated in the lower crust (∼24 km) beneath the western MHA. Interferometric Synthetic Aperture Radar (InSAR) observations indicate no surface rupture during the 2023 event. Consequently, no causative fault is readily apparent. It is feasible that rupture occurred on other structures that are exposed at the surface in the High Atlas Mountains, such as the Tizi n’Test fault (e.g., Domènech et al., 2015), or on blind, shallow crustal faults related to thick‐skinned deformation (e.g., El Harfi et al., 2006).

The event hypocenter estimated by the U.S. Geological Survey (USGS) National Earthquake Information Center (NEIC; Table 1, Fig. 1) places the event in the western MHA at a depth of 19 km (Table 1). The NEIC location is primarily controlled by teleseismic observations; the closest station being WM.TIO, roughly 100 km from the epicenter. Only three other stations are within 500 km, with all azimuths between ∼20 and 100°. Therefore, the estimated epicentral location and depth have relatively large uncertainty. Although the location’s formal horizontal uncertainties are ∼3.1 km, these do not account for all sources of uncertainty, for example, location bias introduced by unmodeled velocity structure; and therefore, the epicenter uncertainty is likely much larger. Improving the hypocenter location estimate requires a more accurate crustal velocity model and/or local and regional arrival time data. Consequently, teleseismic waveform modeling and InSAR data provide the best estimates of average slip depth, whereas InSAR data combined with finite‐fault modeling provide constraints on event location and fault kinematics.

Only two aftershocks above NEIC’s routine cataloging threshold (M 4.5) were observed in the week following the event: an M 4.9, which occurred ∼20 min after the mainshock, and an M 4.6 on 14 September 2023 (Fig. 1). NEIC published three additional events: an M 4.0 and two M 4.2 s; but no other events met the NEIC publication criteria, given the lack of arrival time observations. Consequently, aftershock locations were not useful to provide further details on the geometry of the causative fault(s).

$Mww$ (W Phase), $Mwb$ (body wave), and $Mwc$ (centroid/surface wave) moment tensor estimates (Table 1) consistently estimate an $Mw$ 6.8 oblique thrust solution. These solutions closely match the $Mw$ 6.9 solution produced by the Global Centroid Moment Tensor (Global CMT) Project (Dziewonski et al., 1981; Ekström et al., 2012). Waveform characteristics and $Mww$ and $Mwb$ estimates indicate a >90% double‐couple source, suggesting a likely a relatively simple source time function (STF). $Mww$⁠, $Mwb$⁠, and $Mwc$ estimate the depth of the moment centroid to be 30.5 km, 24.0 km, and 26.3 km, respectively, which agree with the Global CMT depth of 27.8 km (Dziewonski et al., 1981; Ekström et al., 2012). NEIC considers $Mwb$ to best constrain depth because of the relatively highly frequencies (12.5–80 s) modeled compared to $Mww$ and $Mwc$ (>100 s), and because depth phases may be observed following the P wave. We further improve the centroid depth estimate, relying on the NEIC SynDepth software (Yeck and Patton, 2021), which fits higher frequency (1.6–50 s) teleseismic P waveforms via a grid search through apparent STF and depth space. We observe clear depth dependence on fits to individual stations (Fig. 2). The mean estimated depth using a stack of 64 selected stations is 25 km ± 0.6 km (Fig. 2d), consistent with moment tensor modeling. We consider this our best event centroid depth, but note that there is some small uncertainty due to unmodeled velocity structures and complexities from variations in local surface elevation. Regardless, modeling unequivocally shows that rupture primarily occurred in the mid‐to‐lower crust.

To provide the best estimate of source location, depth, and slip distribution, we supplemented seismic observations with geodetic measurements of ground displacement derived from InSAR observations. Following the approaches described in Shea and Barnhart (2022), we generate coseismic interferograms using SAR images acquired from the European Space Agency Sentinel‐1a instrument, with postevent images acquired on 11 September (descending path 154), 15 September (ascending path 45), and 16 September (descending path 52). We process the interferograms with the JPL/Caltech/Stanford InSAR Scientific Computing Environment (ISCE) processing package (Agram et al., 2013), downsampling and then inverting the coseismic displacements to estimate a best‐fitting plane with uniform slip (Table 1), which gives constraints on the location and depth, fault geometry, fault dimension, and slip magnitude and direction (Shea and Barnhart, 2022). The interferograms confirm that the earthquake rupture was blind, as indicated by a lack of surface rupture in the surface displacement patterns. Accordingly, we allowed the inversions to search the fault geometry space (strike, dip, and rake) for both north‐ and south‐dipping focal planes.

The interferograms exhibit distinct evidence of a region of uplift (indicated by negative displacements or decreasing distance between the ground and satellite) with a more subtle region of subsidence (indicated by positive displacements or increasing distance between the ground and satellite) south of the uplift lobe (Figs. S1–S3, available in the supplemental material to this article). The southern edge of the uplift lobe is characterized by steeper displacement gradients than the northern edge, which indicates that subsurface slip occurred closer to the southern edge of the uplift lobe. Geometrically, this relationship implies that the causative fault dips to the north. The surface displacement patterns are fit well by both south‐ and north‐dipping focal planes; however, the steep displacement gradient patterns are slightly better fit by the north‐dipping, steeper focal plane. The InSAR‐based inversions estimate the depth of the shallowest slip at ∼14.5 km depth for a north‐dipping focal plane (dip = 67°) and ∼24.0 km for a south‐dipping focal plane (dip = 24°). Both north‐ and south‐dipping focal plane solutions place the source fault in the middle crust of the High Atlas Mountains and south of the seismic epicenter (Table 1).

We combined teleseismic and InSAR observations (Figs. S1–S3) to better constrain the spatiotemporal event slip distribution using the Wavelet and simulated Annealing SliP (WASP) finite‐fault model (FFM) method (Goldberg et al., 2022). The first version of the FFM was completed prior to any InSAR acquisitions and included only teleseismic broadband waveforms with fault geometries constrained by the $Mww$ moment tensor focal planes. Although the two $Mww$ nodal plane orientations fit the data similarly, we initially selected the north‐dipping plane (225°/69°) due to broad alignment with the trend of faulting in the High Atlas Mountains (Styron and Pagani, 2020). This seismically constrained model consisted of 3 km × 3 km subfaults, with 15 subfaults in both the along‐strike and along‐dip directions. Slip was permitted between rake angles of 49° and 89°, corresponding to ±20° from the average rake vector given by the NEIC $Mww$⁠. In general, the slip model of an earthquake of this size is considered poorly constrained with only teleseismic observations. Nonetheless, the model showed realistic slip amplitude and dimensions of slip compared to established scaling laws (Blaser et al., 2010). The model consisted of one slip asperity roughly 25 km × 24 km (along‐strike × along‐dip), with slip concentrated between ∼10 and ∼35 km depth and maximum slip of ∼1.7 m (Fig. 3a).

We updated FFM solutions with each of three subsequent InSAR acquisitions (v.2, 11 September; v.3, 15 September; v.4, 16 September). The introduction of InSAR line‐of‐sight observations indicated a discrepancy between the hypocenter location and the causative fault, which resulted in a spatial offset between the zone of transition from the uplift lobe to the subsidence lobe seen in the observations compared to the inversion results. Given the uncertainty in the hypocenter location, we opted to pin the fault orientation to the location of the InSAR‐inferred centroid, ensuring that the fault plane position would be consistent with the surface displacement observations. This modification means that the centroid is also treated as the initiation point of rupture. Although it may be an imperfect assumption, our observations are insensitive to rupture onset location. The static InSAR observations have no temporal resolution, and the WASP framework includes temporal shifting of teleseismic waveforms to account for imperfect seismic velocity models, so teleseismic waveforms are also insensitive in shifts to the rupture onset time and location. Our final seismogeodetic model (Fig. 3) estimates a slightly smaller source area (20 km × 24 km) with higher maximum slip (2 m) compared to the early teleseismic‐only model—a typical consequence of including regional datasets (e.g., Goldberg et al., 2022). The simple STF lasting ∼12 s (Fig. 3b) in our FFM is similar to the STF produced by the GEOSCOPE SCARDEC method (Vallée et al., 2011).

Although uncertainty exists in each technique leveraged to characterize the earthquake, each technique also gives unique insights into the event. It is clear that slip in this earthquake primarily occurred at depths in the mid‐to‐lower crust, whether one uses waveform modeling of the event as a point source (as the case in $Mww$⁠, $Mwb$⁠, $Mwc$⁠, and SynDepth) or estimates discrete slip as in our FFM. The nodal planes in the various moment tensors are consistent with oblique (strike‐slip + thrust) slip on either a gentle south‐dipping fault or steep north‐dipping fault. Combining these observations with available fault mapping, several causative faults are plausible (Fig. 4). To identify potential causative faults, we use the Global Earthquake Model (GEM) fault database (Styron and Pagani, 2020) fault traces and preferred dip values. We assume constant dip with depth (i.e., no ramp‐flat or listric fault geometries). The two main systems—the North and South Atlas faults—are proposed low‐dip thrust faults (dip ∼ 20°–25°). This idealized dip value is similar to the east‐striking and gently south‐dipping nodal plane in the moment tensor solutions, and would be consistent with the North Atlas fault. However, the North Atlas fault projection is shallower than the current source location (Fig. 4a). Furthermore, this south‐dipping plane is not preferred in the InSAR and FFM inversions. The west‐striking, steeply north‐dipping plane in the moment tensor solutions contrasts with the published geometries of the South Atlas fault (Fig. 4b). The geometry of this nodal plane indicates a steeply dipping fault within the High Atlas range. An alternative to the range‐bounding faults, Tizi n’Test fault, or other unmapped mid‐range faults, are near the surface projection of the north‐dipping fault plane and the surface projection of the FFM preferred plane (Fig. 4c). The Tizi n’Test fault is proposed to be a reactivated normal fault (basin forming fault in the Triassic), with conflicting views of its Pleistocene and younger activity (Sébrier et al., 2006; Delcaillau et al., 2011; Domènech et al., 2015; Fekkak et al., 2018). It is debated if the Tizi n’ Test fault is presently active. Because the subsurface geometry of this fault is unconstrained, we do not project the trace to depth.

Given these constraints, several options remain for the causative fault. A low‐angle detachment underlying the range could have ruptured if the causative fault is the gently‐south‐dipping plane. Such low‐angle faults have been inferred in the region (Hafid et al., 2006; Sébrier et al., 2006). If the causative fault is the north‐dipping plane, the fault could be an unmapped and/or blind structure within the High Atlas range (e.g., Domènech et al., 2015). Understanding that thrust belts rarely have constant dip with depth, it is geologically plausible that the causative fault may have been a ramp (north‐dipping) or a flat (south‐dipping) of a larger range‐bounding structure like the North or South Atlas faults.

Although events of this size are not unexpected in the MHA, the Al Haouz earthquake is noteworthy, because it nucleated in the lower crust and occurred on a blind thrust fault in a region with well‐described surface faults that are a recognized seismic hazard in the western MHA. The recent seismic hazard models of North Africa rely on moment rate contributions from off‐fault (probabilistic) and on‐fault (deterministic) seismic sources. For example, the epicentral region in the GEM North African probabilistic seismic hazard model is defined by a distributed source zone and the North and South Atlas range‐bounding fault systems (Poggi et al., 2020). The distributed source zone in the epicentral region has the maximum magnitude (⁠$Mmax$⁠) of 6.9, derived from the largest earthquake known within the zone + 0.5 magnitude units. $Mmax$ on the range front fault systems ranges from 7.2 to 7.3 (Poggi et al., 2020). Others estimate $Mmax$ from single‐segment faulting to be M ∼6.1–6.4 along the North and South Atlas faults (Sébrier et al., 2006), although they note that the surface faulting may not represent the extent of faulting at depth. Because the centroid represents the average location of slip, this necessitates the fault slipped both above and below this central depth, and therefore it is clear that the mid‐to‐low crust is seismogenic. Given the large uncertainty in deterministic estimates of $Mmax$ and seismogenic depths, it highlights the necessity for probabilistic seismic hazard zones in regions of low seismicity rates and relatively slow‐moving faults.

The combined use of teleseismic waveform data and InSAR data provided accurate estimates of the mainshock’s location, depth, magnitude, mechanism, and kinematic rupture. Modeling of the two closest publicly available stations showed that numerous aftershocks did occur, but they could not be used in our analysis because of a lack of hypocentral certainty. A more complete understanding of the M 6.8 Al Haouz earthquake sequence would benefit from modeling of the local seismic data to better understand the rupture evolution of the fault(s). This highlights the benefits of having open seismic data following international data exchange standards (Suarez et al., 2008), which allows for engagement of the broader seismic monitoring community, which may help resolve details of earthquake sequences that otherwise would not be possible. Coordinated seismic monitoring between regional, national, and international agencies may provide the best estimates of source parameters that are critical for emergency response applications through long‐term seismic hazard research, assessment, and engineering studies.

The U.S. Geological Survey (USGS) event webpage contains up‐to‐date information of USGS event characterization (https://earthquake.usgs.gov/earthquakes/eventpage/us7000kufc/executive, last accessed December 2023), including origin information, moment tensor modeling, and finite‐fault modeling and data. Details on the impact of the Agadir earthquake can be found at https://earthquake.usgs.gov/earthquakes/eventpage/iscgem878424/impact (last accessed September 2023). The facilities of EarthScope Consortium were used for access to waveforms and related metadata for modeling performed outside of the National Earthquake Information Center’s (NEIC) real‐time system (Patton et al., 2016). These services are funded through the Seismological Facility for the Advancement of Geoscience (SAGE) Award of the National Science Foundation (NSF) under Cooperative Support Agreement EAR‐1851048. The proximal seismic stations used in this study were made available from the Western Mediterranean Broadband Seismic Network, WM, operated by the San Fernando Royal Naval Observatory, Universidad Complutense De Madrid (UCM), Helmholtz‐Zentrum Potsdam Deutsches GeoForschungsZentrum (GFZ), Universidade De Évora (UEVORA, Portugal), and Institute Scientifique Of Rabat (ISRABAT, Morocco; ROA et al., 1995). Copernicus Sentinel data 2023 were retrieved from European Space Agency (ESA) SciHub and Alaska Satellite Facility (ASF) Distributive Active Archive Center (DAAC), and processed by ESA. Interferograms were processed with the Interferometric Synthetic Aperture Radar (InSAR) Scientific Computing Environment (InSAR Scientific Computing Environment [ISCE]) version 2.0.

Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

The authors acknowledge that there are no conflicts of interest recorded.

The authors thank the National Earthquake Information Center (NEIC) 24/7 seismic analysts who are critical to the U.S. Geological Survey (USGS) NEIC’s global earthquake monitoring. Reviews by Walter Mooney and two anonymous reviewers improved this article. The authors are additionally grateful to Editor‐in‐Chief Keith Koper and Ryan Gold, Brian Shiro, and Janet Carter with USGS for their helpful comments.