The 2023 Mw 6.8 Morocco Earthquake: A Lower Crust Event Triggered by Mantle Upwelling?

A M6.8 earthquake struck the High Atlas Mountains in Morocco on 8 September 2023, ending a 63‐year seismic silence. We herein attempt to clarify the seismogenic fault and explore the underlying mechanism for this seismic event based on multiple data sets. Utilizing probabilistic Bayesian inversion on interferometric radar data, we determine a seismogenic fault plane centered at a depth of 26 km, striking 251° and dipping 72°, closely aligned with the Tizi n’Test fault system. Given a hypocenter at the Moho depth, the joint inversion of radar and teleseismic data reveals that the rupture concentrates between depths of 12 and 36 km, offsetting the Mohorovičić discontinuity (Moho) at ∼32 km. Considering a strong link between magma activity and failure in lower crust, we propose that the triggering of the earthquake possibly was mantle upwelling that also supports the high topography.


Introduction
The 2023 Morocco High Atlas earthquake, a catastrophic M6.8 event, marked a seismic reawakening in a region historically quiescent.Occurring on 8 September 2023, this earthquake was the most significant in Morocco since the M5.9 Agadir earthquake in 1960, resulting in approximately 3,000 fatalities and affecting nearly 2.8 million people.Its impact was amplified by a combination of infrequent major seismicity in modern records and inadequate building standards, underscoring a critical gap in regional preparedness and infrastructure resilience.
The diverse geological structures in Morocco include the Rif Mountains, shaped by plate convergence, and the Middle-Atlas, High-Atlas, and Anti-Atlas ranges (Figure 1), sculpted by mantle upwelling processes (Lanari et al., 2023).This mantle-driven activity is central in the genesis of the pronounced topography of the region, with the Atlas system exhibiting peak elevations exceeding 4,000 m.The High Atlas, in particular, is underlain by a crust of 32-40 km thickness (Mancilla & Diaz, 2015;Spieker et al., 2014), a dimension insufficient to naturally support such elevations.According to Airy isostasy, approximately 50 km thickness of crust is required (Missenard et al., 2006;Teixell et al., 2005).The High Atlas represents a doubly vergent intracontinental belt, established during the Cenozoic convergence of the African and Eurasian plates.Two main groups of faults, namely, oblique-slip faults on the orogenic axil and thrusts along the orogenic margins, constitute the transpressional structure in the High Atlas belt.Thermochronological and structural data indicate that abrupt vertical uplift and intense exhumation initiated in the middle/late Miocene contemporary with volcanism (Lanari et al., 2020).These features are typically attributed to processes of inversion tectonics and mantle upwelling (Ellero et al., 2020).Seismically, Morocco is characterized by a low convergence rate of 4 mm/year along the Africa-Eurasia plate boundary (Figure 1a), typically manifesting as earthquakes near the Rif Mountains (Koulali et al., 2011;Peláez et al., 2007).However, the western High-Atlas Mountains display a contrasting tectonic behavior, with a much lower convergence rate of 1 mm/a and sparse seismic activity (Koulali et al., 2011;Lanari et al., 2023).This disparity raises fundamental questions about the seismogenic mechanisms operational in the High Atlas, particularly in the context of the 2023 earthquake.
Inversions of geodetically ground displacements are exclusively capable of determining the most plausible fault geometry at depth (Bagnardi & Hooper, 2018;Wei, Chen, Lyu, et al., 2023;Wei, Chen, & Meng, 2023).In this study, the seismogenic fault geometry is estimated through a fully Bayesian inversion of the co-seismic ground displacements, informed by the interferometric synthetic aperture radar (InSAR) approach.We further incorporate teleseismic data to model the detailed rupture kinematics.Our findings offer new insights into the dynamics of lower crust earthquakes in intraplate settings and challenge existing paradigms of regional seismic hazard assessment.

