Three-dimensional inverse modeling of EM-LIN data for the exploration of coastal sinkholes in Quintana Roo, Mexico

In the Yucatán Peninsula (YP), southern Mexico, cities and towns are settled on a platform of calcareous sedimentary sequence, where karst processes have formed numerous sinkholes, underground water conduits, and caverns. Anthropogenic activities there threaten the only source of freshwater supply, which is in a regional unconfined aquifer; there are no lakes or rivers on the surface. For the sustainable management of this resource in the YP, mathematical tools are needed in order to model groundwater. To determine the geometry of the aquifer, for example the positions of caves, sinkholes, and underground principal conduits, we modified a software to invert three-dimensional electromagnetic low-induction number (3-D EM-LIN) data for a set of profiles at arbitrary angles. In this study we used the EM-LIN geophysical method to explore the Chac-Mool sinkhole system in the state of Quintana Roo, Mexico. We performed inverse modeling in 3-D using the EM-34 instrument for vertical and horizontal magnetic dipoles. The 3-D inversion process yields models that enable us to correlate the path of the underground principal conduits with the subsurface electrical resistivity. In this work we show that inverse modeling of EM-LIN data can give us information about how close to surface the underground water conduits and the location of the boundary between fresh and salty water are.

1 Introduction

The main source of fresh water in the Yucatán Peninsula is a regional unconfined karst aquifer that is constituted by sedimentary limestones (Bauer-Gottwein et al., 2011). Karst aquifers are extremely vulnerable to contaminants because of their high permeability and the peculiar turbulent groundwater flow passing through karst conduits and caves (Worthington, 1999; Parise et al., 2015; Parise, 2019). Rapid population growth and coastal tourism in the state threaten the only source of freshwater supply in the peninsula.

In order to guarantee the sustainable use of this groundwater resource knowledge on the hydrogeological characteristics, such as geometry and position, of caverns and sinkholes and the depth of the freshwater/saltwater mixing zone (halocline), three-dimensional data inversion is needed.

Sinkholes are natural geological features connecting the land surface with underground karst terrains, and they are formed when rainwater dissolves limestone, creating underground voids (Coškun, 2012). Two main groups of sinkholes have been identified in the genetic classification (Williams, 2004; Gutierrez et al., 2008, 2014). The first group comprises solution sinkholes, which are formed by differential corrosion, lowering the ground surface where karst rocks are exposed. The second group comprises subsidence sinkholes, which result from both subsurface dissolution and downward gravitational movement.

Figure 1 Study area: Chac-Mool sinkhole in the state of Quintana Roo, Mexico.

In Quintana Roo many sinkholes, caverns, and underground water conduits have been reported by scuba divers, and the Quintana Roo Speleological Survey has produced an underground map of the Riviera Maya for tourism purposes. However, geophysical techniques have rarely been applied as noninvasive approaches to explore this area (Estrada-Medina et al., 2010; Gondwe et al., 2010; Beauer-Gottwein et al., 2011). Electrical resistivity tomography has proven effective for exploring karst areas (Ahmed and Carpenter, 2003; Chalikakis et al., 2011); however, in the Quintana Roo region the lack of soil on the hard limestone terrain has made placing electrodes a complicated and time-consuming task, raising expenses for data collection. New approaches to geophysical and coastal karst prospecting are therefore needed to develop and maintain sustainability plans in the Yucatán Peninsula (YP).

In this study we aim to explore a novel approach by using electromagnetic (EM) methods at low-induction numbers (LIN) and applying 3-D geophysical inverse modeling (Pérez-Flores et al., 2012) in order to set up a conceptual model of a sinkhole system and gain more knowledge on the geomorphology of the site. The methodology and results could be useful tools for the management of the Quintana Roo coastal zones, which is important for tourism and requires accurate information for prospect plans of development.

We did not find references on the use of EM-LIN in karst systems, but we found that the direct current (DC) and aero-TDEM (Time Domain Electromagnetic Method) were used for the Sian Ka'an Biosphere Reserve by Supper et al. (2009). These authors performed EM-34 measurements, but they did not do any further processing, like performing geophysical inversion.

1.1 Study area

This research was carried out in the Yucatán Peninsula (YP), an area largely dominated by karst landscape (Bauer-Gottwein et al., 2011). From the geological point of view, the YP is constituted by a sequence of calcareous sediments (Bonet and Butterlin, 1962) and is characterized by its flat landform (no topography) and the absence of surface rivers. A review of the YP karst aquifer is well described by Bauer-Gottwein et al. (2011), and an extended description of coastal cave development is given by Smart et al. (2006).

Our study area covers the Chac-Mool sinkhole and is 20 km south of Playa del Carmen in the state of Quintana Roo (approximately 20^∘30^′46.37^′′ N and 87^∘14^′49.32^′′ W) (Fig. 1). The area extends to 1 km² and is fully covered by dense vegetation. The ground presents high secondary porosity. Annual precipitation there is around 1200 mm, and topography is a flat surface with a slope of 9 m a.s.l. (above sea level) within 20 km from the shoreline (CNA, 2016). The hydraulic gradient in the southern part of Playa del Carmen was estimated at 58–130 mm km⁻¹ (Beddows, 2004). Due to its proximity to the coast (2 km), the study area is penetrated by seawater. Water intrusion is dependent on tides and rainfall (Beddows, 2004). Chac-Mool is a sinkhole complex where two underground water conduits presumably connect the Little Brother sinkhole and the Air Dome sinkhole. The underground river pathways in some sections have been documented on maps made by scuba divers (Quintana Roo Speleological Survey), but other sections and vertical components remain unknown. The entire rock matrix is possibly saturated with fresh and brackish water in pores and small conduits. The apparent conductivity is high because it averages the matrix conductivity (low value) with the seawater conductivity (high value).

Figure 2 EM survey on the Chac-Mool sinkhole. The numbered profiles cross the hidden underground water conduits. White lines mark the sinkholes.

1.2 Electromagnetic survey

In September 2015, we carried out a field survey over the study area. We obtained seven profiles with the EM-34 (Geonics) instrument, which operates within the LIN domain as described in McNeill (1980). The main reason for using the EM-34 is that it can accurately obtain data in a more easy and faster way in terrains with no soil, expediting field work in hard terrains.

The basic principle consists in the transmission of an alternating current of constant frequency (f) through a coil, which generates a primary electromagnetic field (H_p) that induces electrical currents in the conductive bodies embedded in the subsoil (following Faraday's law). A secondary electromagnetic field in the ground (H_s) is then generated by the conductive bodies. These two fields differ in amplitude and phase, and they are detected by a coil (receiver) that is separated by a distance s (m) from the transmitter. The induction number, N, is defined as the quotient between s (m) and the skin depth δ (m): N=s (m)∕δ (m). At LIN (N<1) the imaginary part of H_s∕H_p is a straight line for which the slope is the conductivity of a homogeneous half-space. Because the ground is not homogenous, we say we get an apparent conductivity: ${\mathit{\sigma }}_{\mathrm{a}}=(\mathrm{4}/\mathit{\omega }{\mathit{\mu }}_{\mathrm{0}}{S}^{\mathrm{2}})(H_{\mathrm{s}}/H_{\mathrm{p}})$.

Both loops (source and receiver) are commonly used on the same plane (coplanar). We have two possible arrays, one when both loops are parallel to the Earth's surface (vertical magnetic dipoles or VMD) and the other when both loops are perpendicular to the Earth's surface (horizontal magnetic dipoles or HMD). The separation between loops can be extended to 10, 20, and 40 m in both arrays. For this study, measurements were made along six lines (Fig. 2), and the observation points were spaced every 5 m. Because vegetation in the jungle was dense, we were unable to locate profiles anywhere, and so we took the paths around the Chac-Mool, Little Brother, and Air Dome sinkholes. The distribution of the six profiles covers approximately a rectangular area. Therefore, we performed a 3-D inversion in order to follow the three-dimensional complexity of the sinkholes. For the 3-D inverse modeling we followed the method by Pérez-Flores et al. (2012), but the algorithm they used was designed for profiles that were measured in parallel (0^∘) or perpendicular (90^∘) positions with respect to the other profiles. Later on, we show how we modified the equations for arbitrary angle profiles. The length of the six profiles (1 to 6) varies between 50 and 140 m (Fig. 2).

1.3 Inverse modeling

EM data (apparent conductivity, σ_a) were assumed to be approximately the weighed average of the subsurface electrical conductivity distribution, as described by Pérez-Flores et al. (2012). We associated the apparent conductivity (σ_a) with the true subsurface conductivity (σ) by means of a weighting function (that is, the Green function and electric-field product) using the integral equation formulated by Pérez-Flores et al. (2001):

$$\begin{array}{}\text{(1)}& {\mathit{\sigma }}_{\mathrm{a}}\left({\mathbf{r}}_{\mathrm{2}},{\mathbf{r}}_{\mathrm{1}}\right)\cong -{\displaystyle \frac{\mathrm{16}\mathit{\pi }s}{\mathit{\omega }{\mathit{\mu }}_{\mathrm{0}}m}}\underset{v}{\int }\mathbf{G}\left({\mathbf{r}}_{\mathrm{2}},\mathbf{r}\right)\cdot \mathbf{E}\left(\mathbf{r},{\mathbf{r}}_{\mathrm{1}}\right)\mathit{\sigma }\left(\mathbf{r}\right)\mathrm{d}v,\end{array}$$

where r₁ and r₂ are the positions of the source and the receiver, G is the Green function for a homogeneous medium, and E is the electric field for a homogeneous half-space. Equation (1) is an approximation for the low-conductivity contrasts, and it is very useful for an inversion, where G, E, and σ_a are known and σ(r) is unknown.

For the inversion we had to consider how the magnetic dipoles were used. We obtained the vertical and horizontal magnetic dipole (VMD and HMD, respectively) arrays as described by Pérez-Flores et al. (2012). The method we are using was developed and explained in Pérez-Flores et al. (2012). In the following paragraphs we will make a simple modification to the equations already published in order to accept arbitrary profile azimuths. The integral equation for VMD is

$$\begin{array}{}\text{(2)}& {\mathit{\sigma }}_{\mathrm{a},z}\left({\mathbf{r}}_{\mathrm{1}},{\mathbf{r}}_{\mathrm{2}}\right)\cong -{\displaystyle \frac{\mathrm{16}\mathit{\pi }s}{\mathit{\omega }{\mathit{\mu }}_{\mathrm{0}}{m}_z}}\underset{v}{\int }{\mathbf{G}}_{H_z}\left(\mathbf{r},{\mathbf{r}}_{\mathrm{2}}\right)\cdot {\mathbf{E}}_{H_z}\left(\mathbf{r},{\mathbf{r}}_{\mathrm{1}}\right)\mathit{\sigma }\left(\mathbf{r}\right)\mathrm{d}v.\end{array}$$

For HMD the integral equation in the y direction is given by

$$\begin{array}{}\text{(3)}& {\mathit{\sigma }}_{\mathrm{a},y}\left({\mathbf{r}}_{\mathrm{1}},{\mathbf{r}}_{\mathrm{2}}\right)\cong -{\displaystyle \frac{\mathrm{16}\mathit{\pi }s}{\mathit{\omega }{\mathit{\mu }}_{\mathrm{0}}{m}_y}}\underset{v}{\int }{\mathbf{G}}_{H_y}\left(\mathbf{r},{\mathbf{r}}_{\mathrm{2}}\right)\cdot {\mathbf{E}}_{H_y}\left(\mathbf{r},{\mathbf{r}}_{\mathrm{1}}\right)\mathit{\sigma }\left(\mathbf{r}\right)\mathrm{d}v.\end{array}$$

HMD in the x direction is given by

$$\begin{array}{}\text{(4)}& {\mathit{\sigma }}_{\mathrm{a},x}\left({\mathbf{r}}_{\mathrm{1}},{\mathbf{r}}_{\mathrm{2}}\right)\cong -{\displaystyle \frac{\mathrm{16}\mathit{\pi }s}{\mathit{\omega }{\mathit{\mu }}_{\mathrm{0}}{m}_x}}\underset{v}{\int }{\mathbf{G}}_{H_x}\left(\mathbf{r},{\mathbf{r}}_{\mathrm{2}}\right)\cdot {\mathbf{E}}_{H_x}\left(\mathbf{r},{\mathbf{r}}_{\mathrm{1}}\right)\mathit{\sigma }\left(\mathbf{r}\right)\mathrm{d}v.\end{array}$$

The expressions for ${\mathbf{G}}_{H_z}$, ${\mathbf{E}}_{H_z}$, ${\mathbf{G}}_{H_y}$, ${\mathbf{E}}_{H_y}$, ${\mathbf{G}}_{H_x}$ and ${\mathbf{E}}_{H_x}$ can be consulted in Pérez-Flores et al. (2012). VMD profiles can run at any angle (Eq. 2), but HMD profiles run only in either the y direction (90^∘; Eq. 3) or x direction (0^∘; Eq. 4). Arbitrary direction profiles like those observed around the Chac-Mool sinkhole (Fig. 3) constituted a problem. So, we had to modify Eqs. (4) and (5) in order to accept the arbitrary angle profiles.

Figure 3 Profiles crossing underground water conduits in the sinkhole area (numbered lines). The white rectangle is the 3-D modeled area. White lines mark the sinkhole boundaries. Dark blue lines are the suggested underground water conduits paths.

Using a simple rotation for E and G in terms of their vector components, the y direction for HMD is

$$\begin{array}{}\text{(5)}& {\mathbf{G}}_{H_y}\left(\mathbf{r},{\mathbf{r}}_{\mathrm{2}}\right)=d\widehat{i}+e\widehat{j},{\mathbf{E}}_{H_y}\left(\mathbf{r},{\mathbf{r}}_{\mathrm{1}}\right)=a\widehat{i}+b\widehat{j},\end{array}$$

and the x direction is

$$\begin{array}{}\text{(6)}& {\mathbf{G}}_{H_x}\left(\mathbf{r},{\mathbf{r}}_{\mathrm{2}}\right)=e\widehat{i}+f\widehat{j},{\mathbf{E}}_{H_x}\left(\mathbf{r},{\mathbf{r}}_{\mathrm{1}}\right)=b\widehat{i}+c\widehat{j}.\end{array}$$

When we rotate Eq. (3) 90^∘, it becomes Eq. (4). So, we can find E and G in terms of their rotated components:

$$\begin{array}{}\text{(7)}& {\displaystyle }\begin{array}{rl}& \left(\begin{array}{c}{E}_x\\ E_y\end{array}\right)=\left(\begin{array}{ccc}cos\mathit{\theta }& \mathrm{sen}\mathit{\theta }& \mathrm{0}\\ \mathrm{0}& cos\mathit{\theta }& \mathrm{sen}\mathit{\theta }\end{array}\right)\left(\begin{array}{c}a\\ b\\ c\end{array}\right),\\ & \left(\begin{array}{c}{G}_x\\ G_y\end{array}\right)=\left(\begin{array}{ccc}cos\mathit{\theta }& \mathrm{sen}\mathit{\theta }& \mathrm{0}\\ \mathrm{0}& cos\mathit{\theta }& \mathrm{sen}\mathit{\theta }\end{array}\right)\left(\begin{array}{c}d\\ e\\ f\end{array}\right).\end{array}\end{array}$$

If an HMD profile runs at 0^∘, (E_x, E_y) becomes ${\mathbf{E}}_{H_y}$ from Eq. (3). If the profile runs at 90^∘, (E_x, E_y) becomes ${\mathbf{E}}_{H_x}$ from Eq. (4).

Thus, for an arbitrary angle profile, Eqs. (3) and (4) become a single equation,

$$\begin{array}{}\text{(8)}& {\displaystyle }\begin{array}{rl}{\mathit{\sigma }}_{\mathrm{a}}\left({\mathbf{r}}_{\mathrm{1}},{\mathbf{r}}_{\mathrm{2}}\right)& =-{\displaystyle \frac{\mathrm{16}\mathit{\pi }s}{\mathit{\omega }{\mathit{\mu }}_{\mathrm{0}}m}}\int \left[G_x\left(\mathbf{r},{\mathbf{r}}_{\mathrm{2}}\right)E_x\left(\mathbf{r},{\mathbf{r}}_{\mathrm{1}}\right)\right.\\ & \left.+G_y\left(\mathbf{r},{\mathbf{r}}_{\mathrm{2}}\right)E_y\left(\mathbf{r},{\mathbf{r}}_{\mathrm{1}}\right)\right]\mathit{\sigma }\left(\mathbf{r}\right)\mathrm{d}v.\end{array}\end{array}$$

For terms a–e and f see Pérez-Flores et al. (2012).

For the 3-D inversion, we used Eq. (2) for the HMD profiles and Eq. (2) for the VMD profiles. We used 10, 20, and 40 m as the source–receiver separations for VMD and HMD in every profile. We pooled all data sets and performed a joint inversion to obtain a single 3-D conductivity model. We simulated the heterogeneous half-space as a conglomerate of rectangular prisms. We assumed that conductivity in every single prism was constant, however unknown. Equations (2) and (8) can be written as a linear equation system and in a matrix fashion:

$$\begin{array}{}\text{(9)}& {\mathit{\sigma }}_{\mathrm{a}}=\mathbf{W}\mathit{\sigma },\end{array}$$

where σ_a represents the column vector of apparent conductivities, matrix W contains the weights or products of the Green function and electric field and is partitioned for VMD and HMD, and σ represents the column vector of the real conductivities (unknowns). We used quadratic programming to minimize the following objective function, U:

$$\begin{array}{}\text{(10)}& {\displaystyle }\begin{array}{rl}& U\left(\mathit{\sigma }\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left|\left|{\mathit{\sigma }}_{\mathrm{a}}-\mathbf{W}\mathit{\sigma }\right|\right|}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\beta }{\left|\left|\mathbf{D}\mathit{\sigma }\right|\right|}^{\mathrm{2}}\\ & {\mathit{\sigma }}_{\mathrm{lower}}<\mathit{\sigma }<{\mathit{\sigma }}_{\mathrm{upper}}.\end{array}\end{array}$$

Matrix D represents the first-order spatial derivatives of the contiguous prism conductivities. Parameter β controls the smoothness of the 3-D conductivity model; when it was low, we obtained a rough 3-D model. The first term fits the apparent conductivity data taken at the field. The second term in Eq. (10) contains the spatial derivatives of the conductivity in (x, y, z) direction. The smoothness parameter controls the magnitude of the second term. If zero, only the data were fit and the model used to be very rough; if very large, the model converged into a homogenous half-space. We transformed the Hessian to achieve diagonal unity. This way the smoothness parameter varies in a very narrow window. We tested the values 0.1, 0.01, and 0.001. The 0.1 value yields a smooth model and the 0.001 value a rough model. We began with a smooth value that gave the simplest but the most probable model (according to the Occam's Razor principle), and we lowered the parameter to recover more structure; however, after a certain point the structure turned unreal from the geological point of view. The idea was to recover most of the structure while keeping the simplest and most probable model.

2 Resistivity cross sections on the 3-D model

For the 3-D inverse modeling we used an (x, y, z) grid of prisms, assuming constant conductivity in every prism. We performed the inverse modeling choosing $\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{2.5}$ m in the (x, y) directions because the EM measurements were taken every 5 m; the variable discretization of Δz was chosen to be (0, 2, 5, 8, 12, 18, 25, 35, and 50 m), and β=0.01 was the smoothness factor.

Conductivity is the unknown, but we prefer to show the resistivity (the inverse of conductivity) results. In Fig. 4 we present the 3-D resistivity model after the inversion of whole sets of data. In that figure we present the interpolated resistivity cross sections under the six profiles. Blue indicates resistivity areas and red low resistivity. There are spaces between profiles with no data. The 3-D model for those areas is not so reliable. Therefore, as a first approach, we show the model for the areas for which we had data. There is very good coherence where the model crosses. Figure 4 shows irregular paths for the two underground water conduits, according to the map from the divers (x, y, z). Water table depth in the open sinkholes is 7 m. The water conduits follow very intricate paths. We think that there are narrower river branches that have not yet been mapped by the divers. Interestingly, some paths were marked below the resistivity areas. The upper water level of the subterranean river is probably far from the surface, making the rock mass or resistive mass more structurally stable, or maybe those resistive bodies are air-filled caves over the water table. By resistive mass (RM) we refer to the dry limestone rock between the surface and the ceiling of the cave and the air-filling the cave. We can idealize a typical cave in this area (near the coast), vertically consisting of dry limestone (resistive), an air space (resistive), followed by fresh water (lower resistivity), the halocline (mixing of fresh and salty water), and, at the bottom, salty water (lowest resistivity) surrounded by saturated limestones as bedrock. We cannot distinguish in the RM what is dry limestone and what is air-filling the cave.

Figure 4 Three-dimensional resistivity model for the Chac-Mool sinkhole complex. Here, we show only the distribution of the cross sections where the profiles were run. The red and black irregular lines represent the underground water conduits.

In Fig. 5 we show the six cross sections obtained with the 3-D resistivity model. Cross section (a) corresponds to the profile-1 model, cross section (b) to the profile-2 model, and so on. Every profile is indicated with a white circle, which pinpoints the interpolated (x, y, z) hidden-water conduits. The (x, y, z) locations were obtained from the mapping made by scuba divers. We delineated the inferred cave section with a rectangle, because we could not see details. We assumed the saturated limestone was bedrock, because dry limestone resistivity was larger than 1000 Ω m. In the 3-D-model cross sections, the bedrock looks green everywhere (160–170 Ω m). Only some small sections were blue (1000 Ω m).

Figure 5 Cross sections of the 3-D resistivity model for profiles 1 to 6. Resistivity units are base 10 logarithm. Blue color indicates more resistive areas and red the least resistive areas. Blue numbers indicate the other profile crossings. White circles pinpoint the areas where scuba divers have mapped the underground water conduits. Red circles show the position of an underground river inferred from the model. The square polygon is a broad suggestion of the river tunnels.

From the six resistivity cross sections, we can see that most of the river crossings show a green color over them. This means that the subterranean water conduits are probably close to the surface and the thickness of the RM is therefore thin, meaning RMs in those areas are more vulnerable to sinking, though we did not find evidence of subduction or fracturing on the surface. The cross section for profile 1 (Fig. 5a) shows three crosses: one at x=18 m showing a thin RM and the other two showing a thicker RM. Profile 2 (Fig. 5b) shows a green color, meaning thinner RMs. Profile 3 (Fig. 5c) shows one river crossing that is shallow and another deeper one. We clearly detected a shallower subterranean river (green color) using the EM-LIN equipment, but it is not clear how much deeper it is. We must remember that the white circles are interpolations taken from the diver's map. The deeper river crossing coincides with the location of a large resistivity mass between zero and 20 m; this means that divers had to dive below this resistivity mass (1000 Ω m). Profile 4 (Fig. 5d) shows three crossings with green color. Profile 5 (Fig. 5e) shows three crossings, two are deep (between z=20 and z=30 m) and one is shallow (z=15 m). The deeper crossings are consistent with the reported diving depth and the thicker RM shown by the large resistivity mass. However, at x=25 m the river seems to be 10 m deeper, possibly because of the presence of a huge hard rock (very resistive). Profile 6 (Fig. 5f) shows a shallow river and a deeper one. Resistivities are consistent with the position of the river.

We know that divers swam throughout subterranean water conduits. In Fig. 5 we broadly suggest the location of the river crossings (rectangular polygon). Given the color descriptions in Fig. 5, we can say that blue is an indication of dry limestone RMs or dry limestone and air-filled caves at the top of the water conduit or close to the surface. The green color is so widespread that it surely indicates fresh water (50 to 70 Ω m). Also, the resistivity cross section shows a green color where the subterranean water conduits seem to be shallow. We expected to see a narrow blue coloration and green color over those shallow water conduits, but we did not, because the narrowest source–receiver separation at the EM-34 was 10 m (too large to see surface details). In some way the estimated true conductivity is still an average. Maybe if we use a shorter separation, we could see a thinner blue color for the RM and then a green color for the fresh water. The transition from green to red (yellow) could be the transition from fresh water to salty water. We expect fresh water at the top and salty water at the bottom because of the density.

We drew the river section to emphasize that the resolution of the EM-34 instrument is not good enough to sharply isolate the water conduits from the bedrock. A possible explanation is that the upper sections of unaltered bedrock (limestone) are partially saturated with fresh water (because of the 50–70 Ω m values) and the deeper sections are saturated with salty water (because of the 6–10 Ω m values). So, there are no large horizontal resistivity differences between the river location and the bedrock. It is almost certain that the permeability of the bedrock is as high as the permeability of the limestone at the surface. When it rains, water quickly disappears. Aerial electromagnetics (flying 30 to 50 m over the surface) would yield an even lower resolution (Supper et al., 2009).

In profile 4 (Fig. 5d) there is a green color section close to x=70 m (red square). It is possible that there is a shallow subterranean river close to the surface that has not yet been mapped by the divers.

Isometrics of the 3-D resistivity model

The Chac-Mool sinkhole system is a complex of three small sinkholes (Air Dome, Little Brother, and Chac-Mool). According to divers, there are two underground water conduits. Their vertical variations may cause thinning of the limestone RMs and therefore sinking. According to the cross section in Fig. 6, the EM-LIN equipment cannot sharply distinguish between the subterranean river tunnels and the bedrock, maybe because there is not enough change in resistivity. This means that limestone bedrocks are partially saturated with water and therefore under a process of chemical dissolution. The isometric view of the 3-D resistivity model (Fig. 6) shows the spatial distribution of the three sinkholes in the system, the two proposed water conduits and their paths, and the location of the five EM-LIN data profiles.

Figure 6 Isometric representation of the 3-D resistivity model. Straight lines represent the EM profiles. (a) Blue isosurface representing the bottom topography of the dry limestones. (b) Orange isosurface representing the area where fresh and salty waters meet.

The blue and orange surfaces are equal-resistivity surfaces in the 3-D model. The blue surface shows the contact between dry limestone (resistive) and the fresh water (∼80 Ω m). The resistive layer may contain unaltered limestone and/or air-filled caves. It is very interesting that this layer outcrops where underground water conduits are very close to the surface, maybe because the shortest source–receiver distance (10 m) is larger than the RM thickness. This surface does not show where the sinkholes are, because of the lack of data. We did not manipulate the 3-D model in order to force outcrops of areas with sinkholes. The orange surface represents the contact between the fresh water and the salty water (Halocline).

3 Conclusions

In this research we studied the Chac-Mool sinkhole complex by EM methods at LIN. These methods consist of a source loop and a receiver loop operating in two coplanar arrays VMD and HMD. These two arrays or polarizations view inside the Earth in two different ways. We used both arrays to perform a joint inversion and to obtain a single three-dimensional (3-D) resistivity model. Equations had already been published for a mesh of perpendicular and parallel profiles but not for arbitrary angle profiles (Pérez-Flores et al., 2012). In this research the profiles were taken inside the jungle and we took the advantage of man-made paths; however, these paths were located at arbitrary angles. We modified the existing equations and obtained a more general set of equations.

The 3-D inversion of both VDM and HDM arrays led to a single 3-D resistivity model. The cross sections of this 3-D model show the points where the underground water conduits cross. The areas where the underground water conduits are close to the surface could represent hazard zones because of the possibility of RMs collapsing. We also observed the distribution of fresh and salty waters and the areas where they meet or the transition surface (halocline). Our observations indicate that water conduits might run along tunnels, but the resistivity of those tunnels does not differ sharply from the resistivity of the bedrock, meaning that bedrock could be saturated with water (fresh and salty depending on depth). The isometric view shows that the resistivity isosurface corresponds with the bottom topography of the underground RM. At the center of the area of study this RM seems to be very thick, indicating that this area is safe from sinking. This isometric view also shows the contact between fresh and salty water.

The EM-LIN technique is a fast, efficient, and inexpensive procedure for explorations over hard-rock sinkhole areas. It allows us to obtain the geometry of the underground water conduits and the distribution of fresh and salty water.

4 Discussion

EM-LIN methods can be applied for detecting natural caves under roads or caused by shallow mining or for determining resistivity changes caused by landsliding. This geophysical technique is ideal to detect resistivity changes underground, and it has the advantage that no-electrode grounding is required, making it faster and cheaper than DC resistivity. The equations contained in this paper can be easily computed for the determination of 3-D resistive/conductive bodies buried underground. The disadvantage is that commercial equipments were produced for few source–receiver separations limiting the resolution and the penetration depth. The equipment used for acquiring the data of this paper has no source–receiver separations lower than 10 m, which makes difficult to accurately resolve the very shallow thickness of the dry limestones.

We also point out that this EM method works better when the conductivity of the buried target is very different from the bedrock conductivity. In this paper we observed that the subterranean conduit conductivity was not very different from the bedrock conductivity as we expected. These results encourage us to make an accurate EM-LIN modeling in order to find an explanation of the low resistivity contrasts observed.