**Geophysical Prospecting for Groundwater Resources in Phosphate Deposits (Morocco)**

**Fatim-Zahra Ihbach <sup>1</sup> , Azzouz Kchikach 1,2,\*, Mohammed Ja**ff**al 1,2 , Driss El Azzab <sup>3</sup> , Oussama Khadiri Yazami <sup>4</sup> , Es-Said Jourani <sup>4</sup> , José Antonio Peña Ruano <sup>5</sup> , Oier Ardanaz Olaiz <sup>5</sup> and Luis Vizcaíno Dávila <sup>5</sup>**


Received: 27 May 2020; Accepted: 10 August 2020; Published: 24 September 2020

**Abstract:** 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 flooded 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 Youssoufia 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.

**Keywords:** phosphate mines; geophysical exploration; hydrogeology; phosphate extraction; gantour Basin; Ouled Abdoun Basin; Morocco

#### **1. 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).

**Figure 1.** Location and geology sitting of study areas. (**a**) Geological map of the Gantour Basin (modified after [8]); (**b**) Geological map of the Ouled Abdoun Basin (after [9]); (**c**) Lithological columns of the phosphatic series in Gantour Basin showing the main groundwater aquifers; (**d**) Lithological columns of the phosphatic series in Ouled Abdoun Basin showing the main groundwater aquifers; 1: Paleozoic; 2: Trias; 3: Jurassic; 4: Cretaceous except Maastrichtian; 5: Maastrichtian–Eocene; 6: Uncemented phosphates, 7: Marl; 8: Phosphatic marl; 9: Limestone; 10: Phosphatic limestone; 11: Marly limestone; 12: Silex; 13: Clay; 14: Marly-siliceous limestone; 15: Open-pit mine; 16: Aquifers. (N.B. Lambert kilometric coordinates projection for the north of Morocco).

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.

#### **2. 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–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<sup>2</sup> 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–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.

#### **3. 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.

#### *3.1. 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.

**Figure 2.** Location of geophysical surveys carried out in Youssoufia area and photograph showing a borehole intersected by the electrical resistivity tomography (ERT)-2 profile.

#### *3.2. 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–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.

#### *3.3. 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].

**Figure 3.** Location of geophysical surveys conducted in Khouribga area.

#### *3.4. 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–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.

#### **4. Results and Discussion**

#### *4.1. 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.

**Figure 4.** Model resulting from ERT data inversion calibrated to lithological boreholes columns. The dashed white line corresponds to the piezometric level measured on site. The dashed black line indicates the limit between the Maastrichtian and the Senonian formations. 1: Uncemented phosphate; 2: Limestone; 3: Clay; 4: Upper Maastrichtian–Thanetian (dry uncemented phosphate); 5: Lower Maastrichtian: (wet uncemented phosphate); 6: Senonian: (marl).

‧

‧

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].

**Figure 5.** (**a**) Magnetic resonance sounding (MRS) signal relaxation curves versus time for the various pulse moments injected; (**b**) amplitude of measured MRS signal and noise level; (**c**) model resulting from MRS data inversion and its correlation with hydrogeological column; 1: Marl; 2: Phosphatic marl; 3: Uncemented phosphate; 4: Phosphatic limestones.

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–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<sup>2</sup> , we deduced a total water volume of approximately 22,128,875 +/- 2,251,373 m<sup>3</sup> . 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<sup>3</sup> ). 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.

#### *4.2. 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–46]. These structures are often accompanied by a localized silicification event. The research of Kchikach et al. 2002, 2006, and 2012 [16,47,48] shows that the disturbed areas of the phosphatic series are characterized by a strongly resistive electrical signature.

**Figure 6.** (**a**) Time-domain electromagnetic (TDEM) pseudo-section of experimental resistivity data from TDEM sounding survey, (**b**) interpolated 1D inversion resistivity model.

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.

**Figure 7.** Apparent resistivity map obtained from EM31 data.

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 LowerMaastrichtian 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.

**Figure 8.** Summary illustration synthesizing the main outputs of the present study. 1: UpperMaastrichtian– Thanetian (Dry uncemented phosphate); 2: Upper Maastrichtian: (Dry uncemented phosphate); 3: Lower Maastrichtian: (Wet uncemented phosphate); 4: Senonian: (Compact gray marl); 5: Aquifer water content derived from MRS survey; 6: Maximal investigated depth for each geophysical method.

#### **5. 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.

**Author Contributions:** Conceptualization: A.K. and F.-Z.I.; methodology: A.K., M.J., and D.E.A.; software: F.-Z.I., M.J., J.A.P.R., O.A.O., and L.V.D.; validation: A.K., F.-Z.I., and M.J.; investigation: F.-Z.I., A.K., M.J., D.E.A., J.A.P.R., O.A.O., and L.V.D.; resources: E.-S.J. and O.K.Y.; writing—original draft preparation: F.-Z.I. under the supervision of A.K.; writing—review and editing: A.K., M.J., D.E.A., O.K.Y., E.-S.J., J.A.P.R., O.A.O., and L.V.D.; project administration: A.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by OCP-SA and CNRST Morocco project APHOS, GEO-KCH-01/2017.

**Acknowledgments:** This research was carried out as part of the project ID GEO-KCH 01/2017, funded by OCP. Particular thanks are addressed to the managers of the Youssoufia and Khouribga mining sites for their help and logistical support.

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

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Parallel Simulation of Audio- and Radio-Magnetotelluric Data**

#### **Nikolay Yavich 1,2,\* , Mikhail Malovichko 1,2 and Arseny Shlykov <sup>3</sup>**


Received: 22 November 2019; Accepted: 27 December 2019; Published: 31 December 2019 -

**Abstract:** This paper presents a novel numerical method for simulation controlled-source audio-magnetotellurics (CSAMT) and radio-magnetotellurics (CSRMT) data. These methods are widely used in mineral exploration. Interpretation of the CSAMT and CSRMT data collected over an area with the complex geology requires application of effective methods of numerical modeling capable to represent the geoelectrical model of a deposit well. In this paper, we considered an approach to 3D electromagnetic (EM) modeling based on new types of preconditioned iterative solvers for finite-difference (FD) EM simulation. The first preconditioner used fast direct inversion of the layered Earth FD matrix (Green's function preconditioner). The other combined the first with a contraction operator transformation. To illustrate the effectiveness of the developed numerical modeling methods, a 3D resistivity model of Aleksandrovka study area in Kaluga Region, Russia, was prepared based on drilling data, AMT, and a detailed CSRMT survey. We conducted parallel EM simulation of the full CSRMT survey. Our results indicated that the developed methods can be effectively used for modeling EM responses over a realistic complex geoelectrical model for a controlled source EM survey with hundreds of receiver stations. The contraction-operator preconditioner outperformed the Green's function preconditioner by factor of 7–10, both with respect to run-time and iteration count, and even more at higher frequencies.

**Keywords:** geophysical electromagnetic modelling; CSAMT; CSRMT; CSEM

#### **1. Introduction**

Controlled-source audio-magnetotellurics (CSAMT) is a popular method for mineral exploration [1–5]. In principle, a non-uniform plane wave at a study area is generated by two orthogonal transmitter lines located at a sufficient distance (in the far-field zone). A typical operation frequency range is 0.1–10 kHz, which is roughly the high end of the quasi-stationary regime. Transfer functions between electric and magnetic field components, *E*x, *E*y, *H*x, *H*y, and *H*z are computed from measured data and then inverted in order to reconstruct the subsurface electric conductivity. CSAMT is known to provide better data quality in the AMT dead bands (0.1–5 Hz and 700–3000 Hz) than conventional AMT surveys and it requires relatively inexpensive transmitter operations [2,5,6].

Closely connected to CSAMT is the controlled-source radio-magnetotellurics (CSRMT). This method is gaining popularity in the near-surface geophysics [7–11]. It operates at much higher frequencies (typically, up to 250 kHz–1 MHz) and thus provides excellent spatial resolution. The biggest difference with CSAMT is the displacement currents in the air [10,12] cannot be neglected to accurately compute the electric and magnetic fields, which complicates the numerical simulation dramatically. The authors of [13] demonstrated, however, that the full impedance and quasi-stationary impedances are almost identical in the 1D Earth. It suggests that the CSRMT measurements can be interpreted using an inversion software with the quasi-stationary forward simulation engine (note, though, [14,15]).

Numerical forward modeling plays a paramount role in both electromagnetic (EM) survey design and data interpretation. This demand has been driving the investigation of efficient modeling methods for nearly 60 years. Today, 1D and 2D frequency-domain modeling is relatively established, while research on 3D modeling algorithms is a hot topic.

The common problem in 3D modeling of CSAMT and CSRMT data is the solution of a large system of linear equations. This problem is present in all the discretization methods: Finite difference (FD), finite element (FE), integral equation (IE), and others; although in FD and FE modeling the system matrix is sparse, while in IE modeling the system matrix is dense. Consequently, design and implementation of an economical, scalable, and robust solver for this system is the key to successful 3D modeling and inversion.

At the frequencies from 0 to 1 MHz, the EM fields are mainly diffusive. Thus, the errors introduced by resistivity model resampling are local. Consequently, adaptation of the grid to a particular interface is a minor concern. We thus consider the second-order edge-based finite-difference discretization.

Computer simulation of both CSAMT and CSRMT data poses great computational challenges. The land surveys are frequently characterized by an exceptionally large resistivity contrast of the rocks, e.g., 1000 Ωm for granite and 0.01 Ωm for nickel ore [16], leading to very fine numerical grids needed to resolve the shortest skin depth. At the same time, the computational domain should be extended to include attenuation dictated by the longest skin depth, especially at frequencies above ~1 kHz, at which considerable amount of the electromagnetic field travels in the air and resistive rocks with little attenuation. Furthermore, in CSAMT/CSRMT setups the transmitters are grounded, making their numerical treatment more difficult comparing to the conventional audio and radio magnetotellurics (AMT/RMT), respectively. For these reasons, a CSAMT/CSRMT simulation easily results in a discrete problem with millions of discrete unknowns. Beyond being very large, the system matrix is extremely ill-conditioned. Thus, both preconditioned iterative solvers and sparse direct solvers are commonly considered. Although the use of unstructured tetrahedral grids [17] or hexahedral grids with hanging nodes [18] reduces the size of the discrete problem, an efficient solver still will be required.

Direct solvers are attractive for their universality and robustness [6,19], though they require very long initialization and large memory allocations. Thus, iterative solvers are more typically used in practical modeling and inversion. The commonly applied preconditioners are incomplete factorization completed by divergence correction [20] and multigrid [21,22].

In this paper, we assess performance of the two preconditioned iterative solvers introduced by the authors in [23]. They have been successfully tested on MT and marine controlled-source electromagnetic (CSEM) [23] and land CSEM [24] setups. Here, we applied the two preconditioners to quasi-stationary simulation in the CSRMT frequency range. Both serial and parallel modeling performance were evaluated. The solvers were applied to the secondary field formulation to avoid excessive gridding near the source.

The paper is organized as follows. We review forward modelling and briefly discuss implementation in Section 2. In Section 3, we verify modeling accuracy of our code versus an integral equation code [25]. Next, we present results of numerical simulation of a study area in Moscow syneclise. We compare performance of the two preconditioners in a single-threaded mode. Then we present a complete simulation over a realistic geologic model conducted on a cluster. Final remarks are given in Section 4.

#### **2. Three-Dimensional Simulation of CSAMT**/**CSAMT Data**

#### *2.1. E*ffi*cient Finite-Di*ff*erence Simulation*

In CSAMT/CSRMT, the five components of the electromagnetic filed are measured (*E*x, *E*y, *H*x, *H*y, *H*z) and, after preprocessing, converted to the components of the impedance and tipper. The off-diagonal components of the impedance are employed most commonly, that is [2],

$$\begin{aligned} Z\_{xy} &= \frac{E\_{x1}H\_{x2} - E\_{x2}H\_{x1}}{H\_{x1}H\_{y2} - H\_{x2}H\_{y1}} \\ Z\_{yx} &= \frac{E\_{y1}H\_{y2} - E\_{y2}H\_{y1}}{H\_{x1}H\_{y2} - H\_{x2}H\_{y1}} \end{aligned} \tag{1}$$

where subscripts 1 and 2 refer to one of the two orthogonal transmitter lines. Thus, computation of the components of the impedance or tipper requires simulation of the electric and magnetic fields *E* and *H* at several locations and frequencies.

We consider electromagnetic modeling in 3D heterogeneous isotropic media in the frequency domain. Assuming time dependence of *e* −*i*ω*t* , the electric field, *E*(*x*, *y*, *z*), satisfies the following system of partial differential equations,

$$
\text{i } \text{curl } \text{curl } \mathbf{E} - i\omega\mu\_0 \sigma \mathbf{E} = i\omega\mu\_0 \mathbf{F} \tag{2}
$$

where ω is the source angular frequency, µ<sup>0</sup> is the magnetic permeability of the free space, σ(*x*, *y*, *z*) is the conductivity model, and *F*(*x*, *y*, *z*) is the source current density. This is a linear system of three scalar equations with respect to three scalar entries of *E*(*x*, *y*, *z*). We solve this partial-differential equation system in a bounded domain, completed by zero Dirichlet boundary conditions. The magnetic field *H* is then computed via Faraday's law,

$$\mathbf{H} = \frac{1}{i\omega\mu\_0} \operatorname{curl} \mathbf{E}.\tag{3}$$

Since geological models predominately vary in the vertical direction, we can always split the conductivity model into background and anomalous parts,

$$
\sigma(\mathbf{x}, \mathbf{y}, \mathbf{z}) = \sigma\_\vartheta(\mathbf{z}) + \sigma\_\vartheta(\mathbf{x}, \mathbf{y}, \mathbf{z}). \tag{4}
$$

We further assume that the following double inequality holds,

$$
\rho\_\vartheta \sigma\_\vartheta(z) \le \sigma(x, y, z) \le \beta \sigma\_\vartheta(z), 0 < a \le 1 \le \beta,\tag{5}
$$

which controls the contrast of anomalous conductivity with respect to the background.

Given a non-uniform computational grid with *N<sup>x</sup>* × *N<sup>y</sup>* × *N<sup>z</sup>* cells, we apply the conventional edge-based second-order FD discretization. The FD method approximates the unknown electric field *E* = *Ex*, *Ey*, *Ez* with a finite set of discrete values *Ei*<sup>+</sup> <sup>1</sup> 2 *j k*, *Ei j*<sup>+</sup> <sup>1</sup> 2 *k* , *Ei j k*<sup>+</sup> <sup>1</sup> 2 [26]. Each discrete value is attached to the respective edge of the FD grid, Figure 1.

 =

 =

 <sup>0</sup> (, , ) (, , )

= (, , ) {+

− 0 = 0

Faraday's law,

() +

(, , ).

(), 0 < ≤ 1 ≤ ,

× ×

1 2 , + 1 2 , + 1 2 }

1 0

=

(, , ) =

() ≤ (, , ) ≤

12−21 12−21

12−21 12−21

− (, , )

(, , )

**Figure 1.** Discrete electric fields attached to the edges of the grid.

The actual discrete equations can be found for example in [27]. Let us form the unknown vector *e* of the discrete electric fields *Ei*<sup>+</sup> <sup>1</sup> 2 *j k* , *Ei j*<sup>+</sup> <sup>1</sup> 2 *k* , *Ei j k*<sup>+</sup> <sup>1</sup> 2 introduced above. Now we can write the FD discretization of (2) in a matrix form:

$$\mathbf{A}\,\mathbf{e} = i\omega\mu\_0 \mathbf{f}\_\prime \tag{6}$$

where *A* is a complex symmetric system matrix with at most 13 nonzero entries per row, corresponding to the total conductivity distribution. We denote the size of this system as *n*, *n* ≈ 3*NxNyNz*.

Let *A<sup>b</sup>* be the FD system matrix, corresponding to the background conductivity model. Importantly, this matrix can be implicitly factorized using a fast direct algorithm, and the action of the inverse matrix can be efficiently computed [23,28]. As a result, it can be used as a preconditioner to (6):

$$\mathbf{A}\_b^{-1} \mathbf{A} \ e = i\omega \mu\_0 \mathbf{A}\_b^{-1} \mathbf{f}. \tag{7}$$

We will refer *A* −1 *b* as the Green's function (GF) preconditioner. The complexity of applying the GF preconditioner is *O n* 4/3 , and auxiliary memory required is near 3*n* only. Applying the analysis presented in [23] the condition number of the GF preconditioned system can be estimated as follows,

$$\alpha \text{cond} \left( \mathbf{A}\_b^{-1} \mathbf{A} \right) \le \beta / \alpha. \tag{8}$$

This result implies that convergence of an iterative solver applied to (7) has minor or no dependence on the grid size and cell aspect ratio, as well as the frequency, while it degrades on models with high-contrast bodies.

To minimize the impact of high-contrast bodies, another preconditioner could be constructed. Denote as **Σ** and **Σ***<sup>b</sup>* diagonal matrices, corresponding to full and background conductivity, respectively. Let us define the modified FD Green's operator and diagonal matrices according to the following formulae:

$$\mathcal{G}\_b^M = 2i\omega\mu\_0 \boldsymbol{\Sigma}\_b^{\frac{1}{2}} \boldsymbol{A}\_b^{-1} \boldsymbol{\Sigma}\_b^{\frac{1}{2}} + \mathbf{I},\tag{9}$$

$$\mathbf{K}\_1 = \frac{1}{2} (\boldsymbol{\Sigma} + \boldsymbol{\Sigma}\_b) \boldsymbol{\Sigma}\_b^{-\frac{1}{2}}, \ \mathbf{K}\_2 = \frac{1}{2} (\boldsymbol{\Sigma} - \boldsymbol{\Sigma}\_b) \boldsymbol{\Sigma}\_b^{-\frac{1}{2}}.\tag{10}$$

Using this operator, Equation (7) can be written in an equivalent form as follows (see Appendix A):

$$\hat{\mathbf{e}} = \mathbf{G}\_b^M \, \mathbf{K}\_2 \mathbf{K}\_1^{-1} \hat{\mathbf{e}} + i\omega\mu\_0 \boldsymbol{\Sigma}\_b^{\frac{1}{2}} \mathbf{A}\_b^{-1} \mathbf{f},\tag{11}$$

where ˆ*e* = *K*1*e*. By introducing a new operator, *C* = G *M b K*2*K* −1 1 , we rewrite Equation (11) as follows:

$$\mathbf{e}(\mathbf{I}-\mathbf{C})\hat{\mathbf{e}} = i\omega\mu\_0 \Sigma\_b^{\frac{1}{2}} \mathbf{A}\_b^{-1} \mathbf{f}.\tag{12}$$

We will refer this system a contraction operator (CO) preconditioned system, since it was shown in [23] that k*C*k < 1. Moreover, it can be proved that this preconditioned system has a smaller or equal spectral condition number than that of the GF preconditioner, implying faster or equal convergence of the iterative solvers,

$$
\max(\mathbf{I} - \mathbf{C}) \le \max\{\mathbf{1}/\alpha, \boldsymbol{\beta}\}.\tag{13}
$$

The complexity of applying the CO preconditioner is *O n* 4/3 as well. The preconditioners in Equations (7) and (12) were incorporated into the BiCGStab iterative solver [29]. We used the complex-valued version of the solver with the standard complex-valued dot product.

Practical controlled-source modeling requires excessive gridding near the source location. To avoid this, we preferred secondary field modeling, i.e., the electric field *E*(*x*, *y*, *z*) is represented as a sum *Ea*(*x*, *y*, *z*) + *E<sup>b</sup>* (*x*, *y*, *z*), where *E<sup>b</sup>* (*x*, *y*, *z*) is the response due to background conductivity model known analytically [30,31]. In this case, the secondary field *Ea*(*x*, *y*, *z*) will be the unknown function and its FD discretization is performed. The actual source *F*(*x*, *y*, *z*) in Equation (2) is substituted with the secondary source σ*a*(*x*, *y*, *z*)*E<sup>b</sup>* (*x*, *y*, *z*).

#### *2.2. Applicability of Quasi-Statinary Simulation*

As the first step, we tested applicability of the quasi-stationary simulation in the context of controlled-source electromagnetics. The displacement current in the air must be considered when the electromagnetic field generated by a high-frequency dipole oscillator is measured at large offsets [10,12]. Otherwise observed components of the electromagnetic field cannot be matched to the computed quasi-stationary ones; this is known as the propagation effect. Formally, it is achieved by replacing *i*ωµ0σ term in Equation (1) with *i*ωµ0σ − ω <sup>2</sup>µ0ε0ε *rel*, where ε<sup>0</sup> is the vacuum dielectric permittivity and ε *rel* ≥ 1 is the relative dielectric permittivity. It complicates the 3D numerical simulation considerably. However, in contrast to the individual components, the surface displacement-current impedance was almost identical to the quasi-stationary one, at least in a 1D Earth. To illustrate this point, we considered an X-oriented horizontal electric dipole (HED) located at the origin of Cartesian coordinate system and two receiver stations (Figure 2).

The first station corresponded to relatively short offset (*X* = 450 m, *Y* = 750 m) and the second one imitated a typical CSAMT offset (*X* = 450 m, *Y* = 4500 m). We computed the electromagnetic field for the two cases. In the first case, the model was the one in Table 1. Computations were conducted in the quasi-stationary regime, where the squared wavenumber of *i*-th layer was given by *i*ωµ0σ*<sup>i</sup>* . The air conductivity was set to σ<sup>0</sup> = 10−<sup>8</sup> S/m. In the second case, the squared wavenumbers of each *i*-th layer were defined as *i*ωµ0σ*<sup>i</sup>* − ω <sup>2</sup>µ0ε0ε *rel i* . We set ε *rel* <sup>0</sup> = 1 in the air, and ε *rel i* = 4, *i* = 1, 2, . . . in the Earth. The computations were performed by the 1D code of [13]. Computed curves are presented in Figure 2. At frequencies higher than 7 kHz (the close receiver) and 27 kHz (the distant receiver), the quasi-stationary electromagnetic field (dashed lines) differed considerably from the displacement-current field (solid lines). However, ratios of the horizontal component remained essentially the same.

**Table 1.** The background 1D resistivity model.


 ≥ 1

0 0 −

(, , )(, , )

<sup>2</sup>00

(, , ) (, , )

0

– – Cagniard's apparent **Figure 2.** Impact of the displacement currents on the electromagnetic field and the surface impedance on the top of a layered Earth. The left column of panels illustrates computations relative to a short transmitter–receiver separation (receiver Rx-1). The right column of panels presents computations for a long transmitter–receiver separation (receiver Rx-2). Panels (**a**,**b**) depict geometry of the simulation, the top view. An X-oriented horizontal electric dipole (HED) was placed at the origin. Receiver station Rx-1 was located at (450, 750) m, whereas receiver station Rx-2 had coordinates (450, 4500) m. Panels (**c**,**d**) depict magnitudes of individual components of the electric field. Panels (**e**,**f**) depict Cagniard's apparent resistivity (ρ*a*). Panels (**g**,**h**) show impedance phase. In the legends, QS stands for the quasi-stationary computations.

0

= 1

<sup>0</sup> 10 −8

<sup>2</sup>00

0 −

  0

= 4, = 1,2, …

#### *2.3. Implementation*

Forward modeling of a single source involves the following steps:


Our implementation was in C++. The most computationally demanding steps are the right-hand side (RHS) and secondary field computation. To reduce the computational burden, these steps were parallelized using OpenMP shared memory implementation. Simulation of multiple sources and frequencies was parallelized using the Message Passing Interface MPI. In [32], we demonstrated good scalability of this scheme. Consortium for Electromagnetic Modeling and Inversion's

Ωm

#### **3. Numerical Experiments**

#### *3.1. Code Verification on Simple Models*

Our code in general has been validated and benchmarked against analytical layered-earth plane-wave solutions in [23]. Verification of the code in the low-frequency controlled-source electromagnetics scenario was performed in [32]. Here, we conducted several tests concerning particularities of the CSAMT/CSRMT, mainly high frequencies and high conductivity contrast. **Ωm**

−

We used Consortium for Electromagnetic Modeling and Inversion's software PIE3D for benchmarking. PIE3D was rigorously verified and has been used in production for years [25,33,34]. −∞ ∞

In this test, we used a resistivity model consisting of a 1D background model and a 100 Ωm parallelepiped body. The background model is presented in Table 1. The parallelepiped body had coordinates of the two opposite corners (110, −28, 20) m and (240, 28, 52) m; thus, that was a brick 132 × 56 × 32 m<sup>3</sup> with its top at 20 m (Figure 3). ∞

**Figure 3.** The 3D model used for benchmarking. The grey cube depicts position of the receiver station. The sphere is placed at transmitter position.

− The transmitter was an X-orientated point electric dipole located at (32.54, −553.5) m. The receiver was at (200, 80) m. In this test, we used seven frequencies: 0.1, 0.2, 0.5, 1, 2, 5, and 10 kHz.

Finite-difference computations were performed with large and fine numerical grids. The core domain was the same at each frequency 478 × 158 × 150 m<sup>3</sup> . It included both the anomalous body and receiver location (Figure 4).

**Figure 4.** Geometry of the computation grid. The core domain (black rectangle) covers the anomalous body (red rectangle) and receiver position (small black square). The small black circle depicts position of the transmitter.

The grid step size was identical in all frequencies, *h* = 2 m. The padding part of the grids was different for different frequencies. The overall physical size of the grids varied from 56.0 × 56.0 × 28.1 km<sup>3</sup> at 0.1 kHz to 6.7 × 6.4 × 3.3 km<sup>3</sup> at 10 kHz. The number of the internal edges was between 17 M at 10 kHz and 23 M at 0.1 kHz.

PIE3D ran with grid step size *h* = 2, meaning that the body was discretized into 66 × 28 × 16 cubical cells. Computed electromagnetic responses are compared in Figure 5.

difference software (our code), respectively. "Re" and "Im" mean the real and software. "IE" and "FD" mean the integral difference software (our code), respectively. "Re" and "Im" mean the real and **Figure 5.** Benchmark against PIE3D software. "IE" and "FD" mean the integral-equation software (PIE3D) and the finite-difference software (our code), respectively. "Re" and "Im" mean the real and imaginary parts, respectively. Note that the difference curves are magnified by factor 10.

We observed a close match between individual components. The normalized misfit between IE and FD responses, that is k*uIE* − *uFD*k/k*uIE*k, varied from 0.2% to 0.6%. This is encouraging, taking in mind that the computations were done by two very different approaches with no shared code. ‖ − ‖/‖‖ ‖ − ‖/‖‖

#### *3.2. Conductivity Model of Aleksadrovka*

We evaluated our numerical code on a realistic 3D model using an acquisition geometry from a real geophysical survey. The model mimicked composition of Moscow State University's geophysical test camp near the village of Aleksandrovka in Kaluga Region, Russia. The acquisition geometry is presented in (Figure 6). In this work, we did not compare modelled versus measured data. composition of Moscow State University's composition of Moscow State University's

**Figure 6.** Map of the study area from which the acquisition geometry was taken for the numerical simulation. The boxes depict receiver stations. The receiver lines are numbered from 1 to 8.

The receiver grid consisted of 192 receiver stations and two orthogonal transmitter lines (see Figures 6 and 7). All receivers had azimuth 68.5 ◦ . Coordinates of receiver stations and transmitter lines can be found in the Supplementary Materials File S1. The five components of the electromagnetic field were simulated.

**Figure 7.** Topography of the top of the first layer of limestones. The white squares depict receiver stations. The white lines are grounded transmitters lines. Arrows Ex and Hy indicate orientation of the local coordinate axes at receiver stations. Coordinates of the acquisition system are attached to the electronic version of this article.

The model used in this study was compiled from various geologic studies, geophysical surveys, including seismics, and drilling data. Aleksandrovka area is located in Moscow syneclise and has relatively simple morphology of pre-Quaternary layers, but the shallow part of the subsurface has a more complicated structure due to Moscow glaciation [35]. The reference 1D geoelectrical model is based on a 300 m borehole drilled in the camp for water supply purposes. The lithology column is shown in Figure 8.

— — — — — — — — — — — — **Figure 8.** The lithologic column from the borehole shown in Figure 1. Lithology explanation: 1—Sands (QIV), 2—loams (QIV), 3—Interbedding of sand, clay and limestone (C<sup>1</sup> ), 4—Interbedding of sand and clay (C<sup>1</sup> ), 5—Interbedding of sand and clay with lenses of coal (C<sup>1</sup> ), 6—Interbedding of clay and coal (C<sup>1</sup> ), 7—Limestone (C<sup>1</sup> ), 8—Dense green clay(C<sup>1</sup> ), 9—dolomite (D<sup>3</sup> ), 10—Interbedding of limestone, dolomite and gypsum (D<sup>3</sup> ), 11—Interbedding of dolomite and marl (D<sup>3</sup> ), 12—Limestone (D<sup>3</sup> ).

The top part of the column (above 93 m) is composed of mostly clayey sediments with a resistivity of approximately 20 Ωm except for the shallowest Quaternary sands. The sands have a resistivity of about 500 Ωm. There are two thin (~10 m) layers of limestones and dolomites at the depth of 93 m and 115 m. Their resistivity is estimated to 5000 Ωm. The depth of the top of the first carbonate layer was clearly detected by CSRMT across the area. It is underlain by Famennian rocks that consist mostly of dolomites and marls with the overall thickness of 140 m. This formation has a resistivity of about 10 Ωm. Since this estimate came from AMT data and the Famennian formation contains several thin dolomite layers known to be resistive within the region, this resistivity estimate can be regarded as the lateral resistivity. Below this layer, the well entered hard limestones showing its upper 30 m. We estimated the thickness of this layer as 30 m and resistivity as 1000 Ωm, based on the results of the AMT inversion and numerical simulation experiments. Properties of underlying rocks down to 1 km have been obtained from 1D inversion of a single AMT curve. The depth of the Famennian rocks was estimated as 480 m. A very conductive layer (1.5 Ωm) in the range from 480 to 730 m, as well as a resistive (630 Ωm) basement below 730 m, were introduced by the 1D inversion. The final 1D reference model is presented in Table 2.


**Table 2.** The reference 1D resistivity model.

The 1D horizontally layered reference model was deformed based on the topography of the top of the first layer of resistive carbonates. This reference surface was clearly identified from previous CSRMT survey. The most prominent feature of it is the trough oriented along the NW–SE direction, which has the maximal relative depression of about 60 m (Figure 8). Other layers with thicknesses taken from the 1D model conform to the reference surface.

Finally, we added to the model inhomogeneities that represent composition of the shallow part of the subsurface. The subsurface has a relatively complicated geoelectrical structure because of Moscow glaciation. This part has been characterized previously with a detailed CSRMT survey in the frequency range 5–1000 kHz (the far-field responses only) followed by a conventional 2D plane-wave inversion. Several superficial high-resistive sand bodies were delineated on the interpreted 2D cross-sections (Figure 9).

The biggest body was composed of Quaternary glacial or alluvial sands of 160 Ωm with its top at a depth of 20 m. The eastern part of it has an elongated shape, with approximately 100 m along its short axis and thickness of 20 m. This body was inserted into the 3D model. Finally, a 6-m-thick layer of 19 Ωm was introduced across the top of the 3D resistivity model.

The final 3D model is given in Figure 10a,b together with receivers and transmitters. We believe it represents all main features of the Aleksandrovka site. The conductivity model in the Visualization Toolkit VTK legacy file format (https://vtk.org/wp-content/uploads/2015/04/file-formats.pdf) is available in the Supplementary Materials File S1 attached to this article.

**Figure 9.** An example of resistivity cross-sections employed to create the shallow part of Aleksandrovka conductivity model. The cross-sections were obtained previously as a result of the 2D inversion of a CSRMT survey.

**Figure 10.** (**a**) The 3D resistivity model used in the numerical experiment (the shallowest 6 m were removed for visualization). The receivers are depicted as spheres, the two transmitters lines are marked by boxes. Resistivity is shown in color. (**b**) Structure of the shallow geology. The topography of the top of the limestones is shown in color (blue to red). The sand body buried at a depth of 20 m is shown in solid green. The model is attached to the electronic version of this article.

#### *3.3. Numerical Simulation of Aleksadrovka*

We compared efficiency of the sequential 3D modelling code both for GF (7) and CO (12) preconditioners. The transmitter was a point dipole located at the southern end of transmitter line Tx-2 (see Figure 3), with coordinates (1091.101, 812.101, 0.002) and azimuth 155.6 ◦ . We selected three frequencies, 192, 320, and 576 Hz. Parameters of the FD grids are presented in Table 3.


**Table 3.** Parameters of the finite-difference grids.

−10 −10 The code was running on a single thread of a compute node equipped with double 12-core Intel Xeon E5-2680v3 processors running at 2.5 GHz and 128 GB of memory. The relative accuracy of the BiCGStab was set to 10 −10 . Performance of the two preconditioned iterative solvers is presented in Table 4. The table also includes CPU time needed to prepare the RHS.


**Table 4.** Comparison of Green's function (GF) and contraction operator (CO) preconditioned iterative solvers.

(\*) Did not converge in 5000 iterations.

We observed an acceleration of 7×, 9×, and 10× at the three selected frequencies, respectively, both with respect to the iteration count and wall-clock time. Computational load per iteration for the CO preconditioner was only 7–10% larger compared to the GF one. Thus, a reduced number of iterations almost directly translated to the shorter compute time. The GF solver did not converge at 576 Hz frequency for 5000 iterations, which was set as the limit. We attributed slow converge of the GF solver to the highly conductive 1.5 Ωm layer which produced a strong horizontal resistivity contrast in the model. It should be mentioned that topography slowed the computations considerably since performance of the GF preconditioner deteriorated most severely when the air–ground interface was not a horizontal plane. In this case, the time gain of the CO over GF preconditioner was expected to be even greater than one observed in this test.

It is interesting that the time of computing the right-hand side dominated by computation of the background electric field was comparable to the time spent on the solution of a system of linear equations. This is one of the disadvantages of the secondary-field approach. However, RHS computation was easily reduced by using shorted Hankel transform filters and multithreading (see Section 2.2). In the next tests, we applied the CO preconditioned iterative solver only.

In the second set of computations, the whole survey was simulated. The electromagnetic field was computed at 192 receivers. We used 18 frequencies with range from 192 Hz to 550 kHz. It should be emphasized that the electromagnetic field above 25 kHz was not quasi-stationary for our setup, but impedance computed from individual components still can be used to interpret the real measurements (see [13]). Two transmitter lines were used to compute the two polarizations. During computations, the lines were broken into 26 straight segments (15 for Tx-1 and 11 for Tx-2), which had lengths between 40 and 185 m, and then numerically integrated along the transmitter lines. Thus, the total number of individual forward problems was 468. Parameters of the FD grids at different frequencies are given in Table 5.


**Table 5.** Computational grids at different frequencies.

Simulation was conducted on 26 nodes equipped with double 12-core Intel Xeon E5-2680v3 (24 cores per node) processors running at 2.5 GHz, 128 GB RAM, and InfiniBand interconnect. Each compute node was processing the 18 forward problems in the frequency range 192 Hz–550 kHz. There was message passing between nodes during computations. At each node, all 24 cores were working on each forward problem using OpenMP parallelization. The overall computations lasted for 24.5 h. Duration of the individual computation phases is given in Table 6.


**Table 6.** Performance of the forward simulation with OpenMP parallelization on 24 cores.

We observed that computation of the right-hand side (essentially, the background electric field) was comparable to the time allocated to the solution of a linear system. At frequencies above 150 kHz, the right-hand side dominated the other computations. It was caused by a combination of large FD grids and excellent performance of the CO preconditioned solver at the higher frequencies, at which core domains were very thin and enclosed within the top part of the model with moderate

conductivity contrast. It is, however, of minor concern for inverse problems, because the background electromagnetic filed can be stored on disk. The cost of its computation is amortized across many forward problems solved over the course of the inverse problem.

Apparent resistivity and impedance phase along receive line No. 6 are presented in Figure 11. The two polarizations resulted in very similar cross-sections. Apparent resistivity approached the limit of 18 Ωm at frequencies above 25 kHz, indicating the high-frequency electromagnetic field did not penetrate below the top layer.

**Figure 11.** Results of the numerical simulation. Pseudo cross-sections along receiver line No. 6. Top two panels: XY and YX apparent resistivity. Bottom two panels: Impedance phase components, XY and YX.

Maps of the XY apparent resistivity and impedance phase are presented in Figure 12. The map at 192 Hz had a clear correlation with the topography of the first limestone layer (Figure 8), whereas the 1.5 Hz data clearly indicated the high-resistivity sand body (Figure 10b).

**Figure 12.** Results of numerical simulation. Maps of XY apparent resistivity and XY impedance component at several frequencies.

#### **4. Discussion**

In this paper, we considered an effective approach to modeling the data acquired by controlled-source audio-magnetotellurics (CSAMT) and radio-magnetotellurics (CSRMT) methods, which are widely used in mineral exploration. We assessed performance of two preconditioned iterative solvers for CSAMT/CSRMT simulation introduced earlier in [21]. Further, we presented results of numerical modeling of a real CSRMT survey performed near Aleksandrovka, Kaluga Region, Russia. The 3D model we used was based on 2D and 1D inversion as well as drilling data.

Collectively, this study demonstrated that the CO preconditioner is highly efficient even in extreme conditions of high excitation frequencies and large conductivity contrasts, encountered in CSAMT/CSRMT surveys. In contrast, the GF preconditioner degraded severely at the higher frequencies. This observation has not been published earlier as far as the authors are concerned.

– From our experience, extremely large FD grids (above 50 M discrete unknowns) are easily encountered at high frequencies typical for CSRMT (above 10–50 kHz, depending on the subsurface resistivity). However, these conditions were pathological in the sense that the quasi-stationary regime is not valid anymore, and the displacement currents must be considered. The only reason such computations make sense is the quasi-stationary impedance was similar to the displacement-current one, at least in the 1D Earth, as demonstrated by [13].

From the geophysical standpoint, the data at frequencies above, perhaps, 25 kHz were redundant, because the 3D conductivity model did not contain very shallow inhomogeneities (<10 m). Our purpose was to evaluate numerical grids and simulation time if such high frequencies will be needed in the future.

Finally, we emphasize that the results of our study apply to other geophysical methods that require quasi-stationary electromagnetic modelling, specifically MT, AMT, CSMT, CSAMT, CSEM, borehole methods, and others.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2075-163X/10/1/42/s1, File S1: Aleksandrovka model in the VTK format.

**Author Contributions:** N.Y., investigation, software, and original draft preparation; M.M., software, computations, visualization, and original draft preparation; A.S., formal analysis, visualization, and original draft preparation. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Russian Science Foundation, grant number No. 16-11-10188.

**Acknowledgments:** The authors thank M.S. Zhdanov for the opportunity to benchmark our code with PIE3D software, and M. Cuma for their help in performing such calculations. This work has been carried out using ˇ computing resources of the federal collective usage center Complex for Simulation and Data Processing for Mega-science Facilities at NRC "Kurchatov Institute", http://ckp.nrcki.ru/. We would like to acknowledge the Skoltech CDISE's high-performance computing cluster, Zhores [36], for providing the computing resources that have contributed to the results reported herein.

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

#### **Appendix A**

In this appendix, we show equivalence of Equations (7) and (12). Before moving on, denote **Σ***<sup>a</sup>* = **Σ** − **Σ***<sup>b</sup>* which is a discrete secondary conductivity. We will start with Equation (12). After switching from ˆ*e* to *e* and multiplying by **Σ** − 1 2 *b* from the left, we receive,

$$\boldsymbol{\Sigma}\_{\boldsymbol{b}}^{-\frac{1}{2}}\boldsymbol{\mathsf{K}}\_{1}\boldsymbol{e}=\boldsymbol{\Sigma}\_{\boldsymbol{b}}^{-\frac{1}{2}}\boldsymbol{\mathcal{G}}\_{\boldsymbol{b}}^{M}\,\,\,\mathbf{K}\_{2}\boldsymbol{e}+i\boldsymbol{\omega}\mu\_{0}\boldsymbol{\mathsf{A}}\_{\boldsymbol{b}}^{-1}\boldsymbol{f}.\tag{A1}$$

Next, we substitute *K*<sup>1</sup> and G *M b* and employ commutativity of diagonal matrices,

$$\frac{1}{2}(\boldsymbol{\Sigma} + \boldsymbol{\Sigma}\_b)\boldsymbol{\Sigma}\_b^{-1}\boldsymbol{e} = \left(2i\omega\mu\_0\boldsymbol{\Lambda}\_b^{-1}\boldsymbol{\Sigma}\_b^{\frac{1}{2}} + \boldsymbol{\Sigma}\_b^{-\frac{1}{2}}\right)\boldsymbol{\mathsf{K}}\_2\boldsymbol{e} + i\omega\mu\_0\boldsymbol{\mathsf{A}}\_b^{-1}\boldsymbol{f}.\tag{A2}$$

Now, we substitute *K*<sup>2</sup> and rewrite the matrices,

$$\frac{1}{2}(2\boldsymbol{\Sigma}\_{b} + \boldsymbol{\Sigma}\_{a})\boldsymbol{\Sigma}\_{b}^{-1}\boldsymbol{e} = \left(2i\omega\mu\_{0}\boldsymbol{\Lambda}\_{b}^{-1}\boldsymbol{\Sigma}\_{b}^{\frac{1}{2}} + \boldsymbol{\Sigma}\_{b}^{-\frac{1}{2}}\right)\frac{1}{2}\boldsymbol{\Sigma}\_{a}\boldsymbol{\Sigma}\_{b}^{-\frac{1}{2}}\boldsymbol{e} + i\omega\mu\_{0}\boldsymbol{\Lambda}\_{b}^{-1}\boldsymbol{f}.\tag{A3}$$

It is easy to note that the term <sup>1</sup> 2 **Σ***a***Σ** −1 *b e* cancels out, and we receive,

$$\mathbf{e}\left(\mathbf{I} - i\omega\mu\_0 \mathbf{A}\_b^{-1} \Sigma\_a \mathbf{e}\right) \mathbf{e} = i\omega\mu\_0 \mathbf{A}\_b^{-1} \mathbf{f}.\tag{A4}$$

After factoring out *A* −1 *b* in the left-hand side and noting that *A* = *A<sup>b</sup>* − *i*ωµ0**Σ***a*, we obtain

$$\mathbf{A}\_b^{-1} \mathbf{A} \; \mathbf{e} = i \omega \mu\_0 \mathbf{A}\_b^{-1} \mathbf{f},\tag{A5}$$

which is exactly (7).

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