Methods and Results
The catastrophic M6.8 earthquake that struck the High Atlas region in 2023 provided a unique opportunity to delve into the mechanics of seismic activity in intraplate regions, particularly in areas characterized by mantle upwelling and complex fault systems.Our investigation started with modeling co-seismic deformation with Sentinel-1A SAR interferograms on three tracks (See the detail in Text S2 and Table S1 in Supporting Information S1).The co-seismic deformation maps unveiled by our InSAR analysis offered a vivid representation of the earthquake impact.The observed displacement patterns on the hanging wall of the Tizi n'Test (TNT) fault system fault, highlighted the asymmetrical deformation caused by the seismic event (Figures 2b-2d).The northern uplift (∼15 cm) and southern subsidence (∼5 cm) captured in the InSAR deformation maps were crucial in guiding our understanding of the fault behavior.
The Bayesian inversion approach was employed to estimate the most plausible fault geometry given the observational data.Grounded in Bayes' theorem, this approach allowed us to probabilistically infer the fault geometry and slip parameters, laying the foundation for our detailed analysis of the fault characteristics.Utilizing the quasi-static displacement field from satellite radar interferograms, we adopted a Bayesian inference model (See the detail in Text S3 in Supporting Information S1).This model correlates the observed data with our forward model, which operates within the context of the layered elastic media.We adopted the 1-D velocity model IASP91 (Kenntt & Engdahl, 1991), the top layers of which were replaced by the Crust 1.0 model (Laske et al., 2013), and the forward model was computed with the Pyrocko-GF toolkit (Heimann et al., 2019).
The first-order fault geometry is characterized by a uniform rectangular fault plane.The setting of prior distributions restricts the range of exploration in the parameter space of fault geometry.To adequately consider all fault orientations and dips, we set the prior distributions for fault strike and dip as uniform distributions within 0-359°a nd 1-90°, respectively.This allows us to fairly explore the probability of all combinations of fault strike and dip.To determine the location of the fault plane, we solve for the central position of the fault plane, which represents the average position of fault slip.The results from the U. S. Geological Survey indicate a depth of 27 km for the hypocenter, suggesting that the fault plane may be quite deep.We set the prior distribution for the depth of the fault plane center within a uniform distribution of 1-40 km, to encompass as many possible depths as feasible.For the horizontal position of the fault plane center, we used the epicenter coordinates from the U. S. Geological Survey as the reference position (0 km, 0 km), and then set the prior distribution for the eastward and northward coordinates of the fault center within a uniform distribution of 50 to 50 km from the reference position.In summary, we assigned a fairly large range to the prior distribution of fault geometry parameters (Table S2 in Supporting Information S1).This approach avoids excluding possible fault geometries, enabling us to provide probability distributions for all possible fault geometry parameters.The posterior probability distributions sampled by the sequential Monte Carlo (MCMC) method (Dutta et al., 2021) converge to unimodal probability distributions (Figure S5 in Supporting Information S1).The medians of the posterior samples indicate that the fault plane spans 17 km in width along the dip and 20 km in length along the strike (Figure 2a, Figure S6 and Table S2 in Supporting Information S1).The center of the fault plane is positioned at (8.33°W, 30.97°N) with a depth of 26 km.The orientation of the fault plane presents a strike of 251°and a 72°dip to the northwest.
Extending the fault plane to the Earth's surface indicates that the fault strike nearly corresponds to the active TNT fault (Figure 2a).The dipping angle estimation (72°) from our inversion is also close to that of the TNT fault observed in the field investigation (∼70-80°) (Domènech et al., 2015;Ellero et al., 2020;Lanari et al., 2020).Considering that the earthquake was possibly produced by slip on the TNT fault, we tried to set the surface trace of the TNT fault as the top edge of the seismogenic fault and place the bottom edge of the seismogenic fault at a depth of 40 km.We further conducted a simultaneous Bayesian estimation of the fault slip distribution and the non-planar geometry to investigate whether the earthquake formed the TNT fault through either planar or listric structures (See the detail in Text S4 in Supporting Information S1).The Bayesian estimation indicates, with high probability, that the fault features a planar geometry, with a dipping angle of ∼73°(Figure S8 in Supporting Information S1).Most of the fault slips are distributed within the depth range of 16-32 km.The modeling of slip on the planar fault fits well to InSAR data (Figures 2b-2d).Therefore, the seismogenic fault of the earthquake may serve as a segment in the TNT fault system.Our methodology further encompassed the examination of teleseismic waveform recordings from broadband stations (Figure S9 in Supporting Information S1) within the Global Seismographic Network (GSN), obtained via the IRIS data management center.These records, spanning a complete azimuthal coverage, were instrumental in understanding the far-reaching impacts of this seismic event (See the detail in Text S5 in Supporting Information S1).
Based on the fault geometry from the Bayesian estimation, we employed a joint inversion of the InSAR and teleseismic data to estimate detailed rupture kinematics (See the detail in Text S6 in Supporting Information S1).The hypocenters estimated by the U.S. Geological Survey and the GEOFON project of GFZ Potsdam exhibit a large discrepancy in both horizontal position (Figure 2a) and the depth (19 vs. 27 km).According to our fault plane geometry, their epicenter projected onto our fault plane requires a hypocenter depth larger even than the local Moho depth (32 km, Figure 1c).Therefore, we propose a hypocenter located on the fault near the depth of the interface between upper mantle and Moho (Figure 3a).Our inversion analysis revealed that the earthquake rupture initiated at upper mantle, propagated laterally and shallowly, and then arrested at a depth of 12 km.The final slip distribution spans approximately 40 km in length and 28 km in width, with the slip peak at 1.7 m at a depth of 22-24 km.The cumulative seismic moment of 2.18 × 10 19 Nm, equivalent to Mw 6.82, was released within 18 s, underscoring the profound energy release during the earthquake (Figure 3).
We have derived the kinematic parameters of the earthquake utilizing both InSAR and teleseismic observations.However, it's crucial to acknowledge that our kinematic model is subject to limitations and uncertainties inherent in these data sources.InSAR furnishes us with co-seismic deformation data boasting high spatial resolution, enabling us to better delineate fault slip distribution.On the other hand, teleseismic waveform data, renowned for its high temporal resolution, offers valuable insights into the temporal evolution of slip.Given the relatively deep depth of the earthquake, the resolution of InSAR-derived surface deformation data diminishes as we delve deeper, impacting our ability to precisely map slip distribution.Moreover, the absence of nearby seismic stations for this event restricts the spatial resolution of teleseismic stations, consequently weakening our constraints on fault slip and the earthquake hypocenter.These limitations stem from uncertainties in spatial distance and velocity structure.Although our estimated kinematic model aligns with the available data, it may still inadequately capture the intricacies of the entire rupture process.

Discussion
The 2023 Morocco High Atlas earthquake, a seismic event of magnitude 6.8, provides an invaluable exemplar for understanding seismic hazards in intraplate regions, particularly those influenced by mantle upwelling.This event, occurring in a region characterized by slow deformation rates, underscores the importance of re-evaluating seismic hazard assessments, especially in areas far from plate boundaries.Our comprehensive analysis, integrating geodetic, seismic, and seismicity data, has elucidated the underlying mechanisms of this seismic event.
From a local geological perspective, the geometry of high-angle faults appears more plausible compared to lowangle faults.The gentle-angle model deviates significantly from the surface trace of the North Atlas, (NA) fault and intersects abnormally with another NE-SW striking fault in the paper of (Cheloni et al., 2024).There is also inconsistency with the geological structure that NA fault consists of NNW-directed low-angle thrust surfaces possibly rooted in the AM fault at upper crust (Ellero et al., 2020).Conversely, the fault top of the high-angle model aligns more closely with the surface trace of the TNT fault, as observed in our model (Figure 2a).Cheloni et al. (2024) could not simply exclude the high-angle fault model, based on field observations that the TNT fault is a dextral strike-slip fault with a minor reverse component.The timing of the right-lateral oblique-slip kinematic activity in the axial part of the belt is difficult to constrain because syn deformational deposits are not present there (Lanari et al., 2020).In contrast, active morphotectonic evidences that Precambrian rocks thrust over Quaternary thick slope and fluvial deposits are presented along the TNT fault with vertical offset at least 20 m (Delcaillau et al., 2010(Delcaillau et al., , 2011)).Catastrophic rock avalanches suggest that a 4.5 ka seismic event, linked to the proximity of the major TNT fault, triggered these collapses (Hughes et al., 2014).
Compared to previous studies, we employed the different method for fault geometry inversion.Yeck et al. (2023) estimated the seismic moment tensor with teleseismic waveform data, indicating two possible fault orientations: the high-angle north-dipping fault and the gentle-angle south-dipping fault.They indicated that the high-angle model better fits the InSAR observations.At the same time, Cheloni et al. (2024) employed the simulated annealing algorithm and InSAR observations to search for fault geometry in an elastic homogeneous half-space and analyzed both the high-angle and gentle-angle models.Here we provided probabilistic estimates of fault geometry parameters from a globally optimal perspective within the context of the local layered elastic media.We explored the probability distribution of fault geometry parameters in Bayesian sense.Our prior distributions for fault strike, dip angle, and depth parameters were set very wide, allowing us to explore all possible solutions including both the high-angle north-dipping fault and the gentle-angle south-dipping fault in a unified Bayesian inference framework.We then employed a sequential Monte Carlo method to sample the probability density distribution of parameters (Dutta et al., 2021).This method combines simulated annealing and Markov chain techniques, providing a more reliable exploration of the parameter probability distribution space and allowing us to examine trade-offs between different parameter values (Minson et al., 2013).Our results indicate a clear unimodal probability distribution, with parameter values corresponding to a north-dipping high-angle fault model being highly probable (Figures S5 and S6 in Supporting Information S1).
The earthquake initiated at top of upper mantle with the slip peak at 1.7 m at a depth of 22-24 km in the lower crust poses significant questions about the genesis of seismicity in such environments.Typically, top of the upper mantle or lower crustal earthquakes are considered rare due to the prevailing high temperatures, which are generally thought to inhibit brittle failure.Such seismicity is often attributed to fluid or magmatic processes, as evidenced in regions like the East Africa and Taupo Volcanic Zone (Lapins et al., 2020;Reyners et al., 2007).
In the context of the High Atlas region, the mantle upwelling has been a central factor in shaping the regional topography.Studies suggest that this upwelling has contributed to the high relief of the Atlas orogeny (Lanari et al., 2023;Missenard et al., 2006;Teixell et al., 2005).From receiver functions and shear-wave splitting, Miller and Becker (2014) determined that localized, nearly vertical offset deformation of both crust-mantle and lithosphere-asthenosphere interfaces at the flanks of the High Atlas.Mancilla and Diaz (2015) also provided high resolution Moho topography map beneath the western High Atlas Mountains.It shows that Moho depth is only 32 km between the AM fault and the TNT fault, indicating about 4 km throw between hanging wall and footwall along the TNT fault and the AM fault (Figure 1c).Furthermore, the high-angle faulting on the TNT fault, as revealed by our inversion from the earthquake, indicates a deep-reaching offset to Moho and the top of upper mantle.Lee et al. (2022) incorporated anisotropy as an a priori constraint in teleseismic P-wave travel-time tomography as mantle flow can develop seismic anisotropy.Their improved results are more consistent with the hypothesis that the low-velocity anomalies come from the mantle material of the Canary Plume.Especially, it shows clear low-velocity anomalies across the western High Atlas (near the seismogenic zone) from the depth slice of isotropic P-wave tomography image at 75 km.Based on these geophysical data and our findings, we proposed that mantle upwelling, accompanied by the injection of pore fluids or magma along a preexisting vertical fault, constitutes the seismic dynamics of the 2023 M w 6.8 earthquake in Morocco.This geological setting implies that fluids or magma might intrude through these weakened channels.Our finite source inversion analysis sheds light on the significant role of the TNT fault in this seismic event.The high-angle faulting on the TNT fault, as revealed by our analysis, indicates a deep-reaching offset to Moho (Figure 3).Based on the published geophysical data (Mancilla & Diaz, 2015;Miller & Becker, 2014;Spieker et al., 2014), we posit that TNT fault possibly extend to the interface of the lithospheric mantle and asthenosphere.Then pore fluid and magma associated with active hot mantle upwelling could infiltrate and permeate the fragmented wall rocks along the TNT fault, thereby lowering the threshold for fault failure and triggering the earthquake in the lower crust (Figure 4).
The implications of the 2023 Morocco earthquake extend beyond the immediate seismic hazard.They challenge the traditional perceptions of intraplate seismicity and emphasize the need for continuous monitoring and reassessment of seismic risks in similar geological settings.The event occurrence at a significant depth, along with its association with mantle processes, offers a new perspective on the dynamics of seismic activity in the lower crust.It highlights the importance of integrating diverse geological, geophysical, and seismological data to comprehensively understand and anticipate seismic hazards in complex tectonic environments.Our findings, while offering vital insights into the mechanics of the 2023 earthquake, also serve as a reminder of the latent seismic threats in regions with slow deformation rates.This event underscores the necessity for improved infrastructure resilience and preparedness strategies, especially in areas with infrequent but potentially devastating seismic activity.The High Atlas earthquake thus stands as a testament to the dynamic and often unpredictable nature of Earth's tectonic processes, urging a reevaluation of our approaches to understanding and mitigating seismic risks.

Conclusion
The 2023 Morocco High Atlas earthquake, a seismic event of magnitude 6.8, stands as an exemplary case in understanding intraplate seismicity under the influence of mantle upwelling.This rare lower crust earthquake, far from plate boundaries and in a region with a low convergence rate, challenges long-standing paradigms of seismic hazard.The earthquake, causing substantial loss of life and impacting millions, underscores an urgent need to reassess our understanding of seismic risks, especially in regions where such events are historically sparse.Our comprehensive study, employing an integration of geodetic, seismic, and seismicity data, reveals that the seismic event ruptured at a depth of 12-36 km, with a maximum slip of approximately 1.7 m.This event is distinctive not just for its depth and location but also for its association with the mantle upwelling beneath the Atlas region.This upwelling, essential in sculpting the regional topography, also plays a crucial role in localizing crustal shortening along pre-existing structures, thereby contributing to the genesis of this significant seismic event.
The role of mantle upwelling in facilitating seismic activity in the lower crust, as evidenced by our findings, opens new avenues in seismic hazard assessment.Traditional views have often overlooked such deep-seated processes in intraplate regions, focusing more on plate boundary dynamics.However, the deep crustal origin of the 2023 Morocco event, linked to mantle dynamics, suggests a need to broaden our perspective.This includes considering the influence of mantle processes in seismic hazard models, particularly in regions exhibiting similar geological and tectonic settings.Moreover, our study highlights the critical importance of continuous monitoring and advanced analysis techniques in seismic risk assessment.The integration of state-of-the-art geodetic and seismic methodologies provided an unprecedented view of the earthquake mechanics, offering insights essential for future hazard mitigation strategies.Such strategies are particularly vital for regions with infrequent but potentially devastating seismic activity, as demonstrated by the impact of the 2023 High Atlas earthquake.
In conclusion, the 2023 Morocco High Atlas earthquake serves as a stark reminder of the latent seismic threats in regions with slow deformation rates and complex tectonic settings.It calls for a reevaluation of seismic hazard assessments in similar geological contexts, emphasizing the integration of diverse data sets and the consideration of deeper, mantle-related processes.The lessons learned from this event are crucial not only for Morocco but for global seismic risk assessment, urging us to prepare better for the dynamic and often unpredictable nature of Earth's tectonic processes.Mancilla and Diaz (2015); the nearly vertical offsets of both the crust-mantle and lithosphere-asthenosphere interfaces are cited from Miller and Becker (2014).Mantle upwelling and lateral flow are testified by Lee et al. (2022).This conceptual model is inspired by Lanari et al. (2023).Fluid and magma, associated with a hot mantle upwelling, intrude and pour into fragmented wall rocks through TNT fault.This process reduces the threshold for fault failure, and consequently offsetting the Moho and lower crust and generating the Morocco Earthquake.

Figure 1 .
Figure 1.Seismotectonic background of the M w 6.8 earthquake in Morocco on 8 September 2023.(a) Morocco undergoes the convergence of Eurasian plate and African plate.(b) Tectonic geology in Morocco, including Rif, Middle Atlas, High Atlas and Anti-Atlas Mountains.The blue dashed box outlines the tectonic background for the M w 6.8 earthquake (c).Note that the Moho depth is only 32 km between AM fault and TNT fault in the western High Atlas (Mancilla & Diaz, 2015).The detailed fault characteristics around the seismogenic zone in the field are included in the Text S1 in Supporting Information S1.Epicenters estimated by the U.S. Geological Survey (USGS) (blue star) and GEOFOrschungsNetz (GEOFON) operated by GFZ (black star) are indicated.Abbreviations of faults: TNT, Tizi n'Test; NA, North Atlas; AM, Adassil-Medinet; SA, South Atlas.

Figure 2 .
Figure 2. The fault geometry from Bayesian inversion.(a) The projection of the uniform rectangular fault plane on the Earth' surface.The green line represents the segment intersected by extending the fault plane to the surface.Black and blue stars are the epicenter estimated by the U.S. Geological Survey (USGS) and GEOFOrschungsNetz (GEOFON) operated by GFZ, respectively.The inset map shows the 3D view of the estimated fault plane and the distribution of aftershocks 3 months after the mainshock.(b-d) The fit to the data for the simultaneous Bayesian estimation of fault slip distribution and fault geometry, with assuming the TNT fault as the top edge of the fault.The panels in the first column show the radar line-of-slight surface displacements from ascending track path 45, descending track paths 52 and 154.The corresponding model predictions (the second column), and residuals (the last column) are color-coded on the same scale.The solid black line in each panel indicates the top edge of the fault plane.LOS, line of sight.Positive pixel values indicate motion toward the satellite.

Figure 3 .
Figure 3.The rupture kinematics resulting from the joint inversion of interferometric synthetic aperture radar and teleseismic data.(a) Snapshots of the rupture evolution at different time intervals, with the rupture velocity empirically set at 2.8 km/s (0.8 times the S-wave velocity) (Heaton, 1990).(b) Moment release over time.(c) The final fault slip distribution within the tectonic setting.

Figure 4 .
Figure 4. Schematic model of the seismogenic mechanism leading to the Morocco Earthquake.Moho data are derived fromMancilla and Diaz (2015); the nearly vertical offsets of both the crust-mantle and lithosphere-asthenosphere interfaces are cited fromMiller and Becker (2014).Mantle upwelling and lateral flow are testified byLee et al. (2022).This conceptual model is inspired byLanari et al. (2023).Fluid and magma, associated with a hot mantle upwelling, intrude and pour into fragmented wall rocks through TNT fault.This process reduces the threshold for fault failure, and consequently offsetting the Moho and lower crust and generating the Morocco Earthquake.