INTRODUCTION
Population growth, the development of different economic practices (such as livestock and agriculture), the inadequate use of land, construction in sectors that are not suitable for a certain purpose, the accelerated and growing deforestation in last decades and the constant pollution from various water sources, among other factors, generate every day, the increase in the demand for water supply. The Food and Agriculture Organization (FAO), the specialized agency of the United Nations, notes that in 2010 the world population corresponded to 7000 million inhabitants and the extraction of water exceeded 4000 km3 per year. It is estimated that currently 54% of freshwater is consumed and by the middle of the 21st century the world population will reach 12000 million inhabitants, the demand for water will have doubled and the water reserves of our planet will reach their peak (FAO, 2016). The study area, located on the W side of the Eastern Cordillera on the southern end (edge E) of the Santander Massif and the S sector of the Almorzadero moor, presents a strong problem caused by the shortage of water resources. The present surface water sources are insufficient to meet the requirements and supply the basic needs (Salamanca, 2017). In addition, studies on the existence, location and availability of this resources are incomplete or fragmented and, in general, sometimes there is no information on the water resources of the area; only data about surface water, which is in low volumes, generally contaminated and unfit for human consumption is presented. However, the geographic and geological location of the study area, along with the medium rainfall of the sector, allow some drains that make up the water network to have constant flows and many others intermittent, generating optimal and favorable conditions for the occurrence of groundwater (Tarazona and Rueda, 2018). The geophysical prospecting is presented as the most appropriate technique to be applied, as a quick, economic and non-destructive method that allows the identification and characterization of the subsoil, achieving the construction of descriptive models that allow analyzing the behavior of certain physical properties (such as resistivity), necessary for its subsequent (and correct) interpretation (Loke, 2012). Within the diversity of geophysical prospecting techniques, geoelectric methods stand out. The use of geoelectric prospecting for hydrogeological purposes has been widely documented. Geological and hydrological studies in lagoons reefs in Lake Muri (Cook Islands), was performed by the acquisition of 28 Electrical Resistivity Tomography (ERT) lines, using a SuperStingR8 tomograph (Befus et al., 2014). The studies of resistivity developed in the southwest of Sohag (Upper Egypt) determined the impact of the geological environment on the occurrence of groundwater, by characterizing (qualitatively and quantitatively) descriptive geoelectric profiles (Mahmoud and Tawfik, 2015). The subsurface stratigraphy and the structures around the Dakhla oasis (Egypt) were determined through the interpretation of 10 Vertical Electric Soundings (VES) with the Schlumberger distribution and the IPI2Win software (Mohamaden, 2016). Geoelectric studies have been used for decades for hydrogeological purposes, in the generation of 1D profiles and 2D subsoil models, used in the description and characterization of the subsurface, as well as their geological, geoelectric and/or geophysical conditions (Asfahani, 2007). The effectiveness of these methods consists in the cartography of underground structures that integrate or constitute the subsoil (Kumar, 2012) and determine the presence of aquifers (Chirindja et al., 2016). However, this geophysical method is also widely used in mining, in geotechnical investigations and currently in environmental studies (Loke, 2012). Synthetic data processing is carried out to determine the best strategy for the information gathering and the ideal characteristics for the inverse modeling, without the need of costly trial and error tests in the field (AGI, 2009). Used by several authors (Díaz, 2010), (Hernández, 2012) the processing is applied on synthetic profiles of 2D electrical resistivity, in order to obtain the optimal parameters for the acquisition of data in the field. The main scope of this work aims to build, from applied geophysical prospecting: 15 VES and 4 ERT, the 2D model of electrical resistivity of the subsoil, through the processing of 4 synthetic models for the treatment of information geoelectric for hydrogeological purposes. The results obtained here will be of vital importance, taking into account that it is intended to contribute, improve and deepen the knowledge and description of the subsoil to provide fundamental information on aquifer resources in the perimeter area of Málaga, in the department of Santander. The information above serves as technical - scientific support to the different authorities (environmental, territorial, local and educational), becoming the starting point for future studies of greater detail.
GEOLOGICAL SETTING
The area of interest is located in the NE of Colombia, on the Eastern Cordillera, in the E sector of the Department of Santander, in the jurisdiction of Province of García Rovira, including the urban area of Málaga and its periphery. The study area is distributed between southern Málaga and northern San José de Miranda (Figure 1). The sector is constituted by sedimentary rocks belonging (according to the nomenclature of the Maracaibo Basin) to the Cretaceous Río Negro (Kirn) and Tibú - Mercedes (Kitm) formations and Quaternary deposits (alluvial Qal and colluvial Qc). The Río Negro Formation consists in a very variable set of dark, fissile, laminated shales interspersed with fine-grained quartzitic sandstones and some conglomerates, with a high degree of fracturing. The Tibú - Mercedes Formation lithologically corresponds to gray fossiliferous limestone intercalations, shale and gray sandstones of fine to medium grain. Angular fragments of rock within a sandy matrix (Qc - colluvial deposits) and gravelly, sandy and muddy levels (Qal - alluvial deposits) constitute the Quaternary deposits (Vargas et al., 1976). The structural behavior is dominated by compressive stresses, where inverse faulting, fracturing, complex folding and asymmetric structures are the most outstanding features. The main fault system corresponds to the one generated by the Río Servitá fault, an inverse structure of regional type that presents a sinuous line running N - S and controls the course of the river of the same name (Alcaldía municipal de San José de Miranda, 2003). Meanwhile, the folding system is represented by symmetric and asymmetric structures of NE - SW direction, where the Málaga syncline corresponds to its maximum representative, affected by systems of parallel faults (Málaga - Miranda fault) and transverse (Aguablanca fault) to the Río Servitá fault (Alcaldía municipal de Málaga, 2015).
MATERIALS AND METHODS
A total of four geological models were generated in representative sectors, where geographic and geological conditions were favorable for the accumulation of groundwater; and from these models, 4 synthetic profiles of 2D electrical resistivity were developed, using CorelDRAW (2012 version) and ArcGIS (version 10.1) software and the digital elevation model (DEM) resolution 12.5 by 12.5 m, downloaded from the Vertex platform, which is the remote sensing data portal of NASA (UAF, 2018). In the processing stage about the synthetic models of electrical resistivity was applied direct modeling and inverse simulation. In the first process (direct simulation), performed by the finite element method (Díaz, 2010), the pseudosections of synthetic apparent resistivities were generated, while in the inverse modeling the Smooth Model Inversion method was used (AGI, 2009) to obtain the inverted resistivity sections. The optimal parameters, applied in the real data, were obtained by comparing the initial 2D synthetic models with the inverted resistivity profiles, using a code elaborated in the MATLAB software (version R2017a) by means of the resistivity value point-to-point comparison. In the process of data acquisition, a total of 15 VES and 4 ERT were developed. Out of the 15 VES performed, 11 of them were developed by means of the Schlumberger configuration, where the distance AB/2 varied from 1 to 150 m and the value of MN/2 started at 0.25 m and ended at 10 m. The remaining 4 VES were made through the Wenner matrix, where the distance between the electrodes was modified from 3 to 60 m. For the acquisition of each survey a terrameter SAS 1000 of ABEM was used. On the other hand, the 4 ERT were acquired with 56 electrodes using the Wenner- Schlumberger arrangement, where the separation between electrodes in the T1 survey corresponded to 2 m (110 m of extension in surface), in the T2 it was 10 m (total extension of 550 m), T3 was performed with 3 m (165 m of extension in surface) while the electrodes of the ERT 4 were 5 m apart (275 m of extension in surface) by means of a SuperSting R8/IP/SP tomograph. For each measurement, a GPS Garmin brand was used for georeferencing, a Brunton type compass for the azimuth direction record, a level and a topography tool to perform the required topographical correction, as well as a polyethylene hammer, external battery and tape measure. In the processing stage, the data obtained by each developed VES were processed using the IPI2Win software (Bobachev et al., 2000), and the information acquired by each ERT was processed using the AGI Earthimager 2D software (AGI, 2009), according to the optimal parameters determined in the processing of synthetic models.
RESULTS
Synthetic models
The synthetic models of electrical resistivity 2D were made with the objective to get the optimal parameters for the acquisition of data in the field (configuration of the electrodes) and for the inversion process (number of iteration and the damping factor). According to the above, and starting from the initial configuration of the study area, 4 geological models were generated in representative sectors, and from them, 4 synthetic profiles of 2D electrical resistivity were developed (T1- T4). The representative sectors are: the NE area where the UIS (T1) and Airport (T2) profiles were distributed on the axis of the Málaga syncline to evaluate the influence of this structure on the accumulation of groundwater; the S sector in the Guásimo - Naranjitos region (T3) where a model was created on recent deposits (alluvial type) which were generated by the main drainage (Chorrerón or Malagueña stream); the E sector in the Barzal area (T4) where a model was generated outside the study area to analyze the influence of the adjacent region. The four synthetic profiles mentioned above were generated by the AGI Earthimager 2D software, where it was simulated with 56 electrodes, 10 m apart (total extension of 550 m) and the type of instrument used corresponded to SuperSting R8 (parameters of the tomograph to be used). The arrangement of the electrodes is the first parameter to vary within the synthetic profiles. Dipole- Dipole, Wenner-Schlumberger and Schlumberger were the selected distributions, that means 3 files were created for each generated geological profile. In the process of direct modeling (Forward Simulation), the software generates the division or discretization of the subsoil in rectangular blocks using a grid where each cell has a specific resistivity and after that, executes a virtual acquisition, generating a set of synthetic data, represented by means of a pseudo-profile of apparent resistivity calculated. The apparent resistivity of the model is determined from the potential difference in the nodes of the rectangular grid by applying the finite element method (Figure 2 - 3).
For the inverse modeling, it was necessary to configure the process before its execution, modifying the parameters, according to the objectives of this report. Thus, the number of iterations varied from 1 to 20 while the value of the damping factor was modified from 0.5 to 10. Using the Survey Planner option in the tools menu, the software (AGI Earthimager 2D) generated the inverted resistivity profiles, according to the selected configurations. Each model was processed by inverse simulation by varying the previously described parameters and their respective generated sections were saved in a file format (X, Y, Z) with extension .dat. in a folder called 1 filter. The saved profiles correspond to iterations 5, 10, 15 and 20 with damping factor 0.5, 3, 7 and 10. According to the above, for each developed synthetic model, 16 inverted resistivity sections were created for each configuration, 48 in total.
Determination of optimal parameters
To determine the optimal parameters, the inverted resistivity profiles were compared with the initial synthetic models using a code elaborated in the MATLAB software (version R2017a) by means of the resistivity value pointto-point comparison, for its analysis. The software selects the profile (and its parameters) that represents the best investment. Once the profile and its optimal parameters were selected, and with the aim of reducing and limiting the range in the variation of each of them, a second filter was applied, where the number of iterations and the damping factor are modified one more time.
Four values were added (two units above and two units below) to the value determined in the first filter. That means the second filter consisted of five values for each parameter (Table 1). If the value of 0.5 was selected for dumping, the variation was made two tenths above and two tenths below, for a total of five values. The second filter consisted of applying the same process (direct and inverse modeling) on the original synthetic models, modifying the new values in the number of iterations and in the damping factor.
According to the above, 25 inverted resistivity sections were obtained for each of the three configurations, 75 in total, stored in a file in format (X, Y, Z) with extension .dat., in a folder called 2 filter. These generated pseudoprofiles were compared, again, with the original synthetic models to obtain the optimal parameters of the second filter, using the code created in the software MATLAB (version R2017a) mentioned above. Thus, the iteration number and the damping factor of each configuration are obtained, to be applied in the investment process of the acquired data. To determine the best electrode arrangement (Dipole-Dipole, Wenner-Schlumberger and Schlumberger) to be used in prospecting, the parameters selected in the second filter are compared between them, by means of a second code generated in the MATLAB software (version R2017a). The code compares the original synthetic model and the selected optimal parameters of the second filter for each configuration and determines the least error between them (Table 2), selecting the best matrix to be applied in the field. The selected configuration corresponded to Wenner-Schlumberger.
According to the above, for each geological model, a synthetic profile of resistivity was generated, three pseudosections of apparent resistivity calculated (one for each configuration), 48 inverted resistivity models for the first filter (16 for each matrix) and 75 inverted resistivity profiles for the second filter (25 for each arrangement). That means the code developed in the MATLAB software selects the best model among 123 inverted resistivity profiles determining the best configuration of electrodes and optimal investment parameters.
Data acquisition
The acquisition of data was based on taking information directly from the land, which complemented or modified the data previously collected. It focused on the deduction of geological characteristics (distribution of structures present in the subsoil) and geophysical (variations in the electrical resistivity of these structures) through the use of resistivity methods (VES and ERT) and according to the results obtained in the stage previous methodological process. For the acquisition of each VES a terrameter SAS 1000 of ABEM was used, while the ERT were obtained by means of the SuperSting R8/IP/SP tomograph with a maximum extension of 550 meters. A total of 15 VES were developed in 5 representative sectors of the study area, with the aim of using the data obtained here as support information or control points for each of the ERT carried out. The representative sectors were: in the NE sector where 4 surveys were carried out with Schlumberger distribution (S1 - RUIS1A, S2 - RUIS2A, S3 - RUIS3A and S4 - RUIS4A) over colluvial deposits; the central area with 4 soundings, two of them using Schlumberger type configurations (S5 - RCREA1 and S6 - RCREA2) and the others with two Wenner matrices (S7 - RCREW1 and S8 - RREW2) over alluvial deposits; the N sector with two Schlumberger arrangement (S9 - ROVI1 and S10 - ROVI2) over colluvial deposits; the S zone where 1 VES was developed, with Schlumberger distribution (S11 - RYERB) over alluvial deposits; and the W area with 4 surveys, two of them with Schlumberger arrangement (S12 - ELM1 and S13 - ELM2) and the others two Wenner type (S14 - WELM1 and S15 - WELM2) over colluvial deposits. The distribution can be seen in the map of Figure 1. The ETRs were acquired in the same sectors where the synthetic models were made, however road access, topographic restrictions and private farms determined the maximum extension. These sectors were: the NE area called T1 or UIS and T2 or Airport, the S area called T3 or Guácimo - Naranjitos and sector E designated T4 or Barzal, as shown in the map of Figure 1. All the previous geophysical prospection were developed with Wenner-Schlumberger configuration (matrix obtained in synthetic models).
Data processing and interpretation
The processing of the information acquired in the field corresponded to the treatment applied to the data that was collected or obtained. Within this phase, the true electrical resistivities of the materials present in the subsoil were estimated and calculated, from the theory of investment (Loke and Barker, 1995), by specialized software for this purpose and according to the optimal parameters obtained in the synthetic models. The data obtained by each of the developed VES will be processed using the IPI2Win software (Bobachev et al., 2000). On the other hand, the information acquired by each ERT was processed using AGI Earthimager 2D software (AGI, 2009).
Electrical Resistivity Tomography (ERT)
The processing were carried out with the optimal parameters of filter 2 of Table 1, with a good correlation, where the RMS error varies from 0.94% to 7.23% and L2 started from 0.10 to 5.81.
The ERT 1 designated as UIS was developed on the E flank of the Málaga syncline, orientation N - S. It was performed at a height of 2227 m.a.s.l. in the N margin and 2220.8 m in the S sector, which corresponds to its lowest topographic point. With a separation between electrodes of 2 m (total extension of 110 m) and an estimated depth of 27.9 m, the tomography recorded resistivity values from 6.5 to 1435 Ω*m, with an RMS error of 7.23% and L2 equal to 5.81 (Figure 4A). The processing and interpretation generated a profile with three layers of electrical resistivity. The first one has an average thickness of 5 m with a resistivity that varies from 234 to 1435 Ω*m interpreted as a Quaternary deposit of colluvial type in dry conditions. It underlies a layer of variable thickness from 7 to 13 m with low values of resistivity, from 6.5 to 60.5 Ω*m, analyzed as wet clay material belonging to deposits of recent age of colluvial type. The last layer is at least 16 meters thick (the total thickness is unknown) with resistivities of 60.5-234 Ω*m, corresponding to calcareous material (highly fractured and saturated biomicrite limestones) belonging to the upper part of the Cretaceous Tibú - Mercedes Formation (Kitm). This calcareous material, at a distance between 18 and 24 m from the initial point of the tomography and at 14 m of depth, presents a fracture, which is filled with clayey material in humid to saturated conditions (Figure 4B).
The geophysical prospection 2 named Airport was carried out on the axis of the Málaga syncline, perpendicular to this structure. It presents a NW - SE orientation at a height between 2208.3 m.a.s.l. at the NW and 2209 m.a.s.l. at the SE. It has a total extension in surface of 550 m (separation of 10 m between the electrodes) and total depth of 117 m. The resistivity registers a wide variation in its value, from 1 to 5000 Ω*m, with an RMS error equal to 3.91% and L2 of 1.69 (Figure 5A). A 3-layer profile was generated, analyzed and interpreted from the acquired data. The most superficial zone of the subsoil presents very low resistivities (< 39.7 Ω*m) where the NW sector registers the lowest value of the present study corresponding to 1 Ω*m, in a layer that decreases its thickness from 29 m to disappear at 450-460 m from the initial point of the ERT. The aforementioned is interpreted as fine, sandy and muddy materials, belonging to Quaternary deposits of alluvial type in saturated conditions. On the surface and at a distance between 70 and 80 m from the first electrode there is a small area with a high resistivity (5000 Ω*m), the bed of Los Molinos creek, canalized drainage, dry and used as a material accumulation point of anthropogenic waste in the construction of the Airport runway, origin of the high resistivity value. Underlies these materials, alluvial Quaternary deposits, similar to the upper one, but in humid conditions, with resistivities from 8.4 to 71 Ω*m and 45 m of average thickness. However, this layer is not recorded on the SE side and is only observable at the NW end. The last layer corresponds to a material with resistivities from 71 to 333 Ω*m and thicknesses > 70 m (unknown total thickness) and is defined as highly fractured and saturated biomicrite limestones belonging to the upper part of the Cretaceous Tibú - Mercedes Formation (Kitm) (Figure 5B).
The tomography 3 called Guácimo - Naranjitos was obtained in the S sector of the work area, on recent deposits (alluvial type) generated by the main drainage Chorrerón or Malagueña stream. With a SW - NE orientation, the maximum height corresponded to 2184.3 m.a.s.l. and the lowest level was 2176.1 m.a.s.l. The separation of electrodes was established in 3 m, total extension of 165 m and 39.6 m depth. Registers resistivity values from 10.2 to 514 Ω*m, with an RMS error of 2.19% and L2 equal to 0.53 (Figure 6A). After processing, 3 layers of electrical resistivities were interpreted. The most superficial part of the subsoil corresponds to fine, sandy and muddy material, belonging to recent deposits of alluvial type in dry conditions, with variable thickness (from 5 to 15 m) and with the highest resistivity values of the present prospecting (from 132 to 514 Ω*m); however, it is only noticeable in the central sector and on the NE side of the ERT. After this, a layer with resistivity values of 10.2 Ω*m up to 50 Ω*m and variable thicknesses, from 7 to 20 m, is interpreted as wet alluvial deposits, increasing the degree of humidity towards the SW region. On these materials lies a layer with a thickness > 26 m (the total thickness is unknown) and average resistivity of 72 Ω*m, analyzed as calcareous material (highly fractured and saturated biomicrite limestones) belonging to the upper part of the Cretaceous Tibú - Mercedes Formation (Kitm) (Figure 6B).
The ERT 4 named Barzal was acquired in sector E ofthe work zone with SSE - NNW orientation at a height of 1850.9 m.a.s.l. in the SSE area, corresponding to thetopographically lowest point of the present survey and 1856 m.a.s.l. at the NNW end. The separation of the electrodes was every 5 m, with a total extension of 275 m and an approximate depth of 62 m. The resistivity obtained ranges from 14 to 155 Ω*m, with an RMS error equal to 0.94% and L2 of 0.10 (Figure 7A). The processed data generated 3 layers of electrical resistivity, where the first one has a resistivity lower than 26.5 Ω*m, with a maximum thickness of 16 m, analyzed as a Quaternary deposit of the saturated colluvial type, however it is only observable on the NNW side and not on the SSE. The underlying layer corresponds to moist clay material belonging to deposits of recent age of the colluvial type with a variable thickness, from 18 to 30 m. Underlying lithological material with average resistivities of 46.6 Ω*m and thickness of at least 30 m (unknown total thickness) which is defined as highly fractured and saturated biomicrite limestones belonging to the upper part of the Cretaceous Tibú - Mercedes Formation (Kitm) (Figure 7B).
Vertical Electric Soundings (VES)
A total of 15 VES were developed in 5 representative sectors of the study area, with the aim of using the data obtained here as support information or control points for each of the ERT carried out. The data obtained by each of the developed VES were processed using the IPI2Win software, with a good correlation, the error varies from 0.66 to 5.49%.
To exemplify the processing and interpretation, the VES called WELM1 is presented. It was acquired on the right margin at kilometer 2.5 of the road that leads from the Málaga town to the urban area of Molagavita, 70 m before the intersection with the La Magnolia stream. The coordinates correspond to X =1232640, Y =1147445, Z =2448 m.a.s.l. This geoelectric prospecting was carried out with the Wenner matrix where the separation of the electrodes “a” started at 6 m and ended at 60 m, with a total extension in surface of 180 m. Using the IPI2Win software, a 4-layer profile with an error of 1.17% was obtained (Figure 8).
The interpretation of the VES called WELM1 is related to colluvial deposits (with high water content), and calcareous material (biomicrite limestones) with different fracturing conditions and saturation degrees, belonging to the upper part of the Cretaceous Tibú - Mercedes Formation (Kitm) (Table 3).
To visualize the spatial and in depth-variations of the resistivities of the subsoil, correlations were made between the different VES, developing 2D geoelectric profiles from the union of two 1D models, which are presented in Figure 9.
DISCUSSION
The present work reached a maximum depth of 117 m relating the subsoil with the physical properties of the materials that make it up. In the study area, there are two types of soil with a high water content and thicknesses below 3 m: the first type is transported, located on the flanks of the Málaga syncline and the other is residual type, on the core of the structure. The alluvial and colluvial deposits, distributed heterogeneously, have different degrees of humidity, from saturated with resistivity below 20 Ω*m, to dry with records above 150 Ω*m.
The calcareous material corresponds to saturated and highly fractured biomicritic limestones where the resistivity value varies from 20 to 400 Ω*m and thicknesses from 20 to 70 m with intercalations of fractured and humid biomicritic limestones with resistivity records from 400 to 600 Ω*m, where the total thickness is unknown and compact and dry micritical calcareous material where the resistivity value is above 600 Ω*m and unknown thickness (Table 4).
The hydrogeological interest layer corresponds to the calcareous material, saturated and highly fractured biomicritic limestones with secondary porosity (thicknesses from 20 to 70 m and resistivities from 20 to 400 Ω*m) belonging (according to the nomenclature of the Maracaibo Basin) to the upper part of the Cretaceous Tibú - Mercedes Formation (Kitm).
It is recommended to carry out hydrogeological studies to determine the nature of groundwater in the delimited area of this study, as well as its diffusion, movement, regime and reserves. Also, it is necessary to perform an exploratory well to execute a pumping test in order to analyze the hydraulic properties of the hydrogeological interest layer as well as the physical, chemical, organoleptic, bacteriological and radioactive properties of the same. This must be done with the objective of determining the conditions of use, regulation and evacuation or extraction of groundwater present in the work area to improve the social conditions of the community and boost the economic growth of the same.
CONCLUSIONS
In general, the geophysical information of the subsurface and geoelectric models generated allows a more robust characterization of the subsoil in the perimeter area of Málaga town, in the department of Santander (Colombia). The methodology applied is useful to establish a representative geoelectric model of a study area and it could be applied in a zone with absence of information about the subsurface and with water capacity. In the present case of study, the results indicate that the construction or implementation of 2D synthetic profiles improves the processing of real resistivity data in a geoelectric prospecting. The different resistivities obtained in the processes of interpretation and analysis of the data acquired in the field, allow infer that the optimum and favorable conditions for the occurrence of groundwater are presented. The hydrogeological interest layer corresponds to calcareous material, biomicrite limestones saturated and fractured, with resistivities from 20 to 400 Ω*m and thicknesses of 20-70 m, belonging to the upper part of the Cretaceous Tibú - Mercedes Formation (Kitm). The results obtained contribute to the knowledge of the characteristics of the subsoil in the study area and determine the presence of a confined aquifer isolated from the atmosphere by an impermeable layer constituted by clay - silty, materials highly weathered and low energy deposits, localted on the axis of the Málaga syncline, under the urban zone. It is suggested to carry out a project to avoid the water losses of the sewage system, as well as any discharge that may contaminate the aquifer, which has a high potencial exploitation, to complement the supply to the aqueduct urban.