INTEGRATION OF GRAVIMETRIC AND FIELD GEOLOGICAL INVESTIGATIONS FOR STUDYING TECTONICS AND GEODYNAMICS IN THE CENTRAL PART OF THE ATLAS INTRACONTINENTAL FOLD SYSTEM (AZILAL PROVINCE, MOROCCO)

1. INTRODUCTION

The Central High Atlas of Morocco is one of the most characteristic intracontinental mountain belts in North Africa [Mattauer et al., 1977]. Its current structure is the result of several tectonic phases that overlapped in time, with Paleozoic inherited faults repeatedly influencing the style and location of the later deformations. The region was first shaped as a Mesozoic rift basin linked to the opening of the Central Atlantic, and then inverted strongly during the Cenozoic Alpine orogeny [Frizon de Lamotte et al., 2008]. This inversion gave rise to the formation of the NE-SW fold and thrust belt that is largely rooted in the Paleozoic basement discontinuities [Beauchamp et al., 1996; Fekkak et al., 2018; Laville, Piqué, 1992; Teixell et al., 2003].

The structural evolution of the High Atlas is profoundly influenced by inherited structures, including the Paleozoic basement fabrics and the Triassic evaporite layers. These pre-existing weak lithospheric zones were reactivated during multiple tectonic events, initially during the Mesozoic extension and later during the Cenozoic compression. The present geometry of the chain clearly reflects this inversion of earlier extensional structures [Beauchamp et al., 1996, 1999; Skikra et al., 2025]. This mechanism is not unique to Morocco: similar tectonic scenarios were described in the Pyrenees, Iberia, and Tunisia [Bouaziz et al., 2002; Rosenbaum et al., 2002; Teixell, 1998], showing that structural inheritance is a common feature in intracontinental orogeny.

Within this regional framework, the Azilal Province constitutes a key sector of the Central High Atlas (Fig. 1). Its geological architecture is particularly complex, being controlled by inherited basement faults, Triassic evaporitic levels acting as décollements, and the superposition of several extensional and compressional phases [Missenard et al., 2007]. Previous studies highlighted large-scale folds, thrust contacts, normal faults and diapiric structures such as the well-known Talmest diapir, which exemplifies the role of Triassic salt in accommodating deformation [Ettaki et al., 2007; Ibouh et al., 2001; Martín-Martín et al., 2017]. Additional observations in the Tansrift area have revealed copper mineralization associated with the red beds, further illustrating the interplay between tectonics, diapirism, and fluid circulation [Ibouh et al., 2011]. Complementing these surface observations, geophysical investigations are considered crucial for resolving the subsurface geometry and assessing the structural control of inherited discontinuities and evaporitic décollements.

Therefore, geophysical methods are a very useful complement to geological mapping. The Bouguer anomaly maps help to identify structural uplifts and subsiding basins, while processing techniques such as the Horizontal Gradient Magnitude (HGM), CET (Centre for Exploration Targeting at the University of Western Australia) edge detection and Euler deconvolution provide critical constraints on the geometry of subsurface faults and allow depth estimates of crustal features to be made [Cooper, Cowan, 2006; Holden et al., 2008; Reid et al., 1990; Sassioui et al., 2023]. These methods have already been applied successfully in other Moroccan basins, for example, Tadla and Missour, where they revealed the influence of the Triassic evaporites and basement uplifts on basin development [Nait Bba et al., 2019; Najine et al., 2006].

In this study, we apply an integrated approach in the Azilal Province – a part of the Central High Atlas – by combining gravity data with detailed field observations. The Bouguer and residual anomalies are compared with surface geology. The gravimetric lineament analysis, systematically comparing these datasets, aims to determine the degree to which the inherited Paleozoic basement faults are present there. The Triassic evaporitic décollements controlled the regional polyphase tectonic evolution. Ultimately, this research seeks to provide new constraints on the role of structural inheritance in intracontinental inversion tectonics, with implications for the geodynamic evolution of the Atlas system and similar orogenic belts.

2. GEOLOGICAL AND TECTONIC CONTEXT

The Azilal Province, located in the Central High Atlas of Morocco, is part of an intracontinental orogenic system shaped by the complex interaction between the inherited basement structures and the subsequent Mesozoic-Cenozoic tectonic events [Mattauer et al., 1977; Michard, 1976; Piqué et al., 2002]. The region provides valuable insights into crustal reactivation, tectonic inversion, and sedimentary basin development in an intraplate setting.

The geological basement consists mainly of the Precambrian and Paleozoic metamorphic units, including schists and quartzites which are exposed locally in the southern part of the province. These units, affected by the Hercynian orogeny, later served as structural templates for Mesozoic rifting and Cenozoic inversion phases [Beauchamp et al., 1999; Laville, Piqué, 1992]. The Triassic record marks the onset of rifting, with thick red beds, evaporites (gypsum and halite), and basaltic flows related to the opening of the Central Atlantic [Beauchamp et al., 1996].

The Jurassic sedimentation was dominated by shallow marine carbonates, mainly limestones and dolomites accumulated on extensive platforms during a long period of relative tectonic stability. These platforms are particularly well exposed in the Azilal region, as illustrated in the geological map of the study area (Fig. 2). During the Late Jurassic and Early Cretaceous, the basin’s development was marked by a transition to mixed carbonate-siliciclastic sedimentation, indicating local tectonic disturbances and syn-sedimentary faulting [Frizon de Lamotte et al., 2009].

The Cenozoic evolution of the region is characterized by the deposition of thick fluvio-lacustrine successions in syn-orogenic basins. These deposits document the inversion of former extensional faults into compressive structures during the Alpine orogeny, giving rise to the ENE-WSW trending folds and associated thrusts [Hafid, 2000; Teixell et al., 2003]. Recent studies [Moussaid et al., 2023; Oujane et al., 2025] emphasize the persistence of active deformation in the Azilal sector, evidenced by uplifted structures and ongoing seismic activity. Altogether, the Azilal region offers a comprehensive geological archive for understanding the structural inheritance and tectonic inversion processes in the Central High Atlas. It remains a key area for exploring the dynamics of intracontinental deformations and the evolution of complex fold-and-thrust belts in the High Atlas, Morocco.

3. DATA AND METHODOLOGY

The gravity data used in this study were extracted from the WGM2012 global gravity model [Bonvalot et al., 2012], which integrates satellite, terrestrial, and marine observations to provide a consistent representation of the Earth’s gravity field. We worked with the complete Bouguer anomaly provided by WGM2012 (reduction density 2.67 g/cm³) at 2'×2' (~3.7 km at this latitude). Because the global Bouguer product relies on a ~1–2 km DEM for terrain corrections, it is best suited to regional signals. We therefore restrict our interpretations to the first-order, crustal-scale features rather than to small individual structures. Grids were projected to WGS84 UTM 29N for analysis, resampling for display (e.g., to 500 m) does not add the resolution. With regard to regional-residual separation and filtering to emphasize map-scale structural variations while minimizing short-wavelength terrain leakage, we separated regional and residual components using upward continuation rather than a polynomial detrend. Specifically, we continued the Bouguer field by 5, 10, and 15 km upward and used a 10 km upward-continued grid for edge detection and lineament mapping, which suppresses wavelengths shorter than ~20 km. Results are shown only where features are stable across the 5–15 km continuation range. This approach focuses the analysis on wavelengths that we can resolve with the WGM 2012 grid in mountainous terrain. To characterize the structural architecture of the study area, a combination of gradient-based, edge-detection, and depth-estimation techniques was applied within the Oasis Montaj environment. The workflow involved (1) the HGM to highlight lateral density contrasts, (2) CET grid analysis to extract lineaments automatically, and (3) Euler deconvolution to estimate depth and geometry of the causative sources. Structural field measurements, obtained for folds, thrusts and normal faults, were used to validate and refine the geophysical interpretations. Although the 2'×2' spatial resolution (~3.7 km) of the WGM2012 model reduces the probability of detection for very small or shallow faults, this study focuses on crustal-scale and major regional structures. Furthermore, integrating multiple filtering techniques and field observations increases the reliability and accuracy of the results.

3.1. Horizontal gradient gravity magnitude

The HGM method was applied to the residual field to enhance lateral density contrasts such as major faults, basement steps, and basin margins [Blakely, Simpson, 1986; Salem et al., 2008]. HGM is defined as:

(1)

High HGM values correspond to zones of abrupt density change, typically associated with structural discontinuities, while low-gradient regions represent more homogeneous domains. The technique is robust, stable, and less sensitive to noise, making it particularly effective in tectonically complex regions [Blakely, Simpson, 1986; Salem et al., 2008]. In Morocco, similarly applied HGM method enabled successful delineation of fault systems, basin margins, and major tectonic contacts [Hafid et al., 2006; Khattach et al., 2013], supporting its validity for this study. The HGM map generated here clearly outlines the main structural trends and provides valuable constraints for subsequent analyses.

3.2. Lineament extraction via CET grid analysis

Automatic lineament extraction using the CET grid analysis plugin was developed by the Centre for Exploration Targeting [Holden et al., 2008] in Oasis Montaj. This tool enhances directional textures within gridded datasets through multiscale filtering and variance-based edge detection, allowing the automatic extraction of linear features that likely represent faults, fractures, or shear zones. The algorithm computes local orientation tensors and directional contrasts to identify continuous alignments that may correspond to structural discontinuities.

This approach is particularly suited for mountainous and structurally heterogeneous regions, where traditional visual interpretation may fail to capture subtle linear trends. In this study we used multiscale filtering to enhance directional textures, and then applied the "Find lineaments" tool with a contrast threshold of 0.70 (70 %), minimum segment length of 5 km, and a gap-closure/simplification tolerance of 2.0 km. CET analysis yielded a dense and coherent network of lineaments, which were further analyzed using rose diagrams to determine the dominant orientation families. The resulting lineament map displays several preferred directions that coincide with the main tectonic patterns of the study area. Comparable studies in the Moroccan and North African contexts [Bencharef et al., 2022] have demonstrated the effectiveness of CET analysis for mapping regional fracture systems and subsurface fault zones, confirming its reliability within the tectonic framework considered herein.

3.3. Euler deconvolution for source depth estimation

The 3D Euler deconvolution method [Reid et al., 1990] was applied to estimate the depth and spatial distribution of subsurface density sources responsible for the observed gravity anomalies. This technique is based on Euler’s homogeneity equation, which relates the gradients of the potential field to the location and depth of the causative bodies through the Structural Index (SI):

(2)

where N is the structural index, B – the regional field, and (x₀, y₀, z₀) – the coordinates of the source.

In this study, Euler solutions were calculated with SI=0 (contact-type sources), a window size of 10×10 cells, and a tolerance of 15 %. Solutions with depth uncertainty >20 % or low signal-to-noise ratio were discarded. These parameters provide a good balance between resolution and stability, minimizing spurious results but preserving genuine subsurface structures. The derived Euler solutions delineate clustered zones that correspond to deep-seated fault systems and major lithological boundaries. The depth estimates were cross-validated with the HGM and CET results, ensuring internal consistency among the different processing techniques. The interpretation of clustered Euler solutions helped to define the geometry, depth, and continuity of key tectonic features.

3.4. Field trip for validation and interpretation

A field campaign across the Azilal segment documented the main fold and fault systems used to validate geophysical interpretations. We sampled the northwestern part of the Central High Atlas, to validate and interpret geophysical findings. Field observations focused on identifying and characterizing fault zones, measuring structural orientations, and assessing their consistency with geophysical lineaments. The data were collected using GPS, a geological compass, a geological hammer, and QField for visualization, complemented by detailed photographic documentation. These observations provided ground-truth evidence to support the interpretation of the gravity results.

4. RESULTS

4.1. Bouguer and residual anomalies

The Bouguer anomaly map of the Azilal Province (Fig. 3, a) shows clear lateral variations, with values ranging from about –17.5 mGal to +52.3 mGal. Broadly spreading positive values dominate the central and northern domains, whereas a wide negative domain is typical of the southeast. The long-wavelength patterns therein are most simply explained by variations in basement depth and the thickness proportion of relatively low-density Mesozoic cover. Positive anomalies are thick and uplifted; having localized contributions from denser mafic rocks (Triassic basalts, Jurassic gabbros) and compared the Bouguer anomaly map (Fig. 3, a) with the geological map (see Fig. 2), we observe that positive anomalies are concentrated in the central and northern parts of the area dominated by dense Jurassic carbonate ridges and uplifted basement blocks. Negative anomalies broadly coincide with the areas of a thicker sedimentary cover and a deeper basement, especially in the southeastern sector. In the saltinfluenced domains such as Talmest, gravity lows may also reflect the concentration of lowdensity evaporites along the diapiric structures and adjacent minibasins rather than along simple sedimentfilled depocenters. By contrast, the Ouaouizaght locality corresponds to gravity highs, consistent with the uplifted Jurassic carbonate units and a relatively shallow basement. Dense mafic units (Triassic basalts, Jurassic gabbros) could locally enhance the gravity field where they are laterally extensive, although their contribution is likely muted at the 2' WGM2012 resolution. Overall, the long-wavelength gravity patterns are consistent with first-order structural segmentation of the province.

The residual anomaly map (Fig. 3, b), evaluated on the 10 km upward-continued field to suppress short-wavelength terrain leakage, emphasizes gradients and edges associated with major structural boundaries at wavelengths ≥~20 km. Residual positive zones broadly coincide with folded carbonate highs and basement culminations, whereas residual negative zones align with minibasinal domains and diapiric corridors. Importantly, density reductions in salt-influenced systems may localize along salt walls and diapirs due to evaporite concentration, while adjacent minibasins can be comparatively denser; hence, residual lows near the mapped diapiric structures need not imply sediment-filled depocenters. In the Talmest area, the pronounced residual low coincides with the mapped diapiric belt and is consistent either with (i) low-density evaporite concentrations along the salt structures or (ii) a locally deeper basement beneath the adjacent minibasins. By contrast, the Ouaouizaght locality sits within a residual high, consistent with shallow, uplifted carbonate units and a relatively shallow basement. Overall, the Bouguer and residual patterns partition the province into uplifted carbonate-dominated blocks and regionally intervening basinal/diapiric domains.

4.2. Lineament extraction (CET and HGM)

The extraction of lineaments from gravity data clarifies the main structural orientations in the Azilal Province. The CET grid analysis (Fig. 4, a) highlights discrete and segmented lineaments distributed across the area. The rose-diagram of the extracted lineament directions confirms three principal families: a dominant NE-SW set, with subordinate E-W, NE-SW, and NW-SE trends. In addition to straight traces, CET method also reveals curvilinear features, interpreted as fold hinges or influenced zones, particularly in the southeastern part of the province.

The HGM map (Fig. 4, b) enhances sharper and more continuous edges as compared to the CET results. The HGM-derived lineaments track major lateral density contrasts and, within tolerance mapping, are broadly coincident with the mapped thrust boundaries and extensional faults, as well as with salt-related structures, diapir/walls and adjacent minibasins on the geological map. The rose-diagram emphasizes the predominance of NE-SW lineaments, with a secondary NW-SE set. The major NE-SW-striking thrust contacts align with this dominant trend.

The continuity and density of HGM edges support the view that basement-controlled structures segment the province into tectonic blocks. Because HGM delineates edges of high-density bodies, the persistence of these edges over long distances and their alignment with inherited basement grain (Variscan) are consistent with deep-seated features rooted in the Paleozoic basement that partitioned the sedimentary cover into discrete blocks.

Taken together, CET efficiently identifies localized, curvilinear features, whereas HGM emphasizes the continuity of major linear structures. Their agreement yields a coherent tectonic framework in which gravityderived lineaments broadly correlate with mapped geological boundaries and the distribution of density-contrasting facies. Given the 2' WGM2012 resolution, we interpret these patterns at regional scale and avoid assigning shortwavelength features to specific small faults unless corroborated by field evidence.

4.3. Euler deconvolution

Euler solutions (Fig. 5) cluster into three depth bands at the scale of the WGM2012 grid: <1 km, 1–2 km, and 2–4 km. The shallow solutions (<1 km) are only considered indicative, given the 2' resolution and potential short-wavelength terrain leakage; they are commonly plotted along the mapped fold hinges and diapiric belts but they are not quantitative.

Intermediate-depth solutions (1–2 km) form continuous belts that align with the principal lineaments extracted from the HGM and CET analyses (typically within a few hundred meters). We interpret them as upper-crustal contacts and fault zones that bound carbonate highs and adjacent minibasins. Absolute depths should be regarded as order-of-magnitude estimates consistent with the dataset resolution.

The deepest solutions (2–4 km) are concentrated in the southern sector and coincide with the major NE-SW gravity edges. We interpret these clusters as basement-rooted faults or steps in the Paleozoic basement, reactivated during the Atlas inversion. Their distribution supports a basement-controlled segmentation that influenced basin geometry and localization of salt structures.

To enhance robustness, we retained only clustered solutions persisting across reasonable parameter choices and discarded isolated picks or those with large depth uncertainties. Within these bounds, depth variations are small relative to the separation between the three bands, so we use Euler primarily to map relative depth domains rather than precise source depths.

4.4. Field Observations

Field observations (Fig. 6, 7) were used to ensure that surface structures correspond to gravity-derived lineaments. Major faults in the area, including the Azilal-Azourki fault, the Talmest fault zone, and the North Atlas fault, have orientations ranging from NE-SW to E-W and align, within mapping tolerances, with the principal lineaments highlighted by the HGM and CET analyses. Their geometry and along-strike continuity are consistent with basementinvolved structures inherited from the Jurassic rifting and reactivated during the Atlas inversion; we note this as a working interpretation pending additional subsurface constraints.

The Jurassic formations display NE-SW-trending folds and fracture sets. Contractional structures predominate in the studied outcrops, including reverse faults and faultpropagation folds (see Fig. 6, b; Fig. 7). Minor extensional faults do occur locally within the Upper Lias-Dogger series, but they are of secondary importance and likely represent relict rift-phase features. Despite a small scale, the orientations and spatial distribution of these structures closely match the NE-SW lineaments mapped from the HGM/CET results, providing surface evidence consistent with the geophysical interpretation.

The uplifted, high-elevation Jurassic carbonate blocks coincide with the Bouguer gravity highs, consistent with structural culminations above relatively shallow basement. Taken together, the field and gravity datasets show that structural inheritance left a strong imprint on later deformation and basin segmentation which likely influenced fluid-migration pathways during the subsequent Atlas inversion.

Fig. 8 summarizes the orientations of the main fault systems, comparing rose diagrams from gravity-derived lineaments (Fig. 8, a) and geological field measurements (Fig. 8, b). Both datasets show a dominant NE-SW trend, consistent with the regional structural grain inherited from the Paleozoic basement faults and reactivated during the Mesozoic extension and Alpine compression. Subordinate NW-SE and E-W to ENE-WSW families are also present. The close agreement between the gravity and field datasets (orientation differences typically ≤15° at the WGM2012-resovable regional scale) supports the robustness of the integrated interpretation and indicates that inherited structures controlled both surface deformation and deeper crustal architecture.

5. DISCUSSION

Field observations, combined with gravity results, provide new constraints on how inherited structures shaped the present architecture of the Azilal Province. The orientation and geometries of the mapped faults, together with the deformation affecting the Jurassic strata, indicate that pre-existing rift-related discontinuities exerted a persistent control at a later tectonic phase. This influence is especially evident around the Talmest diapir, where a pronounced residual gravity low and clusters of intermediate-depth Euler solutions coincide with the diapir surface trace. This spatial correlation suggests that salt mobilization and diapir rise were facilitated by differential subsidence along the inherited fault zone, with the feeder zone rooted within the Triassic-Liassic interval. Likewise, the Bouguer gravity highs over the uplifted Jurassic limestone plateaus are consistent with the competent blocks above relatively shallow basement steps. Overall, these relationships highlight that tectonic inheritance and lithological heterogeneities exerted a first-order control on the subsequent deformation, basin segmentation, and localization of diapirism in the study area [Ibouh et al., 2001; Martín-Martín et al., 2017; Missenard et al., 2007].

The Bouguer and residual anomaly maps show marked gravity contrasts that broadly follow major lithological units and structural domains. Dense Jurassic carbonates tend to produce positive anomalies, whereas low-density Triassic evaporites and Jurassic marls contribute to negative anomalies. However, these contrasts do more than mirror density differences: they also reflect rheological layering and its role during both rifting and subsequent inversion. Competent carbonate successions behaved as relatively rigid blocks, accommodating deformation by block uplift and flexure and yielding prominent gravity highs. In contrast, the mechanically weak Triassic evaporites localized strain, acted as effective décollements, and enabled vertical flow. Residual lows near the Talmest diapir coincide with the diapir corridor and adjacent minibasins and can reflect both the concentration of low-density evaporites along the salt structures and local variations in the basement depth [Ayarza et al., 2005]. These relationships suggest that lithology-dependent rheology strongly influenced fault propagation, fold geometry, and the location of diapirs during the tectonic evolution of the Azilal region.

The distribution of these anomalies on the updated geological map (Fig. 9) underscores the joint control of lithological heterogeneity and tectonic segmentation. Comparable gravity patterns, documented in the region [Boutirame et al., 2018; Saura et al., 2025; Najine et al., 2006; Vergés et al., 2017] and other sectors of the Central High Atlas [Assoussi et al., 2025a, 2025b; Ayarza et al., 2005; Bouzekraoui et al., 2024], confirm that the anomaly signatures identified in Azilal Province reflect regional-scale tectonic patterns rather than local artifacts. Comparable patterns elsewhere in the High Atlas show low residual anomalies tracking diapiric corridors and adjacent minibasins, as well as highs coinciding with carbonate ridges and basement culminations.

The analysis of lineaments derived from the CET grid and HGM filters (see Fig. 4) highlights three dominant structural orientations, which we classify according to their inferred age, depth range, and tectonic role.

The NE-SW trending faults (Variscan inheritance), corresponding to reactivation of basement discontinuities detected primarily in the deep Euler tier, controlled both Mesozoic rift architecture and Cenozoic thrust localization. This interpretation is supported by previous studies in the region [Frizon de Lamotte et al., 2000; Laville, Piqué, 1992; Michard et al., 2010; Perez et al., 2019; Si Mhamdi et al., 2021]. The NW-SE trending faults (Jurassic transfer zones) corresponded to intermediate-depth structures that accommodate oblique extension during rifting and later acted as lateral ramps during inversion. They are interpreted as paleo-transfer zones that accommodated oblique deformation during the Jurassic rifting [Beauchamp et al., 1996; El Kochri, Chorowicz, 1996; Hafid, 2000; Lanari et al., 2020].

E-W trending structures (Alpine compression) are predominantly shallow and represent Cenozoic thrust fronts and fold axes superimposed on the inherited fabric. This dual structural network is consistent with previous studies in the region [Nait Bba et al., 2019] and in the Phosphates Plateau [Ibouh et al., 2001], as well as with those of inversion settings in Tunisia [Bouaziz et al., 2002]. The coexistence of these orientations explains the segmentation of depocenters, the lateral offset of fold axes, and the focusing of diapirism, particularly at Talmest.

Euler deconvolution results (see Fig. 5) provide additional insights into the crustal significance of the identified fault systems. The Results section highlights three depth groups whose geological implications can be further refined. The shallow tier (<1 km) is consistent with near-surface structures such as fold hinges and diapir belts and treated only as indicative, supporting the role of salt tectonics in basement segmentation [Saura et al., 2025; Hudec, Jackson, 2007; Jackson, Talbot, 1986]. The intermediate tier (1–2 km) aligns with HGM/CET lineaments and likely represents inherited faults within the upper sedimentary pile that were active during the Early-Middle Jurassic rifting and reactivated later, during compressions [Hafid, 2000; Haji et al., 2024; Frizon de Lamotte, Leturmy, 2014]. The deepest tier (2–4 km), concentrated in the southern sector, plausibly reflects basement-rooted faults or steps that acted as long-lived zones of weakness. Given the grid resolution and method assumptions, we use Euler primarily to map relative depth domains rather than precise source depths. Absolute depth values should be treated with caution given uncertainties (±10–15 %) as depicted in related literature. And the consistent clustering of solutions into three depth tiers provides a robust first-order framework for distinguishing shallow salt-related structures from deeper basement-involved faults.

Field evidence further confirms that the Upper Liassic marls display signs of reactivation during the Alpine compression, as described in [Ettaki et al., 2007; Frizon de Lamotte et al., 2011], The Talmest diapir provides a particularly illustrative example: it lies at the intersection of NE-SW and NW-SE fault systems, coinciding with a residual gravity low and intermediate Euler depth solutions. This spatial correlation strongly suggests that diapirism was guided by the reactivation of pre-existing discontinuities, with Triassic evaporites acting as efficient décollement levels. Such observations support the close relationship between inherited fault geometries, evaporite mobility, and subsequent deformation.

The updated geological map (Fig. 9) synthesizes gravity-derived lineaments with some unmapped faults segments field-verified based on HGM and CET maps, the three-category depth classification of structures based on Euler deconvolution, and explicit delineation of the Talmest diapir feeder zone based on residual anomaly patterns. These additions provide a more complete picture of fault network controlling basin segmentation and fold localization in the Azilal Province.

Altogether, our updated structural map (Fig. 9) provides new constraints that are consistent with recent studies documenting multiphase strain and ongoing tectonic activity in the Azilal sector [Casas-Sainz et al., 2023; Saura et al., 2025; Moussaid et al., 2023; Román-Berdiel et al., 2023]. Building on these complementary datasets, we propose a polyphase tectonic evolution involving: (1) Jurassic rifting and subsidence influenced by basement-rooted faults and Triassic evaporites, (2) localized diapirism and basin segmentation during continued extension, and (3) Cenozoic Alpine inversion, which reactivated the earlier extensional faults, thrusts and folds. Similar tectonic scenarios have been described in other sectors of the High Atlas belt [Escosa et al., 2021; Si Mhamdi et al., 2024; Skikra et al., 2025].

From a methodological perspective, this study demonstrates the value of integrating multi-scale gravity transformation with field geological observations for characterizing complex inversion systems. The workflow presented here combines the observed-field analysis, residual-anomaly mapping, gradient-based lineament detection and Euler depth estimation, thus providing a transferable and applicable approach to other intracontinental fold belts where the subsurface data remain limited.

This integrated discussion demonstrates that the Azilal Province is a representative case of inversion tectonics, where the inherited basement structures and evaporitic décollements together controlled the style of deformation and the geodynamic evolution. Beyond the regional framework, these findings contribute to a broader understanding of intracontinental mountain belts worldwide.

6. CONCLUSION

This study uses the geophysical and field data to examine the role of inherited structure in controlling tectonic architecture in the Azilal region in Morocco from rifting to inversion. The Bouguer and residual anomaly maps reveal significant density contrasts, with positive anomalies linked to uplifted carbonate ridges and negative anomalies corresponding to subsiding depocenters filled with the Triassic evaporites and Jurassic marls. Lineament extraction using CET grid analysis and HGM filters identifies three dominant fault systems: NE-SW, NW-SE and E-W, which reflect inherited basement structures. The depth solution indicates the upper-to-mid-crustal extension of the major fault within the region. Future quantitative analysis is planned to further validate these interpreted geometries and refine the depth estimates. We conclude that the Azilal Province provides a representative case for investigating tectonic inversion of inherited basement structures and evaporitic décollements which jointly control the deformation style and geodynamic evolution. Beyond the Moroccan context, these results contribute to a broader understanding of intracontinental orogens and carry implications for natural resource exploration, since the structures that governed tectonic evolution also guided fluid migration and mineralization.

7. ACKNOWLEDGEMENTS

The authors are grateful to the International Gravimetric Bureau (BGI) for providing access to the data. Special thanks are due to A.V. Timofeev and the anonymous reviewer for their thoughtful criticism and constructive suggestions, which significantly strengthened the manuscript.

8. CONTRIBUTION OF THE AUTHORS

Souad Assoussi was responsible for conceptualization, methodology, software, data curation, validation, visualization, and writing the original draft. Slimane Sassioui contributed resources, validation, and review and editing of the manuscript. Brahim Oujane and Fadoua Saadaoui handled resources, data curation, methodology, validation, visualization, and review and editing. Lahcen Ousaid, Zineb Aafir, Hanane Ait Haddou, and Amnay Mezzan contributed to validation and visualization. Youssef Hahou provided supervision and validation.

9. DISCLOSURE

The authors declare that they have no conflicts of interest relevant to this manuscript.