Geoelectrical Tomography Investigating and Modeling of Fractures Network around Bittit Spring (Middle Atlas, Morocco)

Direct current Resistivity (DCR) method was carried out to characterize the hydrogeological connection between the Tabular Middle Atlas (TMA) and the Sa¨ıs Basin. The TMA is one of the most important aquifers in northern Morocco that supplies the deep aquifer of the Sa¨ıs Basin. Electrical resistivity tomography (ERT) survey was focused on the Bittit area that is one of the most important outlet discharges, and it is representative of the relations between the TMA and the Sa¨ıs Basin. The high resolution capabilities of the electrical tomography were used to deﬁne the geological draining features in the framework of water paths from the TMA to the karstic springs. The synthetic data were calculated for the similar model expected in ﬁeld data inversion and inversion result of these synthetic data used as a guide for the interpretation of the inverse data resistivity sections. Joint interpretation of geophysical, geological, structural, and synthetic simulation data allowed identifying a conductive horizontal shallow layer overlying two subvertical families of fractures of NE-SW and NW-SE directions. This result leads to propose hydrological behaviour of water from the Tabular Middle Atlas and the Sa¨ıs Basin at the Bittit Spring, which takes into account for both horizontal ﬂows through stratiﬁcation joints or karst and through subvertical fractures.


Introduction
Water resources and management is one of the major environmental issues in Morocco. The climate changes (decrease of rainfalls) and the increase of water needs for human activities (farming, industry, individual needs) make critical the management of water resources. The precipitations are mainly located in mountains while water requirements are located in the plains: a good understanding of water transfer between these areas is absolutely necessary. The Middle Atlas is a mountain long of about 350 km from south-western to north-eastern Morocco, located between the Rif and High Atlas. The mountain has cold winters with persistent snow above 2000 m. The western border is well watered (from 1000 to 1500 mm annually) because of the meteorological perturbations coming from the Atlantic Ocean. The Tabular Middle Atlas (TMA) located at the NW border of the Middle Atlas is made up of Liassic limestone and dolomite forming a complex karstic system locally highly fractured (Figure 1(a)). It is one of the most important reservoirs in northern Morocco because of the climatic context and of its structure. The water flows essentially along open fractures and karst cavities, with very low matrix permeability [1]. Although it represents a large groundwater potential, it is still badly exploited because of the lack of data about fracture network, aquifer geometry, and recharge area. It makes difficult to define optimally the potential areas of drilling. This reservoir is connected to large aquifer located the in Saïs plain, at north of the TMA. This aquifer can supply water to two large cities (Fez and Meknes) and small villages. The Liassic reservoir of TMA dips in the Saïs Basin under a thick Tertiary marly series, forming an important confined aquifer (Figures 1(a) and 1(b)). Groundwater circulation from TMA to the Saïs Basin is complex due to fracturing and 2 International Journal of Geophysics to karstic Quaternary travertine overlying Lias carbonates. The geometry of the aquifer, water flow paths, and the complex underground karstic organization are misknown due to lack of high resolution geophysical survey. Therefore, a better understanding of hydrogeological behaviour of fracture and karst networks is fundamental for management and protection of water resources.
The goal of this study is to characterize the transition between the reservoir of the TMA where the water is accumulated and the Saïs Basin where the water is used. In this paper, we focus on the underground imaging and on mapping of fractures network in the Bittit Spring area that constitutes one important rising spring at the junction between the TMA and the Saïs Basin, with a great average annual flow of about 1600 l/s (Figure 1(b)). Therefore, it is considered as a representative area of the hydrogeological connections between the TMA and the Saïs plain. To high-light these issues, we choose a geophysical imaging by electrical Resistivity Tomography (ERT), which is an important prospecting tool adapted to groundwater research. Recently, this technique has been largely applied in geophysical investigation of areas having complex geology [2,3] it allows imaging the geometry of underground structures with high resolution and can detect the vertical and lateral electrical resistivity variations induced by rock properties and water content. After a complete description of the studied area, we described the ERT data acquisition and results. The field data are interpreted by using 2D inversion algorithm after topographic correction and a shape filtering. We performed a synthetic modelling in order to refine interpretations of the ERT. Finally, a 3D view was realized to check the geometrical compatibility between the various detected structures and then to define a model of draining hydraulic features.

Geological and Hydrological Settings
The Middle Atlas and the Rif dominate the landscape of Northern Morocco. The Middle Atlas is separated from the southern ridges of the Rif belt by the Saïs plain. Mesocenozoic sediments covering Paleozoic units mainly form this plain and the Middle Atlas. They are deformed during the Mesozoic extension and Cenozoic compression. At North of Middle Atlas, the Tabular Middle Atlas (TMA) corresponds to an elevated tabular structure fractured during the Cenozoic.
The TMA is a large karstic plateau whose elevation ranges between 1000 m and 2200 m. It is mostly made up of Liassic limestones and dolomites overlying clays and basalts of the Triassic and the Paleozoic formations. The Quaternary is represented by alkaline basalt lava flows coming from the Outgui volcano and going to the Saïs plain and by travertine frequently located at the emergences of springs at the junction with the Saïs Basin.
TMA is affected by major NE-SW flexures clearly visible except under basaltic flows and by an intense fracturing with main NE-SW and NW-SE directions. The discharge of TMA aquifer occurs essentially in a complex system of springs in its northern border at the junction with the Saïs Basin [4].
The Saïs plain is an asymmetrical E-W sedimentary basin of 2100 km 2 average total surface where elevation ranges between 500 m and 700 m. It is limited by the Pre-Rif ridges in the north, by TMA in the south, by the Rharb plain in the west and by the Fez-Taza corridor in the east. The Saïs Basin mainly consists of thick Miocene marls reaching 1 km near Fez, overlying, deep carbonate (limestones and dolomites) Lias aquifer and covered by the Plio-Quaternary lacustrine limestone and sandstone that constitute a water table aquifer of less importance. The recent quaternary is represented by basalt flows, alluvial cones, and travertine on the TMA edge [5].
The geophysical surveys carried out in this basin revealed thick Miocene marls at north, several flexures affecting Pliocene-Miocene deposits, and normal faults affecting the Liassic carbonates. The flexures result from the replay accidents of the Palaeozoic bedrock and are also influenced by the Rif tectonic [6]. Several boreholes and geophysical surveys show variations in thicknesses and depths of Lias that is the most important aquifer [7]. The Lias reservoir is affected by several faults generating discontinuous blocks which are either raised or collapsed; by consequence, this deep aquifer is discontinuous [5]. Two aquifers are present in the Saïs Basin: one shallower (<80 m) in the Pliocene-Quaternary formations and a second deeper (up to 1000 m) in the Liassic carbonates. The monitoring of the deep aquifer in the Saïs Basin displays a regular decrease of the water level due to the drought and the overpumping [7].
The Liassic formations are the loci of the water storage and of the circulation both in the TMA and in the Saïs plain. However, groundwater circulations between these two structures are complex because faulting, karst connection, and dolomite geometry control it. Up to now, water flow paths and underground karsts organization remain unknown.
In order to characterize this transition between the TMA and the Saïs Basin, we investigated the Bittit Spring area located at 32 km southeast of Meknes, at the junction between the Tabular Middle Atlas and the Saïs Basin ( Figure 2). The Bittit Spring is one of the most important in term of flowing with an average annual discharge of about 1600 l/s. Amraoui, 2005 [7] concluded that the hydrogeological system of Bittit has significant storage, sustained discharge during the year, and relative slow flow. However, the spring discharge decreases regularly as a result of a long drought.
We carried out a 100 m borehole close to the Bittit Spring to identify precisely the stratigraphy around the spring ( Figure 3). This information will be used to interpret the DCR data. From bottom to top, the borehole showed 57 m of Triassic clays, 33 m of a succession of limestones and marls of Liassic, dolomites of middle Liassic, a large unconformity, and at least a travertine limestone covered by the basalt flow. However, this borehole cannot be generalized on the whole study area, because other boreholes show different lithological facies of Quaternary formations, the presence of Miocene marls, and different depths of the Lias carbonates location. The faulting is probably the main reason of these depth variations.

2D Inversion of DCR Tomography Data.
Electrical prospecting methods are particularly well adapted to geoelectrical discontinuities, faults, and water investigation [8]. Among these methods, Electrical Resistivity Tomography (ERT) allows to conduct fast data acquisition and to obtain underground 2D images with high resolution suitable for fracture identification [9].
Determination of the underground geometrical structure with electrical methods depends on the resistivity contrasts of the underground rocks. Electrical resistivity of rocks varies according to several parameters like porosity, weathering and water content. Chapellier and Mari [10] assigned wide ranges of electrical resistivity for various classes of rocks.
In the absence of laboratory measurements of in situ samples, local vertical electrical soundings can be used to reduce the range of resistivity corresponding to the lithologies of our studied area. In the TMA, Benslimane [11]  Seven ERT profiles were carried out around the Bittit Spring ( Figure 4). Apparent resistivity data were acquired by an ABEM Terrameter SAS1000 and a multielectrode LUND imaging system using a Wenner array configuration of 64 electrodes and 5 m spacing. It allows investigating a depth of about 50 m. The lengths of the profiles vary from 315 m to 395 m with a rollup of the cables. Geological map of the Bittit area provides major NE-SW and NW-SE fault direction that helped to define the direction of the profiles. Some profiles have also been crossed to obtain resistivity information in perpendicular directions.
2D inversion and topographic correction of the pseudosections were performed using the Res2DInv [13] software package. Outliers were removed, and topographic corrections were carried out with measurements of elevation and geographic coordinates for each electrode position through an inclinometer and a GPS. Finally, all the pseudosections were inverted by the least squares method with the same inversion parameters and with 5 iterations, giving root mean International Journal of Geophysics square value less than 6% between the data and the calculated pseudosections.

ERT Inversion
Results. The inversion of the pseudosections data was performed using the Res2DInv software package. The resulting sections display various electrical resistivities ranging from less than 10 Ωm to more than 1000 Ωm ( Figure 5). Observing the inversion results, we have classified the resistivity values into three ranges of resistivities: [ show the presence of vertical structures and heterogeneous and complex aspect. First, the Bittit 7 profile was carried out perpendicularly to the NE-SW profiles in order to verify the presence of the upper conducting layer and its extension towards the TMA. The result shows that this structure seems to be almost continuous from the TMA at the SE to the spring at the NW. We observe also a very resistant block (∼1000 Ωm) with a length of approximately 80 m, limited by conducting vertical structures on the Bittit 5 and Bittit 7 sections. It should be noticed a lack of some data due to the using the roll-up on profiles of Bittit 4, 5, 6, and 7, in a zone width of about 55 m and centered at x = 265 m for Bittit 4, 5, 6 and x = 125 m for Bittit 7, at depths between 40 m to 47.5 m from the ground surface. Therefore, the lateral extent of the resistive blocks and the resistivity value on these profiles must be taken with caution. Nevertheless, we can observe that these resistive blocks are individualised and dislocated by vertical conducting structures on Bittit 2, Bittit 5, and Bittit 7 sections and also by a heterogeneous distribution of the resistivities, especially on Bittit 1 and Bittit 5 sections.
To sum up, the main electrical structures consist of (1) a shallow subhorizontal conductive layer and (2) a deeper resistive layer, discontinuous due to subvertical conducting structures and heterogeneous resistivity media.

Filtering (V/H FF).
The resistivity inversion provides a first image of the geometry of the underground that  smoothes the edge of the geological contacts. In order to limit this effect and to define more precisely the limits of horizontal layers or to improve the detection of the subvertical features, we used the Res2DInv vertical-to-horizontal flatness filter ratio (V/H FF). This filter fixes a relative weight between a vertical flatness filter (fz) and a horizontal flatness filter (fx) in order to take into account the a priori information about the shape of the geological structures if any. By default, the same factor is used for both. Thus, we selected a higher value for the ratio of the V/H FF in order to force the program to produce models that are vertically stretched.  We have applied this filter with a ratio of 4 for all profiles.
In this paper, we show only the results for two of them, representative of the obtained results. Figure 6 shows structures with a higher apparent dip than that observed on the unfiltered sections. In some cases, new vertical structures appear (see Bittit 5 on Figure 6(d)) which were not clearly detected on the first inversion results or which were on the edges of the section and are likely to be confused with edge effects.
We have also applied this filter using a ratio of 0.6 in order to produce models with horizontally stretched structures. This filtering slightly decreases the response of some vertical structures which continue to appear with almost the same dip as in the unfiltered sections. It also removes some vertical structures in the case of Bittit 1, 4, and 6 which became completely horizontal. Then, the vertical structures that persist in both cases are interpreted as vertical geological features of our study area.
Although this filter does not allow a precise quantitative interpretation, for example, of the dip angle or of the structures extension, it enables us to validate the presence of subvertical structures that were finally used on the geological and hydrological interpretation.

Synthetic Simulation
In order to illustrate how a particular geometrical resistivity structures is resolved in the ERT data inversion and then to improve the understanding of the inverted resistivity profiles, a synthetic simulation approach is proposed to model the electrical behaviour of complex structures. In this simulation, the Res2DMod software package [14] was used to compute the apparent resistivity pseudosection thanks to a geometrical model of different resistivity bodies. This method is based on finite difference modeling schemes [15]. The synthetic apparent resistivity data were then inverted with Res2DInv in order to verify how the software recovers the initial model. We used the same parameters than field survey: Wenner array, 5 m electrode spacing and the maximum n-level is 15a = 75 m.
We present two examples related to the two main structural aspects observed on our data sets: one with horizontal layers characterised by the Bittit 4 section, the other one with a vertical structure surrounded by horizontal layers characterised by the Bittit 2 section. In both cases, the final model was constructed after several steps in order to fix the best values of some parameters like number of layers, depth, thicknesses, resistivities. Each time the simulation results were compared to surveyed sections. At the end, we choose the model that gaves a resistivity pseudosection that correctly fit the data. We limit our investigation to simple structures and limited range of resistivities. As a consequence, only the major resistivity anomalies of the data sections are taken into account.

The "Bittit 4" Model.
The inverted section of the Bittit 4 data shows almost tabular layers with contrasted resistivity values ( Figure 5). The first step of this simulation was to determine the number of layers necessary to fit the best with the Bittit 4 pseudosection. The choice of three horizontal layers with contrasted resistivities was based on lithological informations of the study area. It was confirmed by numerical tests that additional layers did not improve the results. Thus, resistivity value of 10 Ωm corresponds to the electrical behaviour of soil when located at the surface and of saturated travertine at depth. Resistivity value of 100 Ωm corresponds to travertine or wet limestone. Resistivity value of 1000 Ωm represents Liassic carbonates and basalt lavas.
The model corresponding to the Bittit 4 profile consists of a travertine layer of about 28 m thick, locally overlain by thin superficial structures, alternating conductive and resistant materials in order to take into account the thin superficial resistivity anomalies observed on the Bittit 4 data section. We showed that this thickness value does not modify significantly the resulting pseudosection when variations of ±2 m are considered. This travertine lies on the non-fractured Liassic limestone of high resistivity. As the computed apparent pseudosection does not fit well the low resistivity anomaly located approximately at 770 m of elevation, that is, 15 to 20 m depth (Figure 7), it was necessary to insert a low-resistivity layer that represents a travertine layer with higher water content.
Many depths have been tested to define the place of this very conductive zone in the travertine layer, form its basement to 10 m depth. Finally the depth of 12 m has appeared to give the best fit. Some thicknesses have also been tested and led to a 4 m thick layer.
The Res2DInv inversion of the synthetic data ( Figure  7(b)) is globally in good agreement with the inversion of the original data pseudosection ( Figure 5), even if some anomalies of small extension of the data are not well reproduced by the model, particularly on the border part of the profile. Nevertheless, this model confirms that the data can be explained by almost continuous horizontal layers, with resistivity values representative of each lithological formation observed in this area. It can also be noticed that the resistivity values and depth of the different layers are well recovered by the Res2DInv inversion, even if the interfaces between layers are smoothed. However, the geometrical interfaces are deformed by the inversion, probably because of edge effect. This must be taken into account in the interpretation of inverted resistivity profiles.

The "Bittit 2"
Model. The inverted section of the Bittit 2 data shows more complex geometry than in the Bittit 4 profile, with subhorizontal layers intersected by one vertical conductive structure ( Figure 5). We choose to model this particular pseudosection because it shows vertical conductive features of major interest in our study and because it is the less heterogeneous one with this kind of structure.
The choice of the number of layers and their resistivity is based on the results obtained for the Bittit 4 section. Their thicknesses were adjusted in order to fit the best with the Bittit 2 data. A vertical feature has been added to simulate a conductive zone such as a fractured structure (Figure 8). Different parameters were tested: the resistivity value, the width, the depth of the top of the fault, and its dip. A resistivity value of 40 Ωm for the fault zone has seemed to be the best value to model the data. Indeed, lower resistivity produced a high-conductivity zone at the base of the inverted model that is not representative of the real data, while greater resistivity values produced a resistant block. A first value of the width of the fault zone was determined from the result of the inverted section of the data and then adjusted by comparison between the model and the data. A variation less than ± 4 m does not significantly modify the results. The depth of top of the fault zone was set to 8 m, underlying the conductive superficial layer, almost continuous along the profile.
Finally, the "Bittit 2" model shows three layers increasingly resistive with depth like that in the Bittit 4 model. The vertical structure offsetting the deep resistant layer can be simulated by a vertical layer of 25 m width with resistivity of 40 Ωm starting from 8 m depth (Figure 8), interpreted as a fractured zone with high conductivity.
International Journal of Geophysics  0  20  40  60  80  100  120  140  160  180  200  220  240  260  280  300 Resistivity model  The aim of these synthetic simulations was to constrain the interpretation of the results provided by Res2DInv. It must be kept in mind that different resistivity models can give very similar pseudosections. We made use of this synthetic simulation to also validate the usefulness of the V/HFF Res2DInv filter in the processing (Figure 9). The inversion of the pseudosection model corresponding to the Bittit 2 data, applying the V/HFF filter with a dumping factor of 4, shows that the geometry of the vertical fault zone is better determined, with a representation of the dip of the fault closer to the true vertical one.
Finally, we have seen that synthetic simulation allowed us to obtain the general geometric structures of the underground to better interpret our data set, even if some parameters like dip of vertical features can not be precisely estimated and to check that resistivity anomalies observed on the inversion result are due to geological structures and not to inversion effects.

Interpretations and Discussion
The results obtained in this study confirm the efficiency of the electrical tomography in hydrogeological application and, in our case, to the location of the water paths between the TMA and the Bittit Spring area. Several resistivity features were observed that can be classified into two sets: horizontal layers located at shallow depth with low resistivity and more complex and contrasted resistivity pattern at depth. These features are interpreted either in terms of lithology or in terms of water content. The borehole associated to the electrical panel provides a strong constrain for our interpretation to determine the resistivity range for main lithologies. Variations inside each main facies were associated to changes in water saturation and/or porosity variations. Localised fracturing or a change in lithology properties induces these local variations. The deep parts (>15 m) of the electrical sections display important lateral variation in resistivity over a short distance. These variations limit blocks with contrasted resistivity: this setup is interpreted as geological blocks limited by subvertical borders. These later structures are mainly located at depth and contrasts with shallow horizontal levels with low resistivity. By consequence, subvertical zones affecting the deep parts of the investigated area are sealed by the tabular units and are interpreted as fault zones. This observation is consistent with surface observations that show no significant faults at surface.
A simple structural map ( Figure 11) was drawn based on correlations of different structures identified on the resistivity sections (shallow low-resistivity layers and lateral resistivity variations). The spatial consistency of the structures was first checked using the geomodeler GOCAD ( Figure 10). In particular, we validated the fit of the resistivity units at the intersection between two resistivity sections, as for instance at the intersection between Bittit 1 and Bittit 2 profiles that revealed the presence of the conductive structure. The different inferred faults were modelled into GOCAD to better constrain the interpolation between the distant resistivity profiles (Figure 10).
To determine the water path, we reported on map ( Figure 11) detected resistivity features with points for deep lateral changes and with grey polygons for low-resistivity layer. This map reveals a large shallow layer of low resistivity that becomes thinner from the eastern part of the survey to the western part where it arises to the surface. It disappears on the north border, close to the Bittit Spring.
The resistivity values are consistent with porous travertine filled by water. This interpretation fit to field observations that pointed out large amount of travertine in the area of the Bittit Spring. The subvertical faulted zones are organised into three trends: ENE-SWS, NW-SE, and E-W directions. These faulted features are deep and are never observed at surface. They often limit blocks with medium resistivity values. Taken into account the resistivity values, these faulted zones do not correspond to very active water paths. They form a network of uplifted and collapsed blocks. The NE-SW direction corresponds to a main trend that affects the TMA. The NW-SE and EW directions were observed on the outcrops of the northern border of the TMA (Figures  11 and 13).
To sum up, the ERT survey showed a subhorizontal conductive zone attributed to porous and saturated travertine that suggests a main drainage layer at shallow depth, mainly located at the south-east of the spring. Furthermore, vertical features seem not to be main water drains. Other geological argues lead to the same conclusions. One of the faults extracted from the Bittit 5 profile could supply the Bittit Spring by extrapolating the fault plane toward NE (Figures 11 and  12). However, the Bittit 5 profile is located on a hill, 50 m higher than the spring. If this fault was draining, water should emerge at the intersection between topography and the fault plane on the hill slope, but no source rises at this International Journal of Geophysics  point. It must be assumed that these subvertical faults are not draining features. All these data confirm the geological evolution of the Bittit region and allow proposing a preliminary local hydrological model. The Bittit area is formed by Liassic basement at a depth around 20-40 m ( Figure 13). This basement is affected by two sets of normal faults generated during two different tectonic events: One associated to a NE-SW to N-S extension and the second one to NW-SE extension. All these geological trends are seen on the electrical sections. These two tectonic phases were previously described in Morocco and belong to a main lower Triassic event [16][17][18][19].
From a hydrological point of view, we have observed water draining zones coming from the TMA to the Bittit Spring that seem to be essentially localized on a shallow subhorizontal level of porous travertine. Locally, no subvertical features seem to be major water paths. Unfortunately, the connection between the Liassic reservoir of the TMA and the deep aquifer of the Saïs plain could not be shown in this study because it seems to be deeper than the investigation depth of this ERT survey. However, it is likely that this connection is governed by complex hydrodynamic behaviour which takes into account both horizontal flow through stratification joints or karsts and flow through subvertical fractures or porous and permeable fractures zones.

Conclusion
This study aims to determine the water drains supplying the Bittit Spring, one the most important rising spring on the Saïs plain. It will constrain the hydrological connection from TMA to Saïs plain. We choose an electrical investigation taken into account for the expected investigation depth (some tens of meters) and the objects to detect (water flow). The obtained results confirm the efficiency of the electrical tomography in this kind of applications.
The TMA represents the water reservoir of the Saïs basin aquifer. One of the major questions concerns the hydraulic connections between these two systems: does it take place through vertical faults or within subhorizontal stratigraphic layers? In order to answer this question, we focused our study to the particular case of the Bittit spring, rising between the TMA and the Saïs plain.
The Electrical Resistivity Tomography investigation carried out in this area consists of seven profiles covering an area of 5 km 2 . The inversion of these pseudosections was performed using Res2DInv software. A Vertical to Horizontal flatness filter was applied to discriminate subvertical resistivity anomalies that can be associated to geological and hydrological features. A synthetic modelling of two of these pseudosections was performed to precise the main structural aspects observed on the data sets, that is, horizontal conductive layers and vertical less conductive features. Then, a 3D view of the results allows proposing an interpretative map of the area, including subvertical features and conductive areas. It reveals the presence of horizontal drainage layer at shallow depth at SE between the TMA and the Bittit spring. The conductivity of this layer should correspond to porous and saturated travertine.
Two subvertical fault families with NE-SW and NW-SE trends respectively were identified and correspond to the known geological evolution of the region. This fractures network highlighted dislocation of the blocks of underground substratum. It was shown that the NE-SW faults can not be related to the spring and thus do not constitute draining features. These observations lead to suppose that the main water paths from the TMA to the Bittit spring are located inside a shallow horizontal layer of porous travertine. The dislocation of the deeper Liassic blocks should prevent water to flow continuously.
This electrical investigation does not allow identifying relations between the TMA and the deep aquifer of the Saïs plain, but it appears that water flow has a complex behavior taking into account all the water drain possibilities (fractures, karsts and stratification joints).
This study initiates a hydrogeophysical research in the Middle Atlas karst in order to improve water resources management and reducing aquifer vulnerability in the region.