1. Introduction

The central region of the Isthmus of Panama is characterized by the lowest precipitation rate during the dry season, affecting negatively various human activities, including the agricultural development of the region. In parallel, these activities are also affected by the El Niño-Southern Oscillation phenomenon. One way to deal with the negative effects of these phenomena in this region is to take advantage of the groundwater in a systematic way. To do so, knowledge about the geological structure of the area where this resource will be exploited is needed.

Groundwater potential is determined by detecting and delineating the subsurface fractured layers that can offer pathways to the groundwater flow system (Berkowitz, 2002), while the aquifer characteristics are controlled by the thickness of the weathered and/or fractured rock (Acworth, 2001).

Located in the northern part of the Estibaná River watershed (central region of the Isthmus of Panama), the San Luis area has a surface geology that suggests the presence of a geological formation consisting of altered intrusive rock, due to which it is classified as an important water reservoir.

Obtaining information on subsurface geological formations through, more conventional, geotechnical drilling has certain drawbacks, such as limited budgets, time constraints, and the need for human resources in the field. This results in a reduced number of tests and therefore in a lack of detailed knowledge about the spatial distribution of the rock and other relevant characteristics. Another option is the use of 2D electrical resistivity tomography methods. In fact, 2D and 3D electrical resistivity imaging have been used in several regions to detect, delineate, and characterize intrusive rock domains (e.g., Mendoza and Dahlin, 2008; Rossi et al., 2017; Soro et al., 2017; Gómez et al., 2019; Hasan et al., 2019).

In Panama, electrical resistivity imaging methods have been employed to detect geological contact zones and faults (Mojica et al., 2017), intrusive bedrock, and groundwater resources (Díaz et al., 2012; Castrellón Romero et al., 2019; Mojica et al., 2019, 2021; Rejiba et al., 2021; Ruiz et al., 2021). The present study therefore aims to delineate the altered intrusive rock domain in the San Luis area through surface geological study and 2D electrical resistivity tomography. We hope these new results will help define future models of groundwater management in the region.

2. Geological context

The San Luis area is located in the northern part of the Estibaná River watershed, on the Azuero Peninsula of Central Panama. The area is dissected by the trace of the Ocú-Parita fault, which crosses the entire peninsula in an east-west direction (Fig. 1).

FIG. 1. Location of the San Luis area in the central zone of Panama. Inset shows a simplified geological map of the surveyed zone and surroundings. OPF: Ocú-Parita fault, PFZ: Pedasí fault zone, RJFZ: Río Joaquín fault zone, ACF: Agua Clara fault, ASF: Azuero-Soná fault (modified from Instituto Geográfico Nacional Tommy Guardia, 2007).

The San Luis area has been characterized as a fragment of the Ocú-Parita fault, which directly influences the Estibaná River catchment channel, particularly at its confluence with the La Villa River channel. Topographically, this sector presents a concave valley morphology with numerous intermediate elevations that suggest the presence of fractured blocks with delayed-stepped vertical displacements, characteristic of grabbed-type structures (Ruiz, 2018¹, 2020²).

In the San Luis area, Late Cretaceous sedimentary sequences overlie fractured crystalline shales and a sterile, porphyritic intrusive body with well-developed subparallel hornblende crystals (Ruiz, 2020). The east-west topography and the presence of sedimentary remnants on the summits of the area suggest that the blocks within this fault zone have experienced vertical movements (Ruiz, 2018). North and south of the San Luis area, the displacement is mainly sub-horizontal, exposing the fractured components to intense erosive processes associated with the Pesé and Tonosí formations (Ruiz, 2020). A simplified geological map of the San Luis area is shown in figure 1.

The intrusive bodies of interest originate from the Loma Montuoso batholith, possibly from a calc-alkaline magma that differentiated into a felsic magma, producing quartz diorites and granodiorites. Based on multi-crystal, incremental heating ⁴⁰Ar/³⁹Ar dating, the age of the intrusive suite is 64.87±1.34 Ma (Mann and Corrigan, 1990). Likewise, Wegner et al. (2011) confirm similar ages to other intrusive bodies recognized on the mineralized belt of the Azuero peninsula, also associated with the Loma Montuoso batholith, which exhibits a predominant ~N70W orientation.

The composition of the intrusive body found in the San Luis area suggests a less felsic magma. The development of subparallel hornblende crystals points towards the presence of abundant water in the magmatic system, unlike the other intrusives known in the zone. Tectonic activity was likely the reason for its late rise.

3. Methodology

To shed light on the subsurface distribution of the geological formations in the study area, two sections were established: one (P1) of 787 m in length, in a southwest-northeast direction, and the second (P2) of 784 m in length, in a southeast-northwest direction. On top of each geological section, 235 m-long profiles were designed to carry two 2D electrical resistivity tomographies (T1 and T2). Figure 2 presents the distribution of the geological (P1 and P2) and electrical (T1 and T2) profiles and a topographic map of the San Luis area.

FIG. 2. Distribution of the geological (P1 and P2; red lines) and electrical (T1 and T2; black lines) on top of the topographic map of the San Luis area (contours every 5 m.).

3.1. Field measurements and parameter settings

The electrical resistivity of the soil is determined by various factors, including soil composition (e.g., mineralogy), particle size distribution, and structure (e.g., porosity, connectivity, pore size distribution, water content) (e.g., Bai et al., 2013). In addition, thermodynamic parameters such as temperature and pressure also influence soil resistivity (e.g., Nover, 2005). This relationship between the electrical resistivity and moisture content of soil and rocks makes the electrical resistivity method a valuable tool in hydrogeological studies (e.g., Reynolds, 1997).

Over the past few decades, the development of sophisticated electronic and computerized systems has led to significant progress in geophysical data recording devices. Our study used a Syscal R1 Plus Switch-48 for resistivity profiling. This instrument consists of two multi-cable systems of 24 channels, each with a maximum distance between the electrodes of 5 m. For carrying out the automatic switching of these electrodes and acquiring profiling data, two parameters were set: (i) acquisition parameters, which include the minimum and maximum number of stacks per measure (in this study, 3 and 6 respectively), the quality value (or minimum standard deviation) of the voltage and electrical current ratio for the stacks (1%), and the time cycle (500 ms); and (ii) geometrical parameters, such as the number of electrodes (48), electrode distance (5 m), number of levels (23), and electrode array (Wenner-Schlumberger). Regarding the electrode array, it was considered a suitable option for achieving both good and vertical resolutions (e.g., Loke, 2021³).

3.2. Scheme of inversion

Each data file was inverted using the Boundless Electrical Resistivity Tomography (BERT) code (Günther et al., 2006). This code uses irregular triangle meshes and the finite-element forward calculation discussed by Rücker et al. (2006). BERT uses the smoothness-constrained Gauss-Newton method, whereby the linearized problem is solved as a function of a Jacobian matrix (which contains the partial derivatives of the model response and the model parameter vectors), the updated model vector, the observations (apparent electrical resistivity values), and the model response vector.

The inversion process was carried out through a constrained smooth minimization operation, using an equation that depended on the error-weighted data fit, the model roughness, and a regularization parameter (Ronczka et al., 2017). According to Günther et al. (2006), model roughness is given in terms of a matrix applied to the model parameter vector. The regularization parameter (represented by l) controls the smoothing constraints, so setting a high parameter could generate a smooth model. In our case study, a value of l=2 was established to obtain electrical anomalies within well-defined limits. During the 2D inversion process, the initial model was modified iteratively to reduce the difference between the model response vector and the observations (Loke, 2021). This data fit was evaluated by calculating the chi-square value (χ²) as a weighted error in each iteration.

3.3. Forward modeling

In order to check the validity of the inversion, a forward modeling was conducted in both profiles through a direct simulation by means of the AGI EarthImager 2D software. The forward solution was obtained by solving a partial differential equation, which contained the scalar electric potential in the Fourier transform domain, the intensity of the electric current in the 2D space through the Dirac delta function, the wavenumber in the transformation of the partial differential equation through the Fourier analysis, and the distribution of the electrical conductivity in the 2D space (Dey and Morrison, 1979; AGI, 2007⁴).

The model response vector (synthetic apparent resistivity values) was generated through the finite difference method implemented in the software. The results of the inversion of the field data from each electrical profile determined the number of 2D structures to choose.

4. Results

4.1. Inversion of the data field

The inverted resistivity images of profiles T1 and T2 are shown in figure 3A and B, respectively, with acceptable χ² values of 0.51 and 0.36. Figure 3A shows a superficial conductive horizon represented in blue and green tones, with a range of calculated electrical resistivity values between 4.07 and 27.19 Ωm and depths that exceed 10 m in some areas. Below this horizon, a high electrical resistivity anomaly is identified in red tones with values between 70.28 and 182 Ωm. This latter anomaly extends laterally and in-depth in almost the whole profile. In between these two electrical anomalies, a yellow-colored horizon is characterized by electrical resistivity values between 27.19 and 70.28 Ωm.

FIG. 3. Inverted resistivity images obtained in profiles T1 (A) and T2 (B). See profile location in figure 2. See text for explanations on the symbols used.

Figure 3B shows a superficial conductive horizon in blue and green tones that extends from ~459 to 569 m along the profile, at variable depths (between ~3 and 12 m), and with calculated electrical resistivity values ranging between 4.07 and 27.19 Ωm. Between positions 334 and ~459 m along the profile, this conductive layer reaches greater depths compared to the right side of the figure. The high electrical resistivity anomaly is shown between positions ~459 and 569 m and below the conducting superficial horizon. This electrical anomaly is represented in red tones, with calculated electrical resistivity values ranging from 27.19 to 182 Ωm.

4.2. Geological sections

The inverted resistivity images obtained above were combined with the surface geological information reported by Ruiz (2018, 2020) to create geological sections along the profiles P1 and P2 (Fig. 4).

The interpretation corresponding to the geological structure presented in these profiles indicates the presence of a warped layer (topsoil) of low electrical resistivity that overlies a fractured body of high electrical resistivity. This body is associated with the altered porphyritic mafic intrusive recognized on the surface.

FIG. 4. Geological sections obtained along the profiles P1 (A) and P2 (B). Gray bands delimit the areas covered by the resistivity images shown in figure 3. See profile location in figure 2.

4.3. Synthetic models

In this study, the theoretical models considered homogeneous structures based on previous findings. The models began with a conductive surface layer representing the residual soil with an electrical resistivity value equal to 10 Ωm. By following the information showcased in figure 4, this superficial layer was set to rest on two important structures: (i) a conductive structure of 15 Ωm, corresponding to volcanic deposits, and (ii) a resistant structure of 200 Ωm that represented the altered intrusive body.

Additionally, a 100 Ωm resistivity layer was established in each synthetic model. This layer was set to represent the transition zone between the conductive, superficial soil layer and the resistant, altered intrusive body.

The same inversion process behind figure 3 was followed to develop the synthetic model by establishing the same number of electrodes and parameters. Hence, two 2D electrical tomographies were created (Fig. 5).

FIG. 5. Synthetic models established for profiles T1 (A) and T2 (B), and 2D inversion results for these models (C and D). See profile location in figure 2. See text for explanations on the symbols used.

5. Discussion

The conductive zones (i.e., 4.07-27.19 Ωm) identified in the T1 tomography (Fig. 3A) correspond to a superficial layer of soil, which is the result of the weathering process of the rock formations in the study area. The resistive zones identified in the T1 and T2 tomographies (70.28-182 Ωm) are associated with the altered intrusive body. The horizon located between the conductive and resistive zones mentioned above, characterized by intermediate electrical resistivity values (27.19-70.28 Ωm), can be associated with a transition zone between both horizons. The beginning of the profile in figure 3A, at position ~359 m, shows the southwestern limit of the resistive horizon or anomaly, which could be associated with a fault zone (Fig. 4A).

The conductive surface horizon (i.e., 4.07-27.19 Ωm) presented in the T2 tomography (Fig. 3B) between positions ~459 and 569 m can be associated with the surface soil layer identified in figure 3A. At the beginning of this profile, between positions 334 and ~459 m, the conductive horizon occupies the maximum vertical extension so that it could correspond to the shallow soil layer mentioned above and followed by a second conductive horizon, possibly a layer of volcanic material, with calculated electrical resistivity values <27.19 Ωm. Finally, the resistive anomaly identified in this second profile is associated with the altered intrusive rock body, whose abrupt end at position ~459 along the profile may likely represent the Ocú-Parita fault (Fig. 4B).

Based on the analysis of the existing outcrops in the area, the simplified geological map, the results of the 2D inversion of observations, and the 2D geophysical modeling, the geological interpretation suggests a sequence of Late Cretaceous sedimentary rocks truncated by the regional-scale Ocú-Parita fault (Fig. 4). The distribution of these rock sequences has been constrained in synthetic models to define some relationships between the results of the data inversion and the response of the model. There is an overall good agreement between the in-depth distribution of the two conductive structures (10-15 Ωm), representing the surface soil layer and the volcanic deposits, the two resistant structures (100-200 Ωm), representing the altered intrusive body and a transition zone, the geological structure of the study area, and the 2D electrical resistivity models obtained from the inversion of the field data.

6. Conclusions

The geophysical study conducted in this work made it possible to identify and delimit an altered intrusive body in the San Luis area of Central Panama. Due to its fractured and altered nature, this body can be considered a natural recharge zone and, therefore, key in setting the ground for future groundwater management models. The interpretation of the geological sequences shown in this work has been represented through theoretical 2D electrical resistivity models, whose model response inversion can be adjusted by considering the inversion results from the observations.

Acknowledgements
We wish to thank the Universidad Tecnológica de Panamá and the Instituto de Acueductos y Alcantarillados Nacionales (IDAAN) for all the support. This research was also supported by the Secretaría Nacional de Ciencia, Tecnología e Innovación-SENACYT, Panama (grant number IOMA17-006). We would like to acknowledge the reviews of this paper provided by R. Pinzón and D. Bertin.

¹ Ruiz, A. 2018. Geología del río Estibaná. Reporte Geológico: 46 p. https://www.tierragrande.co/wp-content/uploads/2024/08/REPORTE-FINALrop.pdf

² Ruiz, A. 2020. Consideraciones tectónicas sobre Estibaná y áreas adjuntas. Reporte Técnico: 19 p. https://www.tierragrande.co/wp-content/uploads/2024/08/consideraciones-tectonicas.pdf

³ Loke, M.H. 2021. Tutorial: 2-D and 3-D electrical imaging surveys. www.geotomosoft.com

⁴ AGI, 2007. Instruction Manual for EarthImager 2D version 2.4.0, Advanced Geosciences Inc.: 100 p, Austin.

Bibliography

Acworth, I. 2001. The electrical image method compared with resistivity sounding and electromagnetic profiling for investigation in areas of complex geology: A case study from groundwater investigation in a weathered crystalline rock environment. Exploration Geophysics 32 (2): 119-128. https://doi.org/10.1071/EG01119

Bai, W.; Kong, L.; Guo, A. 2013. Effects of physical properties on electrical conductivity of compacted lateritic soil. Journal of Rock Mechanics and Geotechnical Engineering 5 (5): 406-411. https://doi.org/10.1016/j.jrmge.2013.07.003

Berkowitz, B. 2002. Characterizing flow and transport in fractured geological media: A review. Advances in Water Resources 25 (8-12): 861-884. https://doi.org/10.1016/S0309-1708(02)00042-8

Castrellón Romero, M.G.; Fábrega, J.R.; Foglia, L.; Mojica, A.; Ruiz, A.; Saavedra, S. 2019. Groundwater model of a fractured rock system as a tool for groundwater management: the Estibaná Sub-Catchment, Azuero Peninsula, Panama. In E-Proceedings of the 38^th IAHR World Congress: 10 p. Panama. https://doi.org/10.3850/38WC092019-1794

Dey, A.; Morrison, H.F. 1979. Resistivity modelling for arbitrarily shaped two-dimensional structures. Geophysical Prospecting 27 (1): 106-136. https://doi.org/10.1111/j.1365-2478.1979.tb00961.x

Díaz, I.; Mojica, A.; Ho, C.; Pinzón, R.; Fábrega, J.; Vallester, E.; Vega, D.; Ogden, F.; Hendrickx, J. 2012. Characterization of shallow groundwater in Eocene sediments of Panama Canal Watershed using electrical techniques. I+D Tecnológico 8 (1): 33-42.

Gómez, E.; Svensson, E.; Dahlin, T.; Barmen, G.; Rosberg, J.E. 2019. Tracking of geological structures and detection of hydrothermal intrusion by geo-electrical methods in the highlands of Bolivia. Journal of South American Earth Sciences 91: 214-226. https://doi.org/10.1016/j.jsames.2019.02.002

Günther, T.; Rücker, C.; Spitzer, K. 2006. Three-dimensional modelling and inversion of dc resistivity data incorporating topography - II. Inversion. Geophysical Journal International 166 (2): 506-517. https://doi.org/10.1111/j.1365-246X.2006.03011.x

Hasan, M.; Shang, Y.; Jin, W.; Akhter, G. 2019. Assessment of aquifer vulnerability using integrated geophysical approach in weathered terrains of South China. Open Geosciences 11 (1): 1129-1150. https://doi.org/10.1515/geo-2019-0087

Instituto Geográfico Nacional Tommy Guardia. 2007. Atlas Nacional de la República de Panamá, Cuarta Edición. Novo Art: 290 p. Panama.

Mann, P.; Corrigan, J. 1990. Model of late Neogene deformation in Panama. Geology 18 (6): 558-562. https://doi.org/10.1130/0091-7613(1990)018%3C0558:MFLNDI%3E2.3.CO;2

Mendoza, J.A.; Dahlin, T. 2008. Resistivity imaging in steep and weathered terrains. Near Surface Geophysics 6 (2): 105-112. https://doi.org/10.3997/1873-0604.2007036

Mojica, A.; Pérez, T.; Toral, J.; Miranda, R.; Franceschi, P.; Calderón, C.; Vergara, F. 2017. Shallow electrical resistivity imaging of the Limón fault, Chagres River Watershed, Panama Canal. Journal of Applied Geophysics 138: 135-142. https://doi.org/10.1016/j.jappgeo.2017.01.010

Mojica, A.; Castrellón Romero, M.G.; Ruiz, A.E.; Saavedra, S.; Salceda, M.; Serrano, E.; Fábrega, J.R. 2019. Hydrogeophysical characterization of intrusive bodies in the tropics: Upper Estibaná sub-catchment case study. In AGU Fall Meeting, Abstract: NS41B-0829. San Francisco.

Mojica, A.; Castrellón Romero, M.G.; Ruiz, A.; Fábrega, J.; Saavedra, S.; Salceda, M.; Serrano, E.; Esquivel, A. 2021. Detection of fractured basaltic bodies for aquifer characterization in South Central Panama using 2D electrical resistivity imaging. In AGU Fall Meeting 2021, Abstract: H35C-1052. San Francisco.

Nover, G. 2005. Electrical properties of crustal and mantle rocks-A review of laboratory measurements and their explanation. Surveys in Geophysics 26: 593-651. https://doi.org/10.1007/s10712-005-1759-6

Rejiba, F.; Mojica, A.; Schmutz, M.; Ruiz, A.; Castrellón Romero, M.G.; Ladrón De Guevara, J.; Saavedra, S.; Fábrega, J.; Finco, C.; Schamper, C.; Cavalcante, L.H.; Llubes, M. 2021. Geophysical characterization of near-surface formations in the La Villa River catchment (Los Santos, Panama). In Near Surface Geosciences, 1^st Conference on Hydrogeophysics. European Association of Geoscientists & Engineers. Bordeaux. https://doi.org/10.3997/2214-4609.202120088

Reynolds, J.R. 1997. An introduction to applied and environmental geophysics. John Wiley & Sons: 796 p. Chichester.

Ronczka, M.; Hellman, K.; Günther, T.; Wisén, R.; Dahlin, T. 2017. Electric resistivity and seismic refraction tomography: a challenging joint underwater survey at Äspö Hard Rock Laboratory. Solid Earth 8 (3): 671-682. https://doi.org/10.5194/se-8-671-2017

Rossi, M.; Olsson, P.I.; Johanson, S.; Fiandaca, G.; Preis Bergdahl, D.; Dahlin, T. 2017. Mapping geological structures in bedrock via large‐scale direct current resistivity and time‐domain induced polarization tomography. Near Surface Geophysics 15 (6): 657-667. https://doi.org/10.3997/1873-0604.2017058

Rücker, C.; Günther, T.; Spitzer, K. 2006. Three-dimensional modelling and inversion of dc resistivity data incorporating topography-I. Modelling. Geophysical Journal International 166 (2): 495-505. https://doi.org/10.1111/j.1365-246X.2006.03010.x

Ruiz, A.; Mojica, A.; Castrellón Romero, M.G.; Fábrega, J. 2021. Paleocauce y reducto lacustre La Mesa en el Plioceno. In Congreso Nacional de Ciencia y Tecnología, No. 18, Asociación Panameña para el Avance de la Ciencia (APANAC): 456 p. Panama.

Soro, D.D.; Koïtaa, M.; Biaou, C.A.; Outoumbe, E.; Vouillamoz, J.; Yacouba, H.; Guérin, R. 2017. Geophysical demonstration of the absence of correlation between lineaments and hydrogeologically useful fractures: Case study of the Sanon hard rock aquifer (central northern Burkina Faso). Journal of African Earth Sciences 129: 842-852. https://doi.org/10.1016/j.jafrearsci.2017.02.025

Wegner, W.; Wörner, G.; Harmon, R.S.; Jicha, B.R. 2011. Magmatic history and evolution of the Central American Land Bridge in Panama since Cretaceous times. GSA Bulletin 123: 703-724. https://doi.org/10.1130/B30109.1