SUMMARY
Constraining the knowledge of deep aquifer structure in the Angad Basin (Northeastern Morocco) remains one of the major challenges for successful borehole drilling projects. This study demonstrates how integrating gravity and electrical data can enhance subsurface geological imaging and provide valuable insights into geological structures relevant to groundwater exploration. Various gravity enhancement techniques, including vertical gradient, horizontal gradient, upward continuation, logistic filter, tilt angle and Euler Deconvolution, have been applied to gravity data in the Angad Basin to identify gravity anomalies and tectonic discontinuities. The resulting gravity maps reveal two elongated depressions trending in the NE-SW direction: the Beni Drar depression in the north and the Oujda depression in the south. Both depressions are delineated by strong gravity gradients, indicating the presence of faults. Multiscale analysis of gravity lineaments highlighted sharper features and identified two principal orientations: N070°–085° and N040°–050°, with the latter being more predominant. Faults striking in the N120°–140° are less prominent. These fault systems play a crucial role in dividing the bedrock into distinct horst and graben structures, aligning with the regional tectonic phases observed in northeastern Morocco. The faults are generally sub-vertical, with estimated depth values ranging from 1000 to 2680 m for 42 per cent of the lineaments. Additionally, the reinterpretation of vertical electrical soundings using 2-D inversion methodology provides insight into subsurface resistivity variations. The 2-D resistivity sections generated from this process illustrate vertical and lateral variations in electrical resistivity over distances of up to 30 km and depths of up to 2 km, revealing the geological layers that form the aquifer. These resistivity sections confirm the presence of multiple faults that significantly influence the structural configuration of the study area, validating the geometry of the two depressions previously identified through gravity data. The 2-D gravity modelling further corroborates and reinforces the findings from the resistivity section inversion, demonstrating that the shape of the residual gravity anomaly curve closely aligns with the morphology of the Jurassic roof structure. The results of this study highlight the complementarity and effectiveness of these two geophysical methods in advancing our understanding of the deep geological structure of the Angad Basin, particularly within the Jurassic limestone, which serves as the region's primary deep aquifer. These findings provide valuable insights for future hydrogeological research in the area.
1 INTRODUCTION
Geophysical investigations play a critical role in studying subsurface geological structures to identify potential resources such as groundwater, minerals and petroleum (Wiederhold et al. 2021; Combes et al. 2022; Abd El-Dayem et al. 2023). The economic development of a country often depends on the effective exploitation of its natural resources. To facilitate resource identification, many countries have conducted extensive geophysical studies, leading to the establishment of comprehensive geophysical data bases (Araffa 2013; Sulaiman et al. 2018; Carrier et al. 2019; Al Deep et al. 2021; Wiederhold et al. 2021; Araffa et al. 2023; Maurya et al. 2024). Advancements in interpretation techniques now allow for the reanalysis of older geophysical data, enabling the detection and characterization of previously unidentified subsurface structures (El Gout et al. 2024). By applying modern computational methods and analytical tools to existing geophysical data base, we can significantly enhance our understanding of subsurface geological formations, potentially uncovering new resource deposits. This approach maximizes the value of past geophysical investigations while reducing the need for costly and time-consuming field surveys. As a result, it supports informed decision-making in resource management and promote sustainable economic growth (El-Awady et al. 2016; Nazar & Ghaeb 2019; Sun et al. 2020; Al-Bahrani et al. 2022).
In recent decades, significant advancements in the interpretation of gravimetric anomalies led to more precise subsurface models. These improvements, driven by advanced computational tools and modelling techniques, allow geophysicists to better analyse density variations in rock formations, leading to enhanced resource exploration and a deeper understanding of the tectonic processes that shaped the Earth's crust. To facilitate the interpretation of geological features, modern methods focus on enhancing and analysing gravimetric anomalies using edge detection, gradient analysis and filtering techniques (Saibi et al. 2019). These techniques help geophysicists delineate geological boundaries and features structures with great accuracy. Edge detection algorithms are particularly useful in highlighting abrupt density contracts within the Earth's subsurface. These contrasts often correspond to boundaries between different geological units and can be identified through changes in intensity in gravity anomaly maps (Ekinci & Yigitbas 2015; Zanella et al. 2020; Ekka et al. 2022; Pham et al. 2022; Ghomsi et al. 2022a, b; Ganguli & Pal 2023; Kamto et al. 2023; Narayan et al. 2023). These boundaries can take various geometric forms depending on the geological context: (i) Linear boundaries, typically associated faults, fractures or structural displacements in rock formations. Such lineaments are crucial in identifying fault networks and structures related to mineral exploration and groundwater reservoirs, and (ii) curved or circular features: They may indicate igneous intrusions (plutons), structural domes or impact craters. Circular anomalies can also be linked to hydrothermal systems, influencing mineralization processes. By integrating these enhanced gravimetric techniques, geophysicists can more accurately map subsurface structures, optimize resource exploration and refine models of geological evolution.
This process is based on the use of horizontal and vertical gradients of gravity anomalies. Blakely and Simpson proposed in 1986 a method of localizing borders by determining the maximum of the horizontal derivative of magnetic or gravity grids, which correspond to the locations of the strongest edges or borders in the image. To estimate the depth of these features, Reid et al. (1990) proposed the Euler Deconvolution method. In 1994, Miller and Singh developed the tilt angle filter, a computational method defined as the arctangent value of the vertical derivative of the potential field to its horizontal derivative, to locate the source-body contact which is indicated by the zero contours (Nasuti & Nasuti 2018). Archibald et al. (1999) apply a tool using a wavelet for edge detection and analysis of potential field data. Khattach et al. (2004) utilized a combination of the horizontal gradient and upward continuation methods to perform a multiscale analysis of gravity lineaments. Ferreira et al. (2013) introduced an enhanced tilt angle method that is based on derivatives of the horizontal derivative. These edges analysis not only provide localization of the geological contacts, but also information about their importance and dip direction. In 2021, Pham et al. presented an innovative approach incorporating the logistic function and horizontal gradient amplitude. This method, acting as a boundary detection filter, exhibits notable enhancements in performance and in its capability to yield high-resolution outcomes.
Electrical resistivity surveying methods are widely used in archaeology, hydrogeology, mineral exploration, engineering, structural and environmental studies (Al-Garni 2009; Poongothai & Sridhar 2017; Loke et al. 2018; Fathi et al. 2019; Lowe et al. 2019; Yusuf & Abiye 2019; Loke et al. 2020; Oyeyemi et al. 2020; Maurya et al. 2024). This technique is particularly favoured in groundwater exploration for assessing groundwater potential and determining optimal borehole sites in fractured confined or unconfined aquifers (Muchingami et al. 2012). The method is intrinsically linked to the chemical and physical characteristics of geological formations, including lithology, porosity, permeability and saturation levels (Aizawa et al. 2011; Bernabé et al. 2011; Johnson et al. 2021; Sawayama et al. 2023). Consequently, rock resistivity data serves as an invaluable resource for developing models of subsurface and stratigraphic structures based on their electrical properties.
Many electrical arrays have been proposed in the literature (e.g. Wenner, dipole–dipole, Schlumberger), each with its strengths and weaknesses. These configurations of electrodes vary in terms of reliability, survey costs, acquisition time and the number of current dipoles used (Loke et al. 2013; Martorana et al. 2017; Oyeyemi et al. 2022). In the Schlumberger array, the potential electrode distance ‘a’ (known as M and N) remains fixed, while the current electrodes (known as A and B) are adjusted to vary the distance ‘s’. The spacing ‘a’ is adjusted as needed to address decreasing measurement sensitivity (Martins et al. 2016). The 1-D vertical electrical soundings (VES) method is particularly favoured for aquifer delineation due to its ability to penetrate deeper into the subsurface. Its popularity is attributed to being cost-effective compared to direct borehole drilling, making it exceptionally useful in situations requiring extensive coverage or probing substantial depths (Riss et al. 2011; Laborgne et al. 2021). The Schlumberger array configuration for VES is widely utilized in electrical resistivity surveys due to its ability to probe deeper and provide better resolution compared to the Wenner array. Additionally, it requires less time to deploy (Alao et al. 2022a, 2022b). However, a significant limitation of this method is its inability to account for lateral variations in subsurface resistivity, which can lead to inaccuracies in the interpreted resistivity and layer thicknesses. To achieve more accurate results, a transformation of a series of single VES data into a 2-D or 3-D survey and interpretation model is necessary for enhancing VES interpretation (Uchida 1991; Khalil & Santo 2011; Riss et al. 2011).
Northeastern Morocco is characterized by a Mediterranean semi-arid climate with irregular annual rainfall and extended dry periods. Groundwater use is increasing to meet domestic, agricultural and industrial demands. The Angad Plain contains both shallow and deep aquifers, which are vital for the socio-economic development of the region. The shallow aquifer, composed of Neogene and Quaternary deposits, is generally subjected to high exploitation rates and is highly vulnerable to pollution and drought. According to the Directorate of Water Research and Planning (2014), the water quality of the Angad aquifer has deteriorated in 62 per cent of the wells surveyed. The main factor contributing to this degradation is nitrate contamination (50 per cent). These affected points are primarily located to the northeast of the city of Oujda, downstream from waste water discharge sites. Overall, the organic and bacteriological quality was satisfactory, except for two wells, one of which is located near a landfill. The mineralogical quality was generally average, with only well (N°IRE 42/12) exceeding acceptable levels of conductivity and chloride for drinking water. Consequently, hydrogeological exploration has increasingly focused on the deep aquifer (Jurassic lime stones and dolomites). However, some deep drilling attempts have been unsuccessful, highlighting the need for comprehensive studies to better understand the deep geological structure of the study area.
Classic structural studies have typically relied on direct observation of outcrops or indirect interpretation of digital tectonic studies based on remote sensing analysis. However, these methods are less effective in the plains, where structures are entirely or partially concealed by Cenozoic and Quaternary deposits. To address this challenge, this study uses, gravity and Schlumberger electrical soundings data of the Angad plain (Northeastern Morocco) as a case study. It explores the effectiveness of combining these two types of geophysical data and enhancing their value through modern techniques to gain new insights into the area's geological structure. The findings will be used to investigate the deep carbonate aquifer, which is experiencing a significant decline in groundwater levels due to intensive exploitation and drought (Hydraulic Basin Agency of the Moulouya, Report 2003). In the study area, VES with electrode separations of up to 10 km allow investigations to a maximum depth of 2000 m. By providing a comprehensive depiction of the geological structure of the Jurassic carbonates within the Angad Plain, this study aims to support sustainable groundwater resource management, which is critical for the region's strategic reserves.
2 GEOLOGICAL AND HYDROGEOLOGICAL SETTING
The Angad Plain, located in the northeastern region of Morocco (Fig. 1a), covers an area of approximately 460 km². It is confined between two orogenic belts: the Beni-Znassen to the north and Jbel Hamra to the south. To the west, it is separated from the Bou-Houria Plain by the Jbels Megrez and Harraza (belts). To the east, it extends into Algeria as the Maghnia Plain (Fig. 1b).
Location (a) and Geological map of the study area (in geographic coordinates) (b), north–south cross-section of boreholes showing a rapid dip of Jurassic rocks beneath the plain (c).
The deformed Paleozoic rocks of the study area crop out in isolated inliers, surrounded by Mesozoic and Cenozoic rocks, in the southern flank of Beni-Znassen and the Jbel Hamra (Oujda). The early Paleozoic is represented by pelites, quartzitic sandstone rocks of Ordovician and the Silurian succession characterized by phtanites and black shales (Hoepffner 1987). The Devonian is characterized by a detrital sequence of argillites and greywackes. The ensemble is deformed by the Eovariscan event; this polyphasic folding that occurred at low-grade metamorphic conditions was mainly oriented in N-S, NNE-SSW and NW-SE with vergence to the west. The area is affected by the second post-Westphalian phase folded NE-SW to E-W (Mrini et al. 1992; Hoepffner et al. 2006).
The Mesozoic succession commenced with the deposition of Late Permian-Triassic red beds (mud-flate, claystones and conglomerates) with limestone intercalation and evaporitic sedimentation, unconformably resting on the Paleozoic basement (Oujidi & Elmi 2000). During this period, an extensional stress field led to the reactivation of older faults, manifesting as normal and strike-slip faults. This synsedimentary tectonic activity is recorded as planes of synsedimentary faults, variation of thickness and angular unconformities. The pre-rift deposits are capped by tholeiitic basalt flows. Subsequently, the geological record transitions to the transgressive lower and middle Jurassic periods, which are primarily composed of marine carbonate platforms characterized by overlying lime stones and dolomites (Sublithographic limestone from the Lotharingian-Carixian ages; basal organic detrital grey limestone attributed to the upper Lias; and brown to beige sandstone with intercalations of limestone in the Dogger epoch). This platform gradually evolves, from the Upper Jurassic, into a varied deltaic sequence, comprising marls, sandstones and calcareous shales. These sedimentary formations provide evidence of shallow-water carbonate environments (Lucas 1942; Du Dresney 1979; Rakus & Valin 1979).
The Neogene, occurring above a regional unconformity (Jadid 2000), is characterized by a lower marine sequence, consisting at the base of conglomerates and lime stones, followed by yellow marls. These facies transition towards the basin margins inton reefal or volcano-detritic facies interbedded with basalt flows. The fossil content of the marls dates this assemblage to the upper Tortonian-Messinian. The Plio-Quaternary marked with continental sediments (molassic horizon of conglomerates, clay and fluvial deposits) and the alkaline volcanism that outcrops at the West and south-west of Oujda yields 40k-40Ar isotopic ages range from 5.6 to 1.5 Ma; this activity suggests that the compressive tectonics was only active within remanent fracture zones (Andries et Bellon 1989). The dominant ENE-WSW fault network inherited from the Hercynian orogeny was reactivated during various times and it is responsible for the architecture of the study area to distinct horst and graben structures (Torbi & Gélard 1994). The final lithostratigraphic series is characterized by the youngest deposits: Quaternary continental sediments. This series represents a continental molassic horizon consisting of conglomerates, clays and fluvial deposits.
The Angad plain, structurally considered as a graben between the Beni-Znassen belts to the north and Jbel Hamra to the south (horsts), contains two hydrogeological resources (Carlier 1971):
The shallow unconfined aquifer flows from south to north through post-Miocene terrains, over an impermeable substratum formed by Miocene marls. These marls outcrop at a few points to the northeast of Jbel Hamra. Their depth position in the plain has been identified through several boreholes, geophysical surveys and extrapolations from flow tests on the groundwater table. The structure of the Miocene marl roof is characterized by a flexure in the centre of the plain, oriented NE-SW, which significantly influences the characteristics of the groundwater table and constitutes a hydrogeological boundary that divides the aquifer into two domains: Angad South and Angad North.
The deep aquifer within the lime stones and dolomites of the Lower Lias (Lotharingian-Carixian) and Middle Lias (Domerian). The Aalenian-Bajocian dolomites are hydraulically connected with the Lotharingian and Carixian formations and, in some areas, are part of the aquifer. Their thickness varies due to synsedimentary tectonic processes and major faults likely play a significant role in controlling groundwater flow pathways. The cross-sectional view of boreholes (Fig. 1c) illustrates a rapid dip of Jurassic rocks beneath the plain. While these rocks are visible at the surface in Jbel Hamra, they extend to depths exceeding 700 m in the plain.
According to a seismic refraction survey, four reconnaissance boreholes were drilled. Only one reached the dolomitic lime stones of the Lias. The others did not penetrate beyond the Middle or Upper Jurassic at depths of 200, 320 and 430 m. Another borehole (1968) drilled for stratigraphic purposes did not reach the base of the Miocene marls at a depth of 700 m near the Oujda-Angad airport (Quang & Simonot 1971; Aqil et al. 2010; Zarhloule et al. 2010).
The deep aquifer is generally confined. Its base is formed by Triassic marls and basalts, while its roof consists of clays and marls from the Upper Jurassic or Miocene. It extends to the eastern boundary of the study area (Algerian territory) and is recharged by direct infiltration of effective precipitation on the horsts (Castany 1998; Lahrach 1999; Zarhloule et al. 2001; Zarhloule et al. 2005). The aquifer exhibits favourable hydrodynamic parameters, with a notably high transmissivity of approximately 0.1 m² s−1 and a storage coefficient estimated at around 10-³ (Zarhloule et al. 2010). The aquifer flows with a slope of approximately 0.5 per cent from south to north. In the Angad Basin, significant aquifer potential exists, but the hydrogeologic properties remain largely unknown due to geological complexities and the lack of deep well investigations. The geophysical study enabled us to delineate the aquifer's geometry and identify optimal sites for drilling.
3 GEOPHYSICAL DATA AND PROCESSING
3.1. Electrical resistivity data
VES is a resistivity surveying method where the apparent resistivity is measured at various electrode spacings around a central point. As the size of the electrode array increases, deeper subsurface layers are formed. Plotting apparent resistivity against the electrode spacing provides insights into variations in subsurface resistivity (Binley 2015; Loke et al. 2020).
The study utilizes data gathered from a survey conducted by the ‘Compagnie Africaine de Géophysique’ in 1977, employing the VES method to delineate the vertical distribution of aquifer units and the subsurface lithologic stratigraphy across the Angad Basin. The data set comprises 256 VES deployed across 17 south–north profiles, employing the Schlumberger array with AB extending up to 10 km (Fig. 2a), with VES points spaced approximately 1 km apart.
(a) Electrical data: location of VES and interpreted profiles (F, K, N and TR). (b) Bouguer anomaly (density correction 2.67 g cm−3) as an overlay over digital elevation model (red dots: gravity stations). Both maps are presented in geographic coordinates.
Given the complexity of the geological terrain in the investigated area, determining both vertical and lateral variations in resistivity is crucial for imaging the subsurface. To address this challenge, we employ a method inspired by Riss et al. (2011) to reinterpret the VES data, generating 2-D resistivity models. Utilizing the Schlumberger array data, we construct a resistivity pseudo-section formatted to align with the Wenner-Schlumberger array for 2-D electrical tomography. The inversion process, aimed at determining resistivities, is executed through RES2DINV (Loke 2000). This process involves initially constructing a resistivity pseudo-section by interpolating apparent resistivities at pseudo-depths equivalent to 0.191 × AB (Edwards 1977).
Triangulation with linear interpolation yields satisfactory results. Subsequently, the pseudo-section is sampled following a Wenner-Schlumberger array sequence at depths specified by Edwards (1977). The resolution is contingent upon the separation distance ‘a’ between virtual electrodes; models are formulated utilizing ‘a’ value of 100 m for comprehensive analysis and 50 m or 20 m to facilitate imaging of the shallowest basin sections. The L1 norm (robust) inversion approach was employed to generate 2-D resistivity models. The blocky or L 1 norm optimization method aims to minimize the sum of the absolute values of spatial variations in model resistivity. This approach favours the generation of piecewise constant models, characterized by distinct regions separated by sharp boundaries. The study will present results from four profiles (F, K, N and TR).
3.2. Gravity data
The data set comprises Bouguer anomaly values acquired in 1967 through surveys conducted by the ‘Fondazione Ing. C.M. Lerici company’ on behalf of the Ministry of Energy and Mines in Morocco. The Bouguer anomaly values were measured at 835 points, with a density of reduction of 2.67 g cm−3 (Fig. 2b). The average measurement density is 0.7 stations km−2. The distribution of gravity data is not uniform due to the challenging topography of the mountainous areas (Beni Znassen and Jbel Hamra), where the density of gravity measurements is significantly lower than in the plains.
To make the data more manageable, the Bouguer anomaly values were interpolated to a regular square grid of 250×250 m cells, for which the Kriging interpolation method was used. Several techniques and methods are used to enhance and analyse gravity anomalies to facilitate the interpretation of geological features and structures.
3.2.1. Separation of regional and residual field
The Bouguer anomalies encompass both residual or local components of gravity, as well as its regional or long wavelength parts, referred to as deep sources. The regional-residual separation of the potential field holds significant importance in interpreting gravity data and discerning geological features accurately. There are several methods for separating the regional and residual components, including filtering, the polynomial fitting method and the upward continuation method.
The Bouguer anomaly map of the study area (Fig. 2b) displays values ranging from −50 to 2 mGal, exhibiting a discernible SE to NW trend. This trend was approximated by fitting a linear polynomial surface (plan) to the Bouguer anomaly map. Subsequently, the identified trend was subtracted to derive residual anomalies, as illustrated in Fig. 3(a). Moreover, to enhance shallow sources, we computed the vertical derivative (Fig. 3b).
Residual anomalies obtained by subtraction of a trend estimated by fitting a first-degree polynomial surface to the Bouguer anomaly map (a). Vertical derivative gravity map (b). Both maps are presented in geographic coordinates and superimposed on the digital elevation model (DEM).
3.2.2. Edges and lineaments analysis
High gravity gradients often indicate the edges of rock units, geological contacts or faults. There are several techniques for edge detection, including: horizontal gradient (equation 1) maxima (Blakely & Simpson 1986), analytic signal (equation 2) maxima (Roest et al. 1992), Tilt Angle (equation 3) (Miller & Singh 1994), multiscale analysis (Boschetti et al. 1999), improved logistic filter (ILF, equation 4) (Pham et al. 2021). These techniques are based on the derivatives of the field in the x, y and/or z directions, which are commonly used for gradient computation or analytic signal.
To identify maxima, the localization process includes comparing a central point with its eight closest neighbours in four directions, all within a 3 × 3 data grid. If the central point has the highest value among its neighbours in all four directions, then a polynomial of second degree is fitted. The maximum of the parabolic curve can be computed analytically, which gives a more accurate estimate of the local maximum than simply comparing the central point with its eight neighbours. If the maximum of the parabolic curve is close to the central point, then it is considered a local maximum.
The tilt angle filter is computed as the arctangent of the ratio of the first vertical to the combined horizontal derivative of the gravity field. It varies between −90° and +90°, with positive values indicating a positive tilt (upward bending) of the gravity field and negative values indicating a negative tilt (downward bending). The zero contours of the tilt angle filter indicate the locations of the source-body contact, where the density contrast is maximum. Therefore, the zero contours of the tilt angle filter can be used to delineate the boundaries of the subsurface density anomalies.
Archibald et al. (1999) introduced a multiscale analysis method of edges. This tool uses wavelet analysis for edge detection in potential field data and can automatically detect gradients or edges at different scales of resolution.
In the presence of a vertical contact between rocks with abrupt changes in density, the gravity field shows a lower value above areas with lower density rocks and a higher value above regions with higher density rocks. Typically, the inflection point of the gravity field is located just above this contact. Moreover, the maxima of the horizontal gradient of gravity anomalies helps in locating contacts. In the case of an inclined contact, although the maxima of the horizontal gradient remains near the contact, it shifts in the down-dip direction (Archibald et al. 1999).
Khattach et al. (2004) utilized a combination of upward continuation and the horizontal gradient method to determine the dip direction of a geological contact. To achieve this, they identified the maxima of the horizontal gradient across multiple heights of upward-continued maps at various elevations. In the context of a dipping contact, it's important to note that the maxima from various levels maps will gradually shift downdip as the continuation height increases. Linear boundaries can form along geological contacts such as faults or flexure. Circular shapes in geology can be indicative of various features, including plutons or domes.
3.2.3. Depth estimation
To estimate depth, two methods were employed: Spectral analysis (Spector & Grant 1970) and Euler Deconvolution (eq. 5) (Reid et al. 1990). For this study, the located 3-D Euler Deconvolution Method (Daniel & Kingsley 2020) was utilized. This method involves selecting data points within a square window centred on the peak of the analytic signal amplitude, yielded more precise and accurate depth estimations compared to the standard 3-D Euler deconvolution. In this study, the following parameters were used: (i) structural index of 0.5, (ii) maximum percentage depth tolerance of 15 per cent and (iii) window size of 11 × 11. The results are presented in Fig. 6(a). The Euler Deconvolution equation is given as:
where (x0, y0, z0) are the coordinates of the magnetic source, F is the total field measured at (x, y, z). B represents the base level. N is the structural index, which relates to the geometry of the source.
4 RESULTS
4.1. Gravity results
Residual and vertical derivative maps (Figs 3a and b) reveal gravity lows and highs. The gravity highs are associated with high-density rocks (such as limestone and Paleozoic rocks), located in elevated topography areas, including Beni-Znassen (to the north), Jbel Hamra (to the south) and Jbel Megrez and Harraza (to the Southwest). In contrast, the gravity lows are observed in the plains, where borehole data indicates the accumulation of low-density rocks (primarily marl) within depressions.
The maps highlight of two depressions elongated in the NE-SW direction: the Beni Drar depression in the north and the Oujda depression in the south which extends into the Beni Oukil plain. These depressions are separated by a boundary or ridge characterized by relatively high gravity values. This boundary is likely associated with the uplifted Jurassic limestone, which may result from geological processes such as tectonic activity, including folding or faulting. Strong gravity gradients, which suggest the presence of discontinuities such as faults, demarcate these areas.
To delineate these geological features three methods were applied: Tilt Angle method (Fig. 4a), Improved Logistic filter (Figs 4b and c) and the horizontal gradient combined with the upward continuation method (Fig. 5). The tilt angle map shows values ranging from −90° to +90° with edges/or boundaries identified at angle zero. Negative values correspond to sources with low density, while positive values indicate sources with higher density. This method effectively defines the geological unit boundaries in the study area. The IL filter (Fig. 4b) enhances the visibility of sharper lineaments, revealing numerous gravity lineaments that were previously obscured by Neogene-Quaternary sediments in the plains. These lineaments predominantly trend NE-SW to ENE-WSW.
(a) Tilt angle map displaying values ranging from −90° to +90°, with edges/boundaries marked at angle zero (black line). (b) Logistic filter map highlighting sharp gravity lineaments (faults), where black dots represent maxima. (c) Impact of faults delimiting the Oujda depression to the south: Three boreholes (see locations in b) reveal a vertical displacement of the Jurassic carbonate exceeding 400 m. The gravity profile indicates that these carbonate rocks deepen as the gravity anomaly decreases.
Multiscale analysis of geological contacts: maxima (coloured dots) of the horizontal gradient of upward continued maps (residual and vertical derivative maps) at various heights are superimposed on the vertical derivative gravity map. Lines F, K, N and TR: correspond to resistivity profiles interpreted below Figs 7–10. Line S: represents the gravity profile used in spectral analysis (Fig. 6b).
Fig. 4(c) illustrates the influence of faults that define the southern boundary of the Oujda depression (south of Oujda city). Borehole data confirm these faults have caused significance vertical displacement of the Jurassic carbonate. Over a distance of 1 km, the carbonate has subsided by more than 400 m, providing clear evidence of fault activity and its impact on the region's geological structure.
To assess the significance and dip of these geological features, the multiscale analysis method was applied to both the residual map (to identify major lineaments) and to the vertical derivative map (to detect minor ones) (Fig. 5). The superimposition of horizontal gradient maxima, obtained at various levels of upward continuation of the residual and vertical gradient maps, delineated lineaments interpreted as faults, with generally exhibit a subvertical dip. The faults marking the boundaries of the depressions are clearly outlined. Two principal fault orientations identified: (i) N070°–085° and (ii) N040°–050° (more predominant). Less prominent faults strike in the N120°–140° direction. The primary fault systems compartmentalize the bedrock into buried depressions (grabens) and ridges (horsts), which align with the regional tectonic phases observed in northeastern Morocco. Given this block configuration, groundwater hydrodynamic flow may be discontinuous, as these faults can function as flow barriers, preferential flow paths or disrupt natural flow patterns, potentially leading to high-pressure discharge from aquifers and the formation of springs.
Recently, El Gout et al. (2024) advanced the understanding of the deep structure of the Jbel Hamra aquifer (south of Oujda City) by examining the relationship between the gravity field and the depth of the deep carbonate aquifer. Their study that: (i) gravity highs indicate the uplift of Jurassic carbonate and Palaeozoic rocks, and (ii) gravity lows suggest that the carbonate platform is buried beneath a substantial sedimentary cover (Miocene marls). Applying this approach can help identify optimal sites for new boreholes and provide critical insights for guiding sustainable groundwater resource management.
The estimated depth values to obtained using the Euler Deconvolution method are presented in Fig. 6(a): 27 per cent of solutions are less than 500 m, 31 per cent range between 500 and 1000 m and 42 per cent between 1000 and 2687 m.
Estimation of sources depths: (a) Solutions of Euler deconvolution superimposed on the vertical gradient map; (b) Radially Average Power Spectrum analysis showing two linear segments with corresponding depths; (c) Power spectrum of profile S (see location in Fig. 5).
These depth estimations align with the results from radially average Power Spectrum analysis (Fig. 6b). Specially, the first linear segment corresponding to deep sources, indicates gives a depth of 2663 m. Additionally, the power spectrum of profile S, which crosses the study area (Fig. 6c; see location in Fig. 5), suggests a depth of 2338 m. These findings provide valuable insights into the thickness of sedimentary deposits within the depressions.
4.2. Resistivity results
Four profiles, F, K, N and TR are presented to illustrate the resistivity characteristics across different zones of the Angad Bassin (Fig. 2a). Profiles F, K and N are N-S in the direction of SEVs. Profile F represents the western zone, while profile K corresponds to the central part, and profile N pertains to the eastern part of the plain. The fourth profile (TR), oriented northwest-southeast (NW-SE) and serves as a transverse cross-section perpendicular to the direction of structures. These profiles provide valuable insights into the subsurface resistivity variations within their respective regions. Figs 7, 8, 9 and 10 display the 2-D resistivity sections derived from the inversion of SEV data. These sections visually represent the resistivity distribution in the subsurface, offering a comprehensive understanding of the geological formations in the Angad Plain. By analysing these resistivity sections, valuable information about the subsurface lithology and structural variations can be obtained. The resistivity values shown in the figures reflect the electrical properties of the different geological formations in the region. The observed resistivity variations highlight the heterogeneity of the subsurface materials, including marl, marly limestone, limestone, sandstones, basalt and other formations.
Integrated Geophysical Analysis of Electrical Profile F: (a) VES Locations and Measurements, (b) Residual Anomaly and Vertical Gradient Gravity Profiles, (c) Resistivity Section (Electrode Spacing: 100 m), (d) Horizontal Gradient Magnitude of Residual Anomaly (HGMR) and Vertical Gradient Gravity (HGMVG), Improved Logistic Filter, (e) Resistivity Section (Electrode Spacing: 50 m). Dashed lines on resistivity sections indicate interpreted faults, black lines correspond to faults detected by HGMR, while red lines indicate faults detected by HGMVG.
High-Resolution Resistivity Section (Electrode Spacing: 20 m), focusing on the portion between x = 10 km and x = 16 km of profile F. The figure includes: (top) measured apparent resistivity, (middle) calculated apparent resistivity and (bottom) inverse resistivity model. Calibration using data from Borehole N°1225 indicates a strong correlation between resistivity layers and lithological layers.
Integrated Geophysical Analysis of Electrical Profile K: (a) VES Locations and Measurements, (b) Residual Anomaly and Vertical Gradient Gravity Profiles, (c) Resistivity Section (Electrode Spacing: 100 m), (d) Horizontal Gradient Magnitude of Residual Anomaly (HGMR) and Vertical Gradient Gravity (HGMVG), (e) Resistivity Section (Electrode Spacing: 50 m, depth: 500 m), (f) Calibration using Borehole data (Borehole N°149), (g) Vertical derivative gravity map overlaying digital elevation model and the resistivity section.
Integrated Geophysical Analysis of Electrical Profile N: (a) VES Locations and Measurements, (b) Residual Anomaly and Vertical Gradient Gravity Profiles, (c) Resistivity Section (Electrode Spacing: 100 m), (d) Horizontal Gradient Magnitude of Residual Anomaly (HGMR) and Vertical Gradient Gravity (HGMVG), (e) Resistivity Section (Electrode Spacing: 50 m, Depth: 500 m).
4.2.1. Resistivity section F
Section F (Fig. 7) spans 24 km and was generated using data from 22 VESs. The maximum electrode spacing (AB/2) used in the VESs was 3000 m, allowing for investigations down to a depth of 1200 m. Three different electrode spacings were employed to generate resistivity sections at varying depths and resolutions. For overall subsurface image, an electrode spacing of 100 m was used (Fig. 7c). To obtain a more detailed view of the subsurface down to 400 m, a smaller electrode spacing of 50 m was employed (Fig. 7e). Additionally, for a specific segment between x = 10 km and x = 16 km, an even finer electrode spacing of 20 m was utilized (Fig. 8). This finer electrode spacing allows for a more focused investigation of features at a smaller scale.
The resistivity values in the Section F vary significantly, ranging from 2 Ωm to over 4000 Ωm. This wide range indicates of the heterogeneity of the geological formations present in the area. Borehole 1225, located at x = 13 000 m with a total depth of 275 m, was used for resistivity calibration (Fig. 8). The first 110 m exhibit relatively low resistivity values (<40 Ωm), attributed to Quaternary sediments, specifically silt deposits. Between 100 m and the borehole bottom at 275 m, resistivity values increase significantly, exceeding 200 Ωm. This highly resistive layer corresponds to limestone and marly limestone of the Upper Jurassic, which outcrops in the western region of the plain (Jbel Megrez). A detail comparison between resistivity values and known lithologies demonstrates a strong correlation, with observed resistivity ranges aligning well with expected values for the identified rock types. This consistency reinforces the reliability of the resistivity measurements. Additionally, the thickness values of the identified layers obtained from the model closely match those derived from the borehole data, further validating the method used.
Fig. 7(c) highlights the geometry of the northern depression (from x = 14 000 to 24 000 m), revealed by gravity data. This geological feature is characterized by highly conductive rocks with resistivity values bellow 40 Ωm, closely associated with sedimentary deposits, specifically Miocene marl. The depth of this formation exceeds 1200 m, indicating a significant geological extent.
Figs 7(b) and (c) show a strong correlation between gravity and electrical data, revealing a clear relationship between these two geophysical parameters. The observations indicate that areas with high gravity values correspond to high resistivity materials, while areas with low gravity values are associated with materials of lower resistivity materials. This correlation provides valuable insights for subsurface characterization, allowing for the identification of regions with contrasting material properties based on their gravity and resistivity signatures. Such information is essential for geological and geophysical investigations, including groundwater exploration, mineral prospecting and structural mapping.
Figs 7(d), (e) and 8 further confirm the correlation between gravity and electrical data, reinforcing the relationship between these geophysical parameters. The maxima observed in the Logistic Filter, the Horizontal Gradient Magnitude of Residual Anomaly (HGMR) and Vertical Gradient Gravity (HGMVG) indicate the presence of faults, as marked by the purple, black and red lines. These faults are evident in the resistivity sections as distinct contrasts in resistivity values. Faults create boundaries between different lithological units, leading to variations in resistivity across the fault zones. The integration of these data sets enhances the understanding of the subsurface structural architecture by delineating the extent and orientation of faults.
4.2.2. Resistivity section K
The ‘k’ section spans a distance of 21 km and was created using data from 21 VESs (Fig. 9). These VES measurements involved a maximum electrode spacing (AB/2) of 5000 m, enabling investigations down to a depth of 2000 m. Data from borehole N°149, with a total depth of 390 m and located at x = 4000 m, were used for resistivity calibration (Fig. 9f). The first 110 m of the borehole exhibit relatively high resistivity values (>200 Ωm), which are attributed to Quaternary rocks, specifically gravel, conglomerates, basalt and silt. Bellow 110 m to the bottom of the boreholes (390 m), the resistivity values decrease significantly to below 30 Ωm, indicating the presence of marl.
This resistivity section traverses the two depressions and reveals three distinct resistivity layers: (i) A discontinuous superficial resistive layer associated with Quaternary materials (Conglomerates, basalt, gravel, clay, silt), (ii) An intermediate conductive layer, with varying thicknesses extending over 1500 m, related to the Miocene marl, (iii) A highly resistive layer in the form of discontinuous blocks, which indicates the presence of significant faults responsible for the uplift of these blocks, corresponding to Jurassic limestone. This uplift of limestone occurs at the boundary separating the two depressions, where it is found at a depth of 750 m. Fig. 9(g) shows the strong correlation and complementarity between the gravity and resistivity data; further supporting the geological interpretations.
4.2.3. Resistivity section N
The section ‘N’ spans a distance of 16 km and was generated using data from 13 VESs (Fig. 10). These VES measurements included a maximum electrode spacing (AB/2) of 5000 m, enabling for investigations to a maximum depth of 2000 m.
The section illustrates the southern depression, where the top of the deep resistive layer (Jurassic limestone) is located at a depth of less than 200 m (Figs 10c and e) extends beyond 2000 m below x = 14 000 m. This layer exhibits a 10° dip towards the north, and based on this dip, the maximum depth is estimated to be around 2500 m. This estimate aligns with the values obtained from spectral analysis, which range from 2300 to 2600 m (Figs 6b and c). These results validate the effectiveness and reliability of the methods used in this study.
4.2.4. Resistivity section TR
The section TR, oriented SE-NW and transverse to the geological structures (Fig. 2a), was generated using 14 SEVs with a maximum electrode spacing (AB/2) of 5000 m, extending investigation capability to a maximum depth of 2000 m. This section (Fig. 11c) illustrates the continuity of the two depressions in depth. The boundary between these depressions, located at x = 17 000 m, is characterized by an uplift of the deep resistive layer, which corresponds to Jurassic limestone. The marl layer reaches a maximum thickness of 1800 m in the southern depression and 1200 m in the northern depression. Additionally, the application of the logistic filter, along with the HGMR and HGMVG, highlights multiple faults that play a key role in the structural of the study area (Fig. 11d).
Integrated Geophysical Analysis of Electrical Profile TR: (a) VES Locations and Measurements, (b) Residual Anomaly and Vertical Gradient Gravity Profiles, (c) Resistivity Section (Electrode Spacing: 100 m), (d) Horizontal Gradient Magnitude of Residual Anomaly (HGMR) and Vertical Gradient Gravity (HGMVG).
5 2-D GRAVITY MODELLING
The inversion of electrical surveys has produced resistivity models that effectively represent subsurface structures and geological variations. However, due to the limited number of boreholes available to validation, we aimed to complement these models by developing density models through gravimetric modelling. We present a 2-D gravimetric model (Fig. 12) corresponding to the northern depression along electrical profile F (Fig. 7). The model layers were defined based on the resistivity section and consist of three layers: (i) a surface layer of marl, (ii) an intermediate layer of marly limestone and (iii) a basal limestone layer. Density values are assigned to each layer to ensure accurate representation. In this study, the modelling process was performed using FASTGRAV software, which calculates the gravity response of subsurface geological structures. The observed gravity data and the calculated response are presented alongside the geological model for comparison.
2-D Gravity model of the northern depression superimposed on the resistivity section (Profile F).
Fig. 12 illustrates both the gravity response of the marl layer alone and that of the complete model. The negative gravity anomaly is primarily attributed to the marl filling within the depression. Additionally, the shape of the residual gravity anomaly curve closely reflects the morphology of the Jurassic roof structure, further supporting the validity of the model.
The results of the gravity modelling corroborate and further validate the findings from the resistivity section inversion, reinforcing the consistency and reliability of the subsurface geological interpretations.
6 SUMMARY AND DISCUSSION
Gravity data reveals variations in low and high values, which can be attributed to the 3-D distribution of rock density. Two distinct gravity lows separated by a high, are observed in the plain. Analysis of resistivity and geological data strongly indicates that these lows correspond to depressions filled with relatively low-density rocks, particularly Miocene marl. In contrast, gravity highs are associated with uplifted areas composed of relatively dense rocks, including Jurassic limestone and dolomite. While the gravimetric maps clearly distinguish the two depressions (Figs 3 and 4a), resistivity sections suggest continuity at depth between these entities in the northeastern part of the study area (Figs 7 and 11). The positive gravity anomaly separating these two zones indicates the uplift of the Jurassic basement at depth (Figs 4a and 9). The results demonstrate a strong correlation between gravity and electrical data, revealing a notable relationship between these two geophysical parameters. Areas with high gravity values correspond to materials with high resistivity, such as carbonates, while areas with low gravity values are associated with low resistivity, such as Miocene marl.
The boundaries of these geological structures exhibit strong gravity gradients, serving as key structural indicators. The application of multiple edge techniques enabled the identification of the major faults responsible for structuring of the study area. In the plain, these faults are concealed beneath recent sedimentary deposits. The results confirm the complementary nature of the methods used in edge detection, including the HGMR, Vertical Gradient Gravity (HGMVG) and the Improved Logistic Filter (Figs 3 and 4b). The latter method offers high-resolution and accurate results; however, in some cases, certain features picks are identified by the other two methods are not detected. Depth estimation from gravity data (Fig. 6), derived using both the Euler method (with 42 per cent of sources ranging from 1000 to 2680 m) and the power spectrum technique (with deep sources depths of 2300 and 2660 m), and correlate well with the depth of the Jurassic basement depth inferred from resistivity sections in the depressions (Figs 7–11). These depth estimates cannot be directly confirmed by drilling data, as the deepest well in the area, located near the northeastern end of resistivity profile N, reaches only 700 m. However, this well confirms that at this depth, Miocene marl is still present, and the Jurassic basement has not yet been reached. To further validate these findings, gravimetric modelling was conducted, reinforcing the results obtained from the resistivity section inversion.
The findings of this study are valuable for deep aquifer prospecting, particularly within the Jurassic limestone and dolomite formations, which are currently exploited in the Jbel Hamra (south) and Beni-Znassen (north) mountains. The groundwater table in these areas has been experiencing a significant decline (≈3 m yr−1 in Jbel Hamra) due to intensive extraction and recurrent droughts. This issue affects a transboundary aquifer shared between Morocco and Algeria, highlighting the need for enhanced cooperation between the two countries to ensure sustainable management. This study presents a comprehensive depiction of the geological structure of Jurassic carbonates within the Angad Basin. These formations, which outcrop in the surrounding mountains, extend to depths exceeding 2000 m within the depressions. The results offer crucial insights into the complex subsurface architecture and distribution of this carbonate formation, supporting future exploration and sustainable management of this deep aquifer.
7 CONCLUSION
This study investigates the structural framework of the Angad Basin by integrating gravity and electrical data. Our approach employs advanced techniques to enhance gravity anomalies and improve subsurface imagining, including regional-residual separation, edge detection methods, depth estimation techniques and modelling. The resulting structural map reveals well-developed NE-SW and ENE-WSW fault systems, represented by several regionally significant faults, while the NW-SE faults are less dominant. Additionally, 2-D resistivity models derived from VESs effectively illustrate both vertical and lateral variations in subsurface resistivity. The geological data further validate the effectiveness of these geological methods in structural analysis.
Integrated analysis of resistivity, gravity and geological data confirms the presence of two depressions: the southern, filled with approximately 2000 m of Miocene marl, and the northern, with about 1200 m. Gravity modelling results align with resistivity section inversion findings, demonstrating that the residual gravity anomaly curve closely mirrors the morphology of the Jurassic roof structure. This structure is primarily influenced by the reactivation of Hercynian faults, which control the region's major fault systems.
This study highlights the strong complementarity between gravity and electrical data, demonstrating their effectiveness in subsurface geological investigations. Their integration provides a robust approach for understanding of subsurface characteristics, improving geological interpretations and supporting resource exploration, particularly for deep aquifers prospecting. The findings offer valuable insights into the complex subsurface architecture and distribution of the Jurassic carbonate formation, serving as a foundation for future groundwater exploration.
In summary, integrating previous geophysical investigations with advanced interpretation techniques provides a powerful, cost-effective alternative to extensive field surveys while maximizing the value of existing data. This approach enhances geological understanding, minimizes duplication of efforts and promotes collaboration research in groundwater exploration and environmental studies.
ACKNOWLEDGMENTS
We thank Dr Kosuke Heki, Editor of Geophysical Journal International, Dr Louise Alexander, Assistant Editor and the reviewers Dr Oyeyemi Kehinde D. and Dr Alison Kirkby, for their valuable comments and suggestions, which significantly contributed to improving our manuscript.
We sincerely thank Pr. Abioui Mohamed from the Department of Geology, Faculty of Sciences, Ibn Zohr University, Agadir, Morocco, for his invaluable assistance, which greatly contributed to improving the English quality of our manuscript.
Authors' contribution: DK: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Supervision, Validation, Visualization, Writing—original draft. REG: Data curation, Investigation, Methodology, Visualization, Writing—review and editing. SZ: Data curation, Investigation, Methodology, Visualization, Writing. AN: Investigation, Visualization, Validation. NN: Investigation, Visualization, Validation. MB: Data curation, Investigation, Visualization.
DATA AVAILABILITY
The gravity data is supported by the Ministry of Energy and Mines, while the electrical data are provided by the Hydraulic Basin Agency of Moulouya. Unfortunately, the data underlying this article cannot be shared.
