1 Introduction
The Moroccan basement holds more than three quarters of the world phosphate reserves (Bouda and Salvan, 1971). The reserves exist in four sedimentary basins, and the most important is the Ouled Abdoun, in central Morocco (Fig. 1A and B). In this basin, the sedimentary cover is about 2000 m thick, has a tabular structure, includes series that range from Permo-Trias to Quaternary (Choubert and Salvan, 1976), and lays unconformably on the Hercynian basement of the western part of the Moroccan Meseta (Michard, 1976; Zouhri et al., 2008). The Ouled Abdoun sedimentary cover is circumscribed in the North and the West by the Paleozoic outcrops of the Meseta, and to the southeast by the Jurassic escarpments of the Middle and High Atlas.
A. Geological map of the study area (Saadi et al. 1965). (a1) Paleozoic (schists and quartzites), (a2) Triassic (red clays and basalt), (a3) Jurassic (c: limestones and dolomites, d: diorite and gabbro), (a4) Cretaceous (red detrital facies), (a5) Tertiary (limestone, marl and phosphatic sands), (a6) Quaternary (alluvium), (a7) fault, (a8) synclinal axis, (a9) anticlinal axis, (a10) seismic line, (a11) interpreted seismic profile, (a12) borehole, (a13) river, and (a14) phosphate deposit. B. Synthetic stratigraphic column of the study area. (b1) Schists and quartzites, (b2) Conglomerate, (b3) red and green slightly anhydritic clays, (b4) doleritic basalts, (b5) anhydride, (b6) dolomite, (b7) marly dolomite, (b8) marls, (b9) gypsum marl, (b10) lumachellic limestones, (b11) dolomite with silex and phosphate debris, (b12) sandy limestone, polygenic conglomerates with clay cement, sandy marl, and (b13) cemented calcareous. C. Stratigraphic column of the phosphatic series of the Ouled Abdou basin. (c1) Marls with gypsum, (c2) phosphatic marls, (c3) marls, (c4) phosphatic limestones, (c5) phosphatic sand, (c6) discontinuous silex levels, (c7) marl and silex, (c8) Thersities slab Gharbi, 1998; Hsissou et al., 1996.
The phosphatic series form the upper part of the sedimentary cover of the Ouled Abdoun basin. Its age ranges from Maastrichtian to Lutetian (Azmany et al., 1986; Boujo, 1976). This series crops out in the northwestern part of the basin, where the phosphate layers are thus exploited in open-pit mines. By contrast, in the South and Southeast, the phosphatic series is covered by the Miocene, Pliocene, and Quaternary continental deposits of the Tadla Plain. Should the phosphate be exploited in these areas, exploration and exploitation methods would thus need to be adapted.
Our study aims to document the distribution and structure of the phosphate series under the continental deposits of the Tadla Plain. Based on the analysis of available gravity, seismic and borehole data, we recover the distribution, depth and overall structure of the phosphatic layer. Doing so, we define the areas with the thickest sedimentary filling and determine the depth of the buried phosphate layer. We anticipate that our results will be useful to guide the future exploration and exploitation of the buried phosphate in the Ouled Abdoun basin.
2 Geological context
The Ouled Abdoun basin belongs to the Moroccan Meseta. It extends over more than 10,000 km2 and is bounded to the north by the Hercynian massif of central Morocco, to the south and east by the Jurassic escarpments of the Middle and High Atlas, and to the west by the Paleozoic outcrops of the Rehamna. Geomorphologically, the basin includes two distinct units, the phosphate plateau and the Tadla Plain (Fig. 1A). Its importance is due not only to its large surface, but also to its high content of marine phosphates (Boujo, 1976; Moutaouakil and Giresse, 1993; Sheldon and Davidson, 1989). Among the Moroccan phosphate basins, it is the one that is known (1917) and exploited (1921) since the longest time.
The sedimentary series of the Ouled Abdoun is traditionally described (Azmany, 1979, Moutaouakil and Giresse, 1993; Sheldon and Davidson, 1989) as comprising five large lithological units. Superimposed from bottom to top (Fig. 1B), these include (i) the Cenomanian, a complex series with alternating marls and limestones, (ii) a Turonian series dominated by dolomites and massive limestones with rare marly intercalations affected by karstic phenomena, (iii) a so-called Senonian series, but probably of Coniacian to Campanian age, which is essentially marly, (iv) a Maastrichtian to Ypresian series made of phosphate beds that alternate with sterile layers, and (v) a Lutetian series mainly composed of marine biocalcarenites that are lithostratigraphically grouped as the Thersities slab (marine gastropods of Hemithersitea and Thersitea sp.; Salvan, 1954).
From a structural point of view, the Ouled Abdoun sedimentary cover has been hardly deformed during the subsequent Upper Cretaceous and Cenozoic tectonic phases (Bouda and Salvan, 1971). Although the edges of the basin have all been uplifted and consecutively eroded during the Alpine orogeny, no major tilting affects the cuesta reliefs in the region. However, the main talwegs that cut the phosphate plateau and drain towards the Oum Rabia River have the same orientations as the main faults that bound the Permo-Triassic grabens exposed in the North and the Northeast of the basin. Ranging between ∼ N20° and N40°, this orientation is also that of the faults that controlled the Jurassic sedimentation and Cenozoic uplift of the Middle Atlas (El Hammichi et al. 2002; Fedan, 1989). Altogether these suggest that ancient NE-striking faults exist under the sedimentary series of the Ouled Abdoun and that some of them might have been reactivated during the Atlas orogeny (Bouda and Salvan, 1971). The reactivated faults might have influenced the Thersities slab structures, and controlled the current hydrographic network. It is therefore possible that some of these inherited faults have controlled the phosphatic series distribution and thickness beneath the continental deposits of the Tadla Plain.
The phosphatic series (Fig. 1C) is formed by alternating layers of sandy phosphate, sterile limestone, and marl-limestone levels with an average thickness of 50 m in the deposits that are currently exploited at the ground surface (Kchikach et al., 2002; 2006). The series begins with phosphate marls (c2) and calcareous levels very rich in bone debris known as the Maastrichtian bone-bed limestones (c3) (Noubhani and Cappetta, 1995; Suberbiola-Pereda et al., 2003). This geological stage marks the beginning of a phosphatogenesis that reached its maximum during the following stages. Above, exists a succession of phosphatic limestone (c4) and sandy phosphate (c5) levels. Montian (Selandian) is represented by sandy phosphates (c5) overlaid by limestone with coprolite and silex nodules (c6) that constitute a benchmark level for mining. The latter is overlaid by alternating regular beds of marly and phosphatic limestones, levels of coarse-grained sandy phosphates, continuous silex horizons (c7), and, sometimes, silto-pelitic levels of Thanetian to Ypresian age (Azmany et al., 1986). The Lutetian is represented at its base by alternating levels of silex, marl, and limestone where the last terms of phosphatogenesis are intercalated. The summit corresponds to the thickest carbonate levels, rich in gastropods, called the Thersites slab (c8) (Salvan, 1954).
The phosphatic series overlays a thick Senonian marl formation (Boujo, 1976; El Assel et al., 2013). In the majority of the Tethysian and Peritethysian basins, the Senonian is generally regressive, whereas the Maastrichtian is rather transgressive (Hardenbol et al., 1998; Salvan and Farkhany, 1982). The Senonian in the Ouled Abdoun basin is commonly composed of very shallow and confined marine environment deposits. Evaporitic deposits, especially gypsum, are frequent (El Assel et al., 2013). Choubert and Salvan (1976) described that, in some places of the Ouled Abdoun basin, a thick Senonian series (more than 250 m) is formed exclusively of marls. The Senonian marls are overlaid by Maastrichtian phosphatic limestones that constitute the base of the phosphatic series. The vertical seismic profiles performed in the Tadla Plain (by oil companies) show a strong contrast of acoustic impedance between the Senonian marls and the Maastrichtian phosphatic limestones. This acoustic contrast forms a significant reflector generally clear on seismic sections.
The phosphatic layers exploited in the Khouribga open-pit mines continue below the Tadla Plain. The plain extends between the phosphate plateau to the north and the High Atlas chain to the south. It occupies an area of ∼3600 km2 in the Ouled Abdoun basin. The Neogene and Quaternary filling of the plain is essentially formed by lacustrine limestones overlain by continental deposits (conglomerates and red clays). The works of El Antaki and El Boustani (1991) and El Mansouri (1993), based on borehole data, showed that the thickness of the sedimentary filling above the phosphatic series in the plain may be up to 300 m.
3 Material and methods
To perform this study, we used the geological, gravity, and seismic data provided by the National Office of Hydrocarbons and Mines (ONHYM) and the Moroccan Ministry of Energy and Mines (MEM) (CGG, 1985; Van Den Bosch, 1981). The data include maps, time migrated seismic sections, and borehole columns.
3.1 Gravity data
We used the gravity data available in the MEM Moroccan geological department, namely a Bouguer anomaly map (d = 2.67 g cm−3, Fig. 2A) established earlier from 3312 ground measurements (Van Den Bosch, 1981). From this map, we first calculated the residual anomaly map by subtracting the regional effect (Fig. 2B). The first trials to eliminate this effect using low-pass filtering or upward continuation did not result in a satisfactory outcome. The best result was obtained using the polynomial regression method (Draper and Smith, 1998).
Gravity maps of the study area: (A) Bouguer anomaly, and (B) residual anomaly. (1) Paleozoic basement outcrop, (2) faults, (3) synclinal axis, (4) anticlinal axis, (5) gravity anomaly (p: positive, n: negative), (6) Tadla Plain boundary, (7) borehole.
To highlight the geological contacts in the residual gravity map, we processed it using a contact analysis method based on two techniques: horizontal gradient and upward continuation using the Geosoft mapping and processing system (Geosoft, 1994).
A (sub)vertical discontinuity separating two formations with different densities generates a significant anomaly in a gravity map that corresponds to a gravity change from lower values in the rocks with lower density to higher values in the denser rocks. A point of significant gravity change (inflection point) in a residual map is thus taken to locate the contact between two types of rocks. We thus sought to identify and locate the most abrupt changes of density in the residual map. For that, we computed the horizontal gradient of the residual map, as those emphasize the inflection points as maxima (Blakely and Simpson, 1986). We located these maxima using a code of “contact analysis” based on the Blakley and Simpson (1986) method that determines the peak points of the gradient amplitude. The code moves a nine-mesh (3 × 3) window on a regular grid of the horizontal gradient. At every step, the value at the central point of the window is compared to that in the eight neighboring cells. The repeating process allows extracting the maximum gradients. This method has been successfully applied in earlier works to identify and map faults and other discontinuities in gravity data (Everaerts and Mansy, 2001; Khattach et al., 2004). Here, we have used it to locate the most significant lateral changes in rock density. We hypothesize that such pronounced changes result from faults dissecting the area and juxtaposing rocks with different densities.
Under this assumption, we then tried to determine the direction of the fault (or contact) dips. This can be done through the upward continuation of the residual gravity map at different elevations (chosen here as 2000, 5000, and 7500 m), and the determination of the horizontal gradient maxima at each elevation. The spatial migration of these maxima as the upward altitude increases indicates the fault dip direction (Fig. 3) (Archibald et al., 1999; Khattach et al., 2004). The values of the upward continuation elevations were selected arbitrarily, yet to satisfy two conditions: (1) to be sufficiently far apart to show the displacements of the horizontal gradient maxima, and (2) to be sufficiently close to prevent that the effect of the gravity contacts (i.e. the faults) be attenuated.
Superimposition of the residual gravity map (background color) on the horizontal gradient maxima calculated for this map (1) and for its upward continuation at altitudes of 2000 (2), 5000, (3) and 7500 m (4). (5) Dip direction, (6) Tadla Plain boundary, and (7) Paleozoic basement outcrop, (8) interpreted seismic profile.
3.2 Seismic data
The seismic data used in this study correspond to 18 seismic reflection profiles, varying in length between 15 and 90 km, implanted following NW–SE and NE–SW directions (Fig. 1A). These profiles were recorded in 1985 by the General Geophysical Company (CGG) on behalf of the ONHYM, in the framework of an oil exploration program conducted in the Ouled Abdoun basin. The seismic survey was performed using a recording time of 2 s and a common 24-order mid-point, which is an acquisition method in multichannel seismic reflection. With this method, seismic traces, corresponding to rays that share the same subsurface reflection point, are collected and stacked in order to improve the signal-to-noise ratio. This is needed because reflections are generally weak and stacking adds a number of traces where the signal is the same but the random noise is different. This results in enhancing the signal and attenuating the noise. The seismic data are presented as time migrated sections at a scale of 1:20,000. A crude approximation of the density of the rocks suggests that the depths reached in the seismic profiles range between 2700 and 3300 m. The data also include a set of 20 borehole columns, distributed across the entire Tadla Plain. Most of these oil drilling reached the Paleozoic basement (CGG, 1985), at depths of 225, 325, 1350 and 1950 m, respectively, in boreholes FBS6, FBS3, DRZ1, and MKS1 (Fig. 2).
First, we examined the seismic sections to identify the reflectors corresponding to the base of the phosphatic series (corresponding to the top of the Senonian) and to the top of the Paleozoic basement (Fig. 4). The two target reflectors were generally clearly identified due to the strong acoustic impedance contrast between the Paleozoic basement (b1, Fig. 1B and , Fig. 5) and the cover (b2, Figs. 1B, 2 and 5), and between the thick Senonian marls (b9, Figs. 1B, 2 and 5) and the dolomitic phosphate limestones of the Maastrichtian that mark the base of the phosphatic series (b11, Figs. 1B, 3 and 5). Then, the seismic sections were calibrated using velocity laws obtained in vertical seismic profiles carried out in prior oil drillings in the Tadla Plain, and the calibrated seismic sections were compared to the borehole data. Our interpretation of the two target reflectors was thus controlled by the borehole columns that intersected these two interfaces. We then digitized the two reflectors on all seismic sections and this allowed us to build two isochrons maps of the base of the phosphatic series and of the Paleozoic top (maps established using the Geosoft mapping and processing system (Geosoft, 1994). We also carefully analyzed each seismic section, especially to identify the possible major faults affecting the sedimentary series, in particular in the Tadla Plain.
Example of seismic section showing the studied corresponding reflectors (top of the basement and top of the Senonian).
Seismic section interpretations along the P1, P2, and P3 profiles. (1) Paleozoic basement, (2) Ante-Eocene cover, (3) Post-Eocene cover, and (4) borehole. OR: Oum Rbia river. TWT: Two-way travel time.
4 Results and discussion
4.1 Gravity
The Bouguer anomaly map reveals that gravity follows an overall regional gradient with values increasing progressively from about −103 mGal in the Southeast to 7 mGal in the Northwest (Fig. 2A). The gradient seems maximum perpendicular to a general NE direction that is also that of the greatest elongation of the basin. The residual map highlights the gravity anomalies better (Fig. 2B). Several major positive (P1 to P10) and negative (N1 to N8) anomalies are revealed, which we interpret as reflecting the lateral variations of the underground rock density. The geological outcrops of the Paleozoic basement at the northeastern and southwestern ends of the study area coincide with the P1 and P5 positive gravity anomalies (Fig. 2B). This suggests that the major gravity anomalies likely reflect variations in the topography of the top of the Paleozoic basement; structural highs and lows in this basement would be marked by positive and negative gravity anomalies. The alignment of the P1, P2, and P3 positive anomalies thus form a major structural high that likely corresponds to an intumescence of the Hercynian basement along a NE–SW direction. The Paleozoic basement that crops out to the Northeast of the central Morocco Hercynian Massif would thus extend further to the southwest, forming a ridge deepening from northeast to southwest up to the town of Fkih Ben Salah. The geological data show that this ridge coincides with an important anticlinal structure in the “Atlas direction,” whose axis has been mapped at the surface between anomalies P1 and P2 (Verset, 1988) (Fig. 2B). Another major positive anomaly P4, thus a structural high, is aligned with the P1–P3 ridge. This suggests a significant decrease in the depth of the Paleozoic basement as one approaches the Hercynian Rehamna Massif (whose eastern limit is 5 km to the west of the P4 anomaly). The P5 and P6 anomalies in the south might be the signature of a third structural high in the Paleozoic basement. As a matter of fact, the Paleozoic basement crops out in the Jebilets massif at the P5 anomaly. We infer that the Paleozoic basement deepens towards the northeast (it is deeper at P6).
West of Khouribga, the positive anomalies P7, P8, P9, and P10 suggest that the Paleozoic basement is relatively shallow under the phosphate plateau. Outside these structural high zones, the Paleozoic basement seems to deepen significantly, as attested to by the negative anomalies in the Khouribga–Oued Zem zone (N7 and N8), between El Borouj and Fkih Ben Salah (N5 and N6), and to the North of the Jebilets massif (N4). The southern edge of the study area hosts dark blue to green negative anomalies, and includes four major zones of negative anomalies, N1, N2, N3, and N4. The Tadla Plain has been shown to be an asymmetrical depression of tectonic origin whose formation was related to the uplift of the High Atlas (Michard, 1976). Its deepest part, and hence major depot-center with thicker sediment content, is in the South, about where the N1 to N3 negative anomalies are observed.
The borehole data confirm the above interpretations. In the areas of positive anomalies, the Paleozoic basement is intercepted at shallow depth. For instance, in the anomaly zone P3, the basement is found at 225 m and 325 m in the boreholes FBS6 and FBS3, respectively (Fig. 2). By contrast, the Paleozoic basement is found at 1350 m (borehole DRZ1) and 1950 m (borehole MKS1) in the area of the N3 negative anomaly, confirming that the Paleozoic basement is deeper beneath the negative gravity anomalies.
In addition, the residual anomaly map shows that strong gravity gradients (≥0.89 mgal/km) separate some of the recognized structural lows and highs. These strong gradients may correspond to contacts or discontinuities, likely faults and/or flexures. We have analyzed these gradients following the methodology described in section 3 (Fig. 3). This analysis identifies contacts the majority of which trend NE–SW and dip towards the southeast. The northeastern trend of the identified discontinuities coincides with the overall structural orientation of the Atlas chain and with the major elongation of the Tadla Plain.
We interpret the identified NE structures as faults dipping to the southeast. Most are several tens of kilometers long, possibly up to ∼100 km. However, the gravity data resolution does not allow recovering either the fault lengths precisely or the precise geometry of their traces. Among these faults, F2 and F3 are likely major normal faults whose down-to-the-SE slip has contributed to the topographic variation revealed in the gravity anomalies (Fig. 3). F1 might be a third major normal fault system that crosses the study area in the North.
4.2 Seismic
While interpreting the seismic reflection time migrated sections, we focused on the phosphatic series structure and its deepening under the Neogene and Quaternary deposits of the Tadla Plain. First, we identified the two reflectors that correspond to the base of the phosphatic series and to the Paleozoic top (Fig. 5). The reflector of the latter could be easily followed along the different sections and from one section to another at their intersection.
Fig. 5 presents a 2D interpretation of the seismic sections along three ∼ NW–SE profiles that cross the entire Tadla Plain. Their location was chosen to best benefit from the available borehole data needed for the seismic data calibrating.
As synthesized in Fig. 5, the seismic sections reveal a progressive thickening of the post-Paleozoic sedimentary series from NW to SE. We interpret this thickening as mainly due to the vertical, normal motion across the F1, F2, and F3 major faults recognized in the gravity data; the normal slip on these faults has given a stairway structure to the sedimentary series, inducing their thickening from NW to SE.
The seismic sections also show a thickening of the sedimentary series from NE to SW in the Tadla Plain. This induces a deepening of the phosphatic series under the Neogene and the Quaternary Tadla Plain cover, along with its thickening along the NE–SW and NW–SE directions.
From the digitization of the two above-mentioned seismic reflectors, we established their isochron maps (Fig. 6A and B). These maps better reveal the overall distribution and structure of the sedimentary series under the Neogene and the Quaternary Tadla Plain cover. The isochron map of the phosphate layer base shows that it crops out or is shallow north of the El Borouj-Kasbat Tadla axis (Fig. 6A). South of this axis, the base of the phosphatic series progressively deepens to reach its maximum depth in the sedimentary deposit–center located west of Beni-Mellal. According to the borehole data, the bottom of the phosphatic series is located at a depth of ∼150 m near Fkih Ben Salah (Borehole BJ115, Fig. 1A) and at a depth of 130 m near the town of Tadla (Borehole BJ109, Fig. 1A). These results corroborate our interpretation of the seismic sections. The Paleozoic top isochron map shows roughly the same pattern. The Paleozoic basement is shallow to the North of the study area and also deepens progressively towards the South. The maximum depth of the basement is also found in the same area as for the phosphatic layer, i.e. in the sedimentary deposit–center in the South. The two reflectors thus show the same structure overall. This suggests that the F1, F2, and F3 faults that control the overall pattern of the Tadla depression are rooted in the Paleozoic basement and were reactivated during the subsequent sedimentary processes.
Isochron maps of the base of the phosphatic series (A) and the top of the Paleozoic basement (B).
We eventually built a structural map of the Tadla Plain. This map integrates the different structures recognized in this study and from available data (Fig. 7), and also shows the isochron map of the phosphatic layer base. Except for the Khouribga mining district area, where the phosphatic layers are exposed, the structural map allows a better understanding of the distribution and deepening of the phosphatic layers under the Neogene and the Quaternary Tadla Plain cover. The map also helps locating the areas where the phosphatic layers are deeper, or sustain pronounced depth changes across the F1, F2, and F3 normal faults. We anticipate that this new vision of the phosphatic layer distribution and structure might be important to guide the future exploitation of the phosphates, especially south of each of the identified normal faults.
Structural map of the phosphate plateau and the Tadla Plain. (1) Faults, (2) synclinal axis, (3) anticlinal axis, (4) Paleozoic basement outcrop, (5) isochron of the phosphatic series base, (6) interpreted structures from the gravity data, (7) interpreted major fault, (8) dipping direction.
5 Conclusion
In this study, we have provided a revised and refined description of the phosphatic series distribution and structure under the Neogene and the Quaternary cover of the Tadla Plain. We showed that, while this series crops out in the phosphate plateau where it is exploited in open-pit mines, the phosphatic layers are progressively deepening southeastwards from the plateau, through a series of steps induced by a few major NE-trending normal faults dipping southeast. The base of the phosphatic series would deepen to more than 500 m in depth near the Moroccan High Atlas chain in the Southeast.
This overall structuration implies an augmentation in the phosphate reserves from northwest to southeast. The isochron map of the base of the phosphatic series suggests that it is the deepest, and hence the phosphate layers thickest, south of the F2 major fault and west of the Beni-Mellal–Fkih Ben Salah axis (Fig. 7).
The present study thus suggests that, in the long term, the open-pit phosphate extraction method currently used in the Khouribga mining district will have to be replaced by other methods better suited to the more complex situation revealed here, where the phosphatic series is buried at significant depth in the plain. South of the NE–SW axis materialized by the F2 fault, the phosphate exploitation should be carried out either by underground mines or by large open quarries well dimensioned to allow the machines to access the lowest phosphate level. The Office Cherifian of the Phosphates group will soon carry a recognition program that might contribute to refine the present results.
Acknowledgments
This research was carried out as part of the project ID GEI-KCH 01/2017, funded by OCP (Office Cherifian of the Phosphates group, Morocco). Particular thanks are addressed to the Ministry of Energy and Mines and to the National Office of Hydrocarbons and Mines, Rabat, Morocco, for the data made available. The comments and suggestions of an anonymous reviewer and the associate editor Isabelle Manighetti were a welcome help to improve the manuscript.