Geophysical Prospecting for Groundwater Resources in Phosphate Deposits (Morocco)

: The Moroccan phosphate deposits are the largest in the world. Phosphatic layers are extracted in open-pit mines mainly in the sedimentary basins of Gantour and Ouled Abdoun in Central Morocco. The purpose of this study was to prospect and evaluate the water potential of aquifers incorporated in the phosphatic series using the following geophysical methods: Magnetic resonance sounding (MRS), electrical resistivity tomography (ERT), time-domain electromagnetics (TDEM), and frequency-domain electromagnetics (FDEM). The objective was, on the one hand, to contribute to the success of the drinking water supply program in rural areas around mining sites, and on the other hand, to delimit ﬂooded layers in the phosphatic series to predict the necessary mining design for their extraction. The use of geophysical methods made it possible to stratigraphically locate the most important aquifers of the phosphatic series. Their hydraulic parameters can be evaluated using the MRS method while the mapping of their recharge areas is possible through FDEM surveys. The results obtained in two selected experimental zones in the mining sites of Youssouﬁa and Khouribga are discussed in this paper. The application of the implemented approach to large phosphate mines is in progress in partnership with the mining industry.


Introduction
Phosphate is a valuable mineral resource for different industries. In Morocco and all over the world, phosphate rocks are being extracted for their phosphorus content. Phosphate is used to produce phosphoric acid (PA) that have various uses in our daily lives. First, PA is mainly used in agriculture for the production of soil fertilizers. It is also used as animal feed supplements. PA is also widely solicited in laboratories because it resists to oxidation, reduction, and evaporation. It has many other industrial applications such as soaps and detergents, water treatment, dentistry, etc.
Moroccan subsoil includes more than three-quarters of the world's phosphate reserves [1,2]. These reserves exist mainly in the sedimentary basins of Gantour and Ouled Abdoun in Central Morocco (Figure 1a,b). In these basins, the age of the tabular structure of sedimentary cover ranges from Permo-Triassic (251.9 My = million years ago) to Quaternary (2.58 My) [3]; it rests unconformably on the Hercynian basement of the western part of the Moroccan Meseta [4,5]. The phosphatic series investigated in this study corresponds to the upper part of the sedimentary cover. Its age extends from the Maastrichtian (72.1 My) to Lutetian (47.8 My) [6,7]. This series outcrops in the northern part of the aforementioned basins where the phosphatic layers are exploited in open-pit mines. To the south, this series is buried under the Neogene (23 My) and Quaternary continental deposits of the Bahira and Tadla plains (Figure 1a,b).
Minerals 2020, 10, x FOR PEER REVIEW 2 of 16 Morocco (Figure 1a,b). In these basins, the age of the tabular structure of sedimentary cover ranges from Permo-Triassic (251.9 My = million years ago) to Quaternary (2.58 My) [3]; it rests unconformably on the Hercynian basement of the western part of the Moroccan Meseta [4,5]. The phosphatic series investigated in this study corresponds to the upper part of the sedimentary cover. Its age extends from the Maastrichtian (72.1 My) to Lutetian (47.8 My) [6,7]. This series outcrops in the northern part of the aforementioned basins where the phosphatic layers are exploited in open-pit mines. To the south, this series is buried under the Neogene (23 My) and Quaternary continental deposits of the Bahira and Tadla plains (Figure 1a,b). The present study was conducted to characterize the aquifers of the phosphatic series and better understand their properties and characteristics. This study is part of an ambitious sustainable The present study was conducted to characterize the aquifers of the phosphatic series and better understand their properties and characteristics. This study is part of an ambitious sustainable development program led by the Office Chérifien des Phosphates (OCP) around phosphatic extraction sites. The drinking water supply of rural agglomerations in the vicinity of mining sites is among the main objectives of this program. In addition, OCP's geo-mining engineers are interested in determining the flooded areas in the phosphatic series which will allow them to predict all the required preparations to extract phosphates in such a hydrogeological context. This study combined the use of at least two different geophysical methods across two selected areas, namely Youssoufia and Khouribga mining sites (Figure 1a,b). Electrical resistivity tomography (ERT) and magnetic resonance sounding (MRS) surveys were carried out in the Youssoufia area ( Figure 1a); while time-domain electromagnetic (TDEM) and frequency-domain electromagnetic (FDEM) measurements were performed in the Khouribga site ( Figure 1b).
The treatment and interpretation of the acquired data made it possible to locate the main horizon of the aquifer contained in the phosphatic series and map on the surface conductive corridors that would constitute potential recharging zones.
The experimental research conducted at the two aforementioned sites allowed for the establishment of a scientific approach for larger hydro-geophysical studies to be conducted on the phosphatic series. It also made it possible to define appropriate devices and sampling steps for the used geophysical methods. The MRS method is a satisfactory tool to explore the aquifers and determine their hydrodynamic parameters in the phosphatic series context. FDEM surveys, which are very easy to implement in field, are the most efficient tool for mapping aquifer recharge zones in the geological context of our studied areas. The ERT and TDEM methods can also be used in areas where MRS data are affected by natural electromagnetic noise. We discuss, in this article, the approach adopted and the main obtained results.

Geological and Hydrogeological Context of the Study Areas
This research was conducted at two OCP mining sites in two distinct sedimentary basins: Gantour and Oulad Abdoun. The Youssoufia site is located in the Gantour Basin. In this basin, the sedimentary series ranges from the Triassic to Quaternary period [7]. It is surrounded from the north and south, respectively, by the Paleozoic outcrops of the Rhamna and Jebiletes Massifs (Figure 1a). The phosphatic series outcrops in the northern basin where they are exploited in open-pit mines in Youssoufia and Benguerir. In the southern basin, this series is buried below the Neogene and Quaternary deposits of the Bahira Plain. In the Gantour Basin, the phosphatic series is well known owing to its exploitation and borehole recognition ( Figure 1c). It starts with the Maastrichtian, mainly formed of fine sandy phosphates, phosphatic marls, and a clay-sandstone complex constituting an impermeable reference layer [7,10]. The Paleocene that is the lower epoch of the Paleogene, which includes the Montian (61.6 My) and the Thanetian (59.2 My), is mainly composed of alternating layers of uncemented phosphates, limestones and/or phosphatic marls, coprolites, and flint nodules. The Lower Eocene (Ypresian (56 My)) is commonly formed of phosphatic marls intercalated by phosphatic limestones, silexite, and siliceous marl [11,12]. Uncemented phosphatic layers become weak and often overlay impermeable clay layers, known as Ypresian clays by geo-mining engineers, in the Gantour Basin. Lutetian (47.8 My) strata mark the end of the phosphatic series. It is characterized by a fossiliferous marly limestone slab, known as the Thersitae slab [7,[13][14][15][16].
The second study area is in the Ouled Abdoun Basin situated approximately 26 km south of Khouribga. This basin includes the largest phosphate deposits in the world. It extends over more than 10,000 km 2 and is bounded on the north by the Hercynian Massif of Central Morocco, on the south and east by the Jurassic escarpments of the Middle and High Atlas, and on the west by the Paleozoic outcrops of the Rhamna [14]. Geomorphologically, this basin contains two distinct units: The phosphate plateau and the Tadla Plain (Figure 1b). The phosphatic series (Figure 1d) is formed of alternating layers of uncemented phosphate, sterile limestone, and marly limestone with an average thickness of approximately 50 m in the deposits under extraction [16,17]. It begins with phosphatic marls and limestone layers that are very rich in bone debris, known as bone-bed limestones, of Maastrichtian age [18,19]. This stage marks the beginning of a phosphatogenesis that reaches its maximum during the following stages. Above, successive layers of uncemented phosphates and phosphatic limestone occur. The Montian is represented by uncemented phosphates, overlain by limestone with coprolites and flint nodules that constitute a reference layer for mining extraction. The last is also overlain by alternating regular beds of marly and phosphatic limestones, coarse-grained uncemented phosphatic layers, continuous flint horizons and, sometimes, silto-pelitic layers from Thanetian to Ypresian in age [6]. The Lower Lutetian consists mainly of alternating flint layers and slightly phosphate marls and limestones that mark the last stages of phosphatogenesis. The phosphatic series is overlain by strong fossiliferous carbonated layers, rich in gastropods, known as the Thersitae Slab [20]. The phosphatic layers outcrop or are beneath a thin quaternary cover in the Northern Ouled Abdoun Basin. They are extracted in the open-pit mines of Khouribga, Daoui, Merah, and Sidi Chennane. To the south and southeast, the phosphatic series is buried beneath the Neogene and Quaternary continental and lacustrine deposits of the Tadla Plain [21].
Hydrogeologically, the main groundwater aquifers in the sedimentary series of the Gantour and Ouled Abdoun basins correspond to karstic Turonian limestones, uncemented phosphatic horizons, and fractured limestones that cover the phosphatic series [8,[21][22][23]. The Turonian aquifer is very deep (greater than 250 m) in the two aforementioned basins, and therefore less recognized and extracted. The aquifer corresponding to the Lutetian fractured limestones of the Thersitea slab is only important for water supply in areas where the phosphate series is buried under the Bahira and Tadla Plains. Along the perimeters of the open-pit mines in both the Gantour and Ouled Abdoun Basins, the piezometric level is generally below the Ypresian limestones.
This study is concerned with the aquifers in the uncemented phosphate horizons of the phosphatic series. Actually, this aquifer system is composed of at least eight superimposed layers of uncemented phosphates; each is overlaying a sterile impermeable interlayer of clay, marl, continuous flint bench, or siliceous phosphatic limestone. These aquifer horizons connect with each other via a network of fractures and microfaults affecting the phosphatic series. They are from bottom to top (Figure 1c,d) the following: (i) Two layers of Maastrichtian uncemented phosphates that constitute the most important aquifer of the Youssoufia and Benguerir mines; (ii) four layers of Montien-Thanetian uncemented phosphates that often overlie siliceous limestone layers. The thickness of these layers changes from one location to another in the studied basins [8]. Finally, the Ypresian phosphatic coarse layers overlie either plastic clays or siliceous marl. They correspond to the Ypresian strata indicated with blue background in Figure 1c.

Materials and Methods
To conduct the present hydrogeophysical study of the OCP mining sites, we chose two zones to perform experimental geophysical surveys. These zones were selected near extraction areas where geological and hydrogeological data are available to calibrate the results of the used geophysical methods. The first zone was located in the Gantour Basin at the "Panel 6" deposit according the OCP terminology (Figure 1a). The second is adjacent to the Sidi Chennane mining deposit in the Ouled Abdoun Basin (Figure 1b). During this experimental phase, we arbitrarily chose to test the MRS and ERT methods at the Youssoufia mine site, while the TDEM and FDEM methods were used at Khouribga. The same methods can be used at any of the OCP mining sites because of the similarity of the geological and hydrogeological contexts, as shown in Figure 1c,d; in addition, the two sites chosen to conduct the present study have the advantage of a low electromagnetic noise level. They are located far from any source of anthropogenic noise generated by human activities, such as power lines, radio antennae, electric fences, motors, pumps, buried conductors, etc.

ERT Profiles
ERT profiles provide two-dimensional (2D) images of subsurface resistivities. The implantation of a large number of electrodes, at regular intervals, along a rectilinear profile, enables one to obtain pseudo-sections of apparent resistivities. Using modern equipment, the acquisition is automatic by interrogating quadrupoles of increasing length, thus making it possible to obtain increasing measurement depths. A data inversion is then performed to obtain a quantitative measurement of the true electrical resistivity at real depths.
An ABEM Terrameter LS 12 with 81 electrodes was used and made available by the Andaluz Institute of Geophysics (AIG) at Granada University, Spain. We adopted an edge gradient acquisition device which allowed for a better assembly between electrodes by combining symmetrical and asymmetrical configurations. Such a configuration provides an image of subsurface resistivities at high resolution and with less artifacts, e.g., distortions of the inversion resistivity model caused by the relatively high noise contamination of the data [24]. Three ERT profiles 700 m long with 5-m electrode spacing were collected ( Figure 2). A total of 6531 apparent resistivity measurements were acquired. Before the inversion of the data, their quality control was thoroughly controlled. No automatic filtering was applied but the data were edited, and all bad datum readings were manually eliminated. This resulted in slight reduction of the total corrected datum point for each profile. The data inversion was performed using the AIG's Res2dinv software (V.358-Geotomo Software). The program uses a smoothness-constrained Gauss-Newton method to obtain a 2D resistivity model [25]. The process is automatic once the inversion parameters are established and a user-defined initial model is not needed. The iterative inversion process was set at a maximum of 10 error-iterations using the L2-norm (root mean square), offering an error, between the apparent resistivities in the field and those generated by the calculated models, less than 5%. The inversion model contained 15 layers and 1692 resistivity blocks. It was refined to half the electrode spacing because of the sensitivity of the gradient array to near surface variation. The sensitivity matrix was calculated using the incomplete Gauss-Newton method that is recommended for large data sets. In general, the first 40 m of the three profiles showed good relative sensitivity values. Between 40 and 60 m depth, the sensitivity was acceptable. The ERT-3 profile showed an area of low sensitivity values due to the presence of the resistive body. Prior to the inversion process, topographic data acquired every 50 m by an OCP surveyor team were added to the resistivity measurement files. An ABEM Terrameter LS 12 with 81 electrodes was used and made available by the Andaluz Institute of Geophysics (AIG) at Granada University, Spain. We adopted an edge gradient acquisition device which allowed for a better assembly between electrodes by combining symmetrical and asymmetrical configurations. Such a configuration provides an image of subsurface resistivities at high resolution and with less artifacts, e.g., distortions of the inversion resistivity model caused by the relatively high noise contamination of the data [24]. Three ERT profiles 700 m long with 5-m electrode spacing were collected ( Figure 2). A total of 6531 apparent resistivity measurements were acquired. Before the inversion of the data, their quality control was thoroughly controlled. No automatic filtering was applied but the data were edited, and all bad datum readings were manually eliminated. This resulted in slight reduction of the total corrected datum point for each profile. The data inversion was performed using the AIG's Res2dinv software (V.358-Geotomo Software). The program uses a smoothness-constrained Gauss-Newton method to obtain a 2D resistivity model [25]. The process is automatic once the inversion parameters are established and a user-defined initial model is not needed. The iterative inversion process was set at a maximum of 10 error-iterations using the L2-norm (root mean square), offering an error, between the apparent resistivities in the field and those generated by the calculated models, less than 5%. The inversion model contained 15 layers and 1692 resistivity blocks. It was refined to half the electrode spacing because of the sensitivity of the gradient array to near surface variation. The sensitivity matrix was calculated using the incomplete Gauss-Newton method that is recommended for large data sets. In general, the first 40 m of the three profiles showed good relative sensitivity values. Between 40 and 60 m depth, the sensitivity was acceptable. The ERT-3 profile showed an area of low sensitivity values due to the presence of the resistive body. Prior to the inversion process, topographic data acquired every 50 m by an OCP surveyor team were added to the resistivity measurement files.

MRS Data
Schematically, the physical principle of the MRS method is based on the fact that the protons that form the hydrogen nuclei of water molecules in Earth's magnetic field have positive magnetic moments that, at equilibrium, are aligned in the direction of this main field. The excitation of these protons by emitting a disturbing magnetic field at a specific frequency, termed a "Larmor frequency," modifies this equilibrium state and causes magnetic precession moments around Earth's magnetic field. After cutting off the exciter field, and during the return to their equilibrium state, a magnetic

MRS Data
Schematically, the physical principle of the MRS method is based on the fact that the protons that form the hydrogen nuclei of water molecules in Earth's magnetic field have positive magnetic moments that, at equilibrium, are aligned in the direction of this main field. The excitation of these protons by emitting a disturbing magnetic field at a specific frequency, termed a "Larmor frequency," modifies this equilibrium state and causes magnetic precession moments around Earth's magnetic field. After cutting off the exciter field, and during the return to their equilibrium state, a magnetic relaxation field is generated by the protons, thus constituting the MRS measured signal. The higher the signal amplitude, the larger the number of protons in resonance; thus, this provides information regarding the subsoil water content. The importance of the precession process in the excitation and relaxation at the field cut-off is also related to the average pore size of the aquifer formation. The relaxation time (decay time of the measured signal) is longer as the solicited protons are those of water less enclosed in the rock, and therefore of an aquifer with a high hydrodynamic potential guaranteeing a sufficient pumping flow rate [26][27][28][29].
First, we measured the ambient electromagnetic noise and assessed the quality of the MRS signal in the geological context of the Youssoufia mine. Then, an MRS survey was conducted using a loop of 100 × 100 m, which is the loop size of the MRS systems employed to realize the survey. We chose this dimension in order to maximize the depth of investigation of the study area. The transmitter/receiver loop was centered on well P6043 (Figure 2) to compare and calibrate the results of this method with real hydrogeological data. Numis Pro from IRIS Instruments was used from the University Cadi Ayyad of Marrakech. MRS data acquisition and inversion were, respectively, performed using the Prodiviner and Samovar software of the aforementioned company [30]. The measured parameter is the relaxation magnetic field created by excited protons when they return to their equilibrium state. The physical parameter that we sought to identify was the water content of the different subsoil layers. The result is a plot showing the hydrogeological model corresponding to the measured signal.

TDEM Data
Electromagnetic methods are based on the diffusion of an electromagnetic field (EM) in the subsoil to determine its electrical resistivity. For the TDEM soundings conducted at the Khouribga mining site, the EM field, termed the primary field, is created by the sudden cutting off of a flowing current in a transmitting coil placed on the ground. The secondary field, related to the induced current, is measured by a receiver coil after the current cut-off. The TDEM sounding curve provides the variation in the apparent resistivity of the subsoil according to the time, in other words according to the depth, which increases with time during the secondary field measurement [31,32]. Therefore, recorded curves can be converted to pseudo-depths [33,34]. The investigation depth depends on the formation resistivities, signal-noise ratio, as well as the measurement loop size. It is classically estimated at 1.5 to 2.5 times the size of the loop, but depends on the local conditions [35]. The Tem-Fast ProSystem, from the geological service of the OCP in Khouribga, was used. The surveys were conducted using a square emission loop of 20 m side pulled on the ground. A total of 84 TDEM soundings were completed at a sampling step of 12.5 m along a profile 1050 m in length ( Figure 3). Obtained data, from each survey, were inverted using the TemRes software based on the Temfast Prosystem technique [35].

FDEM Data
The FDEM method applied in the present study was carried out using the Geonics EM-31 MK-2 system. This technique helps record the electrical conductivity whose variations reflect the facies heterogeneities in the near subsoil. The measurements were performed without any contact with the ground; this makes the field survey very easy to implement and provides a higher acquisition performance than the TDEM technique. The EM31-MK-2 system contains a transmitter and receiver coils fixed at the ends of a 3-m-long support bar and connected to the control box in the middle of the device. The transmitter generates a primary magnetic field at a given frequency. A secondary field is generated and recorded by the receiver when the primary field encounters conducting areas in the ground. The ratio of the vertical component of the secondary field in the quadrature with the primary field is proportional to the apparent conductivity of the investigated soil. The investigation depth of the EM31 method is theoretically equal to 1.5 to 2 times the distance between the transmitting and receiver coils, for a medium-conducting environment [36][37][38].
The EM31 method was used at the Sidi Chennane deposit (Figure 3), in the same area where the TDEM surveys were conducted. Measurements were made along profiles spaced every 10 m. These profiles cover 33 juxtaposed parcels of 100 × 100 m equaling a total length of 36.3 km and 207,746 measured values. A global positioning system (GPS) was connected and synchronized with the EM31 acquisition console, allowing for measurement geolocalization.

FDEM Data
The FDEM method applied in the present study was carried out using the Geonics EM-31 MK-2 system. This technique helps record the electrical conductivity whose variations reflect the facies heterogeneities in the near subsoil. The measurements were performed without any contact with the ground; this makes the field survey very easy to implement and provides a higher acquisition performance than the TDEM technique. The EM31-MK-2 system contains a transmitter and receiver coils fixed at the ends of a 3-m-long support bar and connected to the control box in the middle of the device. The transmitter generates a primary magnetic field at a given frequency. A secondary field is generated and recorded by the receiver when the primary field encounters conducting areas in the ground. The ratio of the vertical component of the secondary field in the quadrature with the primary field is proportional to the apparent conductivity of the investigated soil. The investigation depth of the EM31 method is theoretically equal to 1.5 to 2 times the distance between the transmitting and receiver coils, for a medium-conducting environment [36][37][38].
The EM31 method was used at the Sidi Chennane deposit (Figure 3), in the same area where the TDEM surveys were conducted. Measurements were made along profiles spaced every 10 m. These profiles cover 33 juxtaposed parcels of 100 × 100 m equaling a total length of 36.3 km and 207,746 measured values. A global positioning system (GPS) was connected and synchronized with the EM31 acquisition console, allowing for measurement geolocalization.

Youssoufia Mining Site
The study area belongs to the phosphate deposit of Panel 6, approximately 26 km east of Youssoufia on the main road No 9 linking Marrakech to El Jadida. Figure 4 shows the results of ERT profiles P1, P2, and P3. This figure also shows the synthetic lithological columns of the recognition wells of the phosphatic deposits created by OCP. This allowed us to constrain the models of the calculated resistivities and assign the measured values to the different lithological units of the phosphatic series. We also show on the same figure the curve of the piezometric level measured in the OCP boreholes during the geophysical field surveys. The study area belongs to the phosphate deposit of Panel 6, approximately 26 km east of Youssoufia on the main road No 9 linking Marrakech to El Jadida. Figure 4 shows the results of ERT profiles P1, P2, and P3. This figure also shows the synthetic lithological columns of the recognition wells of the phosphatic deposits created by OCP. This allowed us to constrain the models of the calculated resistivities and assign the measured values to the different lithological units of the phosphatic series. We also show on the same figure the curve of the piezometric level measured in the OCP boreholes during the geophysical field surveys. In the southwest part of the three ERT profiles, the resistivity values are less than 20 ohm‧m in the sub-surface. These low values can be explained by the more or less wet marl and clay facies that we observed in this area during the geophysical surveys. The remainders of the profiles at this depth (i.e., about 10 to 20 m) are represented by higher resistivity values greater than 140 ohm‧m. These values are attributed to the phosphatic limestone and uncemented phosphatic dry layers of the Upper Maastrichtian and Thanetian, between 0 and 25 m depth, in the OCP's lithological wells, but also in the surrounding mining trench. Underneath, clearly conductive areas with resistivities less than 60 In the southwest part of the three ERT profiles, the resistivity values are less than 20 ohm·m in the sub-surface. These low values can be explained by the more or less wet marl and clay facies that we observed in this area during the geophysical surveys. The remainders of the profiles at this depth (i.e., about 10 to 20 m) are represented by higher resistivity values greater than 140 ohm·m. These values are attributed to the phosphatic limestone and uncemented phosphatic dry layers of the Upper Maastrichtian and Thanetian, between 0 and 25 m depth, in the OCP's lithological wells, but also in the surrounding mining trench. Underneath, clearly conductive areas with resistivities less than 60 ohm·m occur. They correspond to the damped or high-water content of the uncemented phosphatic layers of the Lower Maastrichtian and to Maastrichtian impermeable clays (2 m) that constitute a reference level for the extraction of the Youssoufia open-pit mining [10,17]. Conductive layers still occur beneath the Maastrichtian formations until the bottom of the ERT sections. They correspond to the Senonian marls described in the stratigraphic columns of the OCP wells since the depth of 50 m. The most conductive horizon of all the ERT sections is situated between 25 and 35 m depth. It corresponds to the uncemented phosphatic layer of the Lower Maastrichtian that overlies the aforementioned impermeable clays. At some locations, it is difficult to make delimitation according to the depth of the top and the bottom of the real aquifer horizons, corresponding to Maastrichtian uncemented phosphatic layers and marl and/or limestone intercalary, which are probably very wet. Both formations provide comparable electrical signatures, as reflected in the ERT sections by the thick blue section. Therefore, the sensitivity according to the vertical of the ERT method is more or less limited. However, laterally, the ERT sections show moderately resistive zones within the same conductive stratigraphic horizon. This lateral variation in the resistivities can be attributed to a decrease in the water content within the same conductive horizon or to local facies change. These facies variations within the same layer are actually described along the operating trenches.
The piezometric level measured in the P1821, P6123, and P6043 boreholes during the geophysical surveys exactly coincide with the interpretation provided to the ERT sections. We interpolated between the OCP boreholes and extrapolated outside of them, assigning aquifer resistivity values less than 60 ohm·m. The perfect match between the depth of the piezometric level measured in the OCP boreholes and the roof of the conductive horizon was shown by the ERT data inversion (Figure 4). Thus, this method can be used to identify and locate the aquifers in the phosphatic series at the Youssoufia mining site.
During the present study reliable MRS data were recorded after some unsuccessful attempts to perform the measurements with low stacking values, even if the survey area was located far from any electromagnetic source of noise. To obtain data of good signal-to-noise ratio (average of 5.5), it was necessary to take the measurements using 250 stacks. The emitted and the received signals had very close frequencies (difference less than 1.0 Hz). A smooth inversion of the MRS data was performed using Samovar software with a filtering window of 197.8 ms, a time constant of 15.00 ms, an average signal to noise ratio of 2.66, and a fitting error for the free induction decay 1 (FID1) of 9.05%. Figure 5 shows MRS survey data and results including the signal relaxation curves versus time for the various pulse moments injected (Figure 5a). The hydrogeological model resulting from the inversion helps determine the depth of the aquifers and the free water content that estimates the porosity of the rock in a saturated environment. One can clearly see that the amplitude of the measured signal changed according to the excitation pulse and thus the depth. This shows the presence of subsoil water and the variation in water content from one level to another (Figure 5b). We also note that the MRS signal was clearly distinguished from the electromagnetic noise, which was fortunately very low at the studied site. Generally, two aquifer levels can be distinguished (Figure 5c): The first is between 20 and 30 m, and the second between 32 and 45 m. The water content, which reflects the porosity of the layers below the piezometric level, is, respectively, 3% and 6.7%. These values are quite comparable to those provided by pumping tests conducted on the same aquifer outside the study area [22]. The average aquifer thickness derived from the MRS calculated model was approximately 23 m. The depth of the piezometric level measured in well n°6043 at the center of the MRS loop was 21.23 m. The MRS signal's decay time was relatively low (100-150 ms) which meant that, in general, the aquifer was moderately permeable [39][40][41]. The recognition well lithological sections in the phosphatic deposits showed that the aquifer delimited by the MRS data inversion corresponded to the two uncemented phosphatic layers of the Lower Maastrichtian in the study area (Figure 5c). Given that the studied area was approximately 19,242,500 m 2 , we deduced a total water volume of approximately 22,128,875 +/-2,251,373 m 3 . The uncertainty of the volume estimate was about 10% due to the error of average porosity calculated from MRS inversion and the quantification of the unsaturated part of the aquifer. This volume was quite comparable to the that obtained via the direct calculation using the OCP well piezometric data (22,829,600 +/-2,282,960 m 3 ). It is with this purpose that the study was conducted using MRS to better understand the water reserves and hydrodynamic The average aquifer thickness derived from the MRS calculated model was approximately 23 m. The depth of the piezometric level measured in well n • 6043 at the center of the MRS loop was 21.23 m. The MRS signal's decay time was relatively low (100-150 ms) which meant that, in general, the aquifer was moderately permeable [39][40][41]. The recognition well lithological sections in the phosphatic deposits showed that the aquifer delimited by the MRS data inversion corresponded to the two uncemented phosphatic layers of the Lower Maastrichtian in the study area (Figure 5c). Given that the studied area was approximately 19,242,500 m 2 , we deduced a total water volume of approximately 22,128,875 +/-2,251,373 m 3 . The uncertainty of the volume estimate was about 10% due to the error of average porosity calculated from MRS inversion and the quantification of the unsaturated part of the aquifer. This volume was quite comparable to the that obtained via the direct calculation using the OCP well piezometric data (22,829,600 +/-2,282,960 m 3 ). It is with this purpose that the study was conducted using MRS to better understand the water reserves and hydrodynamic parameters of the Maastrichtian aquifer. These results strongly support the use of this method in other areas of the Youssoufia open-pit mines.

Khouribga Mining Site
To facilitate interpretation of the results, we used the same color palette for both resistivity (TDEM method) and resistivity calculated from conductivity (EM31 method). The TDEM sounding data are presented as a pseudo-section obtained by interpolation of the experimental apparent resistivities from one sounding to another using the minimum curvature method [42] (Figure 6a). The same method was used to interpolate the 1D inversion of the 84 performed soundings in order to build the resistivity model illustrated in Figure 6b. This model shows a clearly conductive horizon whose top varies approximately in altitude from 495 to 510 m, and in depth from 24 to 47 m. This conductive horizon is continuous and extends for more than 1 km along the TDEM executed profile. Beneath the 70 m depth, the top of moderately resistant layers is locally observed. They correspond to the compact gray marl with gypsum of the Senonian, which can be seen in the OCP's recognition wells at this depth [43]. Above the conductive horizon, moderately resistant and resistant soils consisting of phosphatic layers with very low water content with intercalated limestone and flint layers occur. At the near-surface, the layers are even more resistant because they are dry and rich in limestone encrustation. Other resistant bodies also cut the phosphatic series. They are attributed to the disturbed areas of the phosphatic series, known as "derangements" by the OCP mining engineers [43][44][45][46]. These structures are often accompanied by a localized silicification event. The research of Kchikach et al. 2002Kchikach et al. , 2006Kchikach et al. , and 2012 shows that the disturbed areas of the phosphatic series are characterized by a strongly resistive electrical signature.

Khouribga Mining Site
To facilitate interpretation of the results, we used the same color palette for both resistivity (TDEM method) and resistivity calculated from conductivity (EM31 method). The TDEM sounding data are presented as a pseudo-section obtained by interpolation of the experimental apparent resistivities from one sounding to another using the minimum curvature method [42] (Figure 6a). The same method was used to interpolate the 1D inversion of the 84 performed soundings in order to build the resistivity model illustrated in Figure 6b. This model shows a clearly conductive horizon whose top varies approximately in altitude from 495 to 510 m, and in depth from 24 to 47 m. This conductive horizon is continuous and extends for more than 1 km along the TDEM executed profile. Beneath the 70 m depth, the top of moderately resistant layers is locally observed. They correspond to the compact gray marl with gypsum of the Senonian, which can be seen in the OCP's recognition wells at this depth [43]. Above the conductive horizon, moderately resistant and resistant soils consisting of phosphatic layers with very low water content with intercalated limestone and flint layers occur. At the near-surface, the layers are even more resistant because they are dry and rich in limestone encrustation. Other resistant bodies also cut the phosphatic series. They are attributed to the disturbed areas of the phosphatic series, known as "derangements" by the OCP mining engineers [43][44][45][46]. These structures are often accompanied by a localized silicification event. The research of Kchikach et al. 2002Kchikach et al. , 2006Kchikach et al. , and 2012 shows that the disturbed areas of the phosphatic series are characterized by a strongly resistive electrical signature. The apparent resistivity map obtained after smoothing the EM31 data, and at regular sampling steps of 1 × 1 m, shows a series of conductive anomalies in the studied area. These anomalies, with conductivity values generally greater than 20 mS/m (Figure 7), consist of altered and fractured zones of limestone encrustation. However, the areas where these encrustations are thicker and less The apparent resistivity map obtained after smoothing the EM31 data, and at regular sampling steps of 1 × 1 m, shows a series of conductive anomalies in the studied area. These anomalies, with conductivity values generally greater than 20 mS/m (Figure 7), consist of altered and fractured zones of limestone encrustation. However, the areas where these encrustations are thicker and less fractured result in conductivities less than 15 mS/m. Conductive anomalies are superimposed with fractured and altered corridors that can be easily delimited in the field. Actually, during the geophysical surveys, these corridors were distinguished in the field by their wet aspect. They constitute potential zones for surface water infiltration towards the identified aquifer in the same areas as shown using the TDEM method. Despite its limited investigation depth, the EM31 technique is a very good tool for mapping groundwater aquifer recharge zones in the phosphatic series. The most productive water wells and boreholes should be, in our estimation, installed in line with these highly conductive anomalies. fractured result in conductivities less than 15 mS/m. Conductive anomalies are superimposed with fractured and altered corridors that can be easily delimited in the field. Actually, during the geophysical surveys, these corridors were distinguished in the field by their wet aspect. They constitute potential zones for surface water infiltration towards the identified aquifer in the same areas as shown using the TDEM method. Despite its limited investigation depth, the EM31 technique is a very good tool for mapping groundwater aquifer recharge zones in the phosphatic series. The most productive water wells and boreholes should be, in our estimation, installed in line with these highly conductive anomalies. Finally, a summary illustration synthesizing the main outputs of the present study was established ( Figure 8). The latter clearly shows that in both Youssoufia and Khouribga areas, the Lower Maastrichtian constitute a high potential water-bearing aquifer. The figure clearly depicts the convergence between the ERT and the MRS findings in Youssoufia. In fact, the conductive character of the Lower Maastrichtian is explained by the water content evidenced by the MRS survey. In the Khouribga site, TDEM data revealed that the Lower Maastrichtian is also clearly a conductor. Since this age (Lower Maastrichtian) is made with the same facies as Youssoufia (uncemented phosphate), we can assume that it may constitute the main reservoir of underground water.  Finally, a summary illustration synthesizing the main outputs of the present study was established ( Figure 8). The latter clearly shows that in both Youssoufia and Khouribga areas, the Lower Maastrichtian constitute a high potential water-bearing aquifer. The figure clearly depicts the convergence between the ERT and the MRS findings in Youssoufia. In fact, the conductive character of the Lower Maastrichtian is explained by the water content evidenced by the MRS survey. In the Khouribga site, TDEM data revealed that the Lower Maastrichtian is also clearly a conductor. Since this age (Lower Maastrichtian) is made with the same facies as Youssoufia (uncemented phosphate), we can assume that it may constitute the main reservoir of underground water.

Conclusions
Minerals 2020, 10, x FOR PEER REVIEW 12 of 16 fractured result in conductivities less than 15 mS/m. Conductive anomalies are superimposed with fractured and altered corridors that can be easily delimited in the field. Actually, during the geophysical surveys, these corridors were distinguished in the field by their wet aspect. They constitute potential zones for surface water infiltration towards the identified aquifer in the same areas as shown using the TDEM method. Despite its limited investigation depth, the EM31 technique is a very good tool for mapping groundwater aquifer recharge zones in the phosphatic series. The most productive water wells and boreholes should be, in our estimation, installed in line with these highly conductive anomalies. Finally, a summary illustration synthesizing the main outputs of the present study was established ( Figure 8). The latter clearly shows that in both Youssoufia and Khouribga areas, the Lower Maastrichtian constitute a high potential water-bearing aquifer. The figure clearly depicts the convergence between the ERT and the MRS findings in Youssoufia. In fact, the conductive character of the Lower Maastrichtian is explained by the water content evidenced by the MRS survey. In the Khouribga site, TDEM data revealed that the Lower Maastrichtian is also clearly a conductor. Since this age (Lower Maastrichtian) is made with the same facies as Youssoufia (uncemented phosphate), we can assume that it may constitute the main reservoir of underground water.

Conclusions
This study clearly shows the efficiency of the MRS method for prospecting groundwater resources and evaluating their importance in the geological context of Youssoufia open-pit mining. The hydrogeological model resulting from the MRS data inversion corroborates with the real hydrogeological data of the study area. The highlighted aquifer corresponds to two layers of high-water content Maastrichtian uncemented phosphate, as shown by the extraction works in the surrounding areas. The water volumes estimated via the piezometric data analysis and MRS method are quite comparable. Despite the average signal amplitude, the exceptional low electromagnetic noise conditions, in the geological context of the OCP open-pit mines, encourage the use of this method to survey and characterize the groundwater aquifers in the phosphatic series.
ERT profiles allowed us to delimit in detail the conductive horizons attributed to the groundwater aquifers. They clearly show that aquifers in the study area are situated below 21 m depth. Compared to the MRS method, the ERT method is less sensitive; the aquifer horizons corresponding to the uncemented phosphatic layers and intercalary of the wet marls and limestones show comparable electrical signatures, thus limiting the precise location of the top and the bottom of the supposedly saturated aquifer. However, this method is sensitive to lateral resistivity variations along the same aquifer horizon. These variations reflect the changes in water content from one zone to another along the horizon.
TDEM surveys at the Khouribga mining site highlighted a conductive horizon whose depth varies from 24 to 47 m. The interpolation between the models of the 84 surveys conducted shows that this horizon is continuous throughout the prospected zone. The integration of the lithological data sections of the surrounding boreholes enables one to attribute this conductance to the Lower Maastrichtian phosphatic layers. These layers overlie compact and impermeable Senonian marls that prevent deeper groundwater infiltration. Therefore, uncemented phosphatic layers of the Lower Maastrichtian constitute the most important aquifer in the geological context of the Khouribga open-pit mines. Conducting larger geophysical surveys, combining the TDEM and MRS surveys, which provide satisfying results in the similar hydrogeological context of the Youssoufia deposits, would enable a better understanding of the Lower Maastrichtian aquifer and evaluation of its hydrodynamic parameters.
Despite its limited investigation depth, the FDEM method can be used as a tool for mapping and delimiting the aquifer potential recharge zones in the phosphate series. The superposition of the conductivity map obtained from the FDEM data on the TDEM results, and on the spatial distribution of the aquifer hydrodynamic parameters deduced from the MRS data, could guide borehole or well installation in the most productive areas for groundwater extraction. The location of these areas would allow OCP mining engineers to implement all the necessary preparations and arrangements to extract phosphates in these flooded areas.

Conflicts of Interest:
The authors declare no conflict of interest.