- Email: [email protected]
3D magnetotelluric characterization of the geothermal anomaly in the Llucmajor aquifer system (Majorca, Spain)
Journal of Applied Geophysics 68 (2009) 479–488
Contents lists available at ScienceDirect
Journal of Applied Geophysics j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / j a p p g e o
3D magnetotelluric characterization of the geothermal anomaly in the Llucmajor aquifer system (Majorca, Spain) C. Arango ⁎, A. Marcuello, J. Ledo, P. Queralt Departament de Geodinàmica i Geofísica, Facultat de Geología, Universitat de Barcelona, C/Martí i Franquès s/n, 08028 Barcelona, Spain
A R T I C L E
I N F O
Article history: Received 15 April 2008 Accepted 23 May 2008 Keywords: Llucmajor aquifer Magnetotellurics Geothermal exploration Geoelectrical structure 3D modeling
A B S T R A C T In the Llucmajor aquifer system (Majorca Island, Spain) some geothermal evidences have appeared. This phenomenon is not isolated to Majorca and it is present in other areas, where it can be associated with structural conditions, especially to the extensional event suffered by the island after the Alpine Orogeny. However, the origin of this anomaly in Llucmajor is not well known, and there is no surface geological evidence of these structural conditions. With the aim of delineating the geoelectrical structure of the zone and identifying the geological structure that allows the presence of this anomaly, an audiomagnetotelluric (AMT) survey was carried out. The AMT data was processed using a Wavelet Transform-based scheme. Dimensionality analysis indicates that the geoelectrical structure is mainly 3D. The 3D model was obtained by trial and error forward modeling, taking accounting of the responses from the determinant of the impedance tensor. The model shows a vertical resistivity distribution with three horizons associated with different units: on the top, a shallow high resistive media related to an unconfined shallow aquifer; in the middle, a conductive layer related to the aquitard, and below it, another resistive media related to the confined deeper aquifer. The intermediate horizon shows a sudden thinning beneath the thermal anomalous zone that can be identified as a weakness zone (fault or fracture) connecting both aquifers. An exploratory well was drilled after the AMT survey and reached almost 700 m in depth. This allowed correlating the resistivity distribution of the 3D model with data logging and lithology obtained from the well, showing a proper agreement between them. © 2008 Published by Elsevier B.V.
1. Introduction The Llucmajor aquifer system is located in the southern part of Majorca Island (Spain). In this area water is extracted from the unconfined aquifer of the system, and used for agricultural purposes. From the exploitation wells the inferred thickness of the aquifer is around 150 m. A geothermal anomaly has been detected in this area, manifested by anomalous high water temperatures in some wells aligned along a NS–EW direction, and reaching a maximum temperature of 50 °C (Fig. 1). Because of the lack of volcanic sources in this region, these high values must be associated to vertical movements of water. Taking the reported thermal gradient in the zone as 46 ± 8 mK/m (Fernàndez and Cabal, 1992), a minimum depth of 700 m is needed to obtain such a temperature. Understanding the origin of this hot water, as well as its possible economical exploitation was the aim of a research program conducted by the IGME (Spanish Geological Survey) and supported by the Balearic Government. To explain this behavior a simple conceptual model with two aquifer units was proposed: one unconfined aquifer
⁎ Corresponding author. Departamento de Geomagnetismo y Exploración, Instituto de Geofísica, UNAM, 04510 México, D.F., Mexico. E-mail address: claudiar@geofisica.unam.mx (C. Arango). 0926-9851/$ – see front matter © 2008 Published by Elsevier B.V. doi:10.1016/j.jappgeo.2008.05.006
on top, and another confined aquifer to supply the hot water. The aquitard between both aquifers must present a discontinuity (fault or fracture) to favor vertical fluid movements. This area is composed of Quaternary silts, and there is no surface geomorphologic evidence of structures that allows the rising of deep groundwater. In the studied zone, two geophysical surveys were reported. The first one employed vertical electrical soundings (IGME, 1985), but a high resistive shallower layer did not allow good penetration of the electrical current, and limited the range of the investigation depth to the upper aquifer. The second geophysical survey is related to the seismic reflection campaign carried out in Majorca between 1984 and 1990 (Benedicto et al., 1993). In this campaign, one of the profiles went through Llucmajor zone, but did not unfortunately have enough resolution to study the origin of the thermal anomaly. Taking this information into account, the magnetotelluric exploration was considered the proper technique for generating a model to explain the origin of the geothermal anomaly. Among the electromagnetic methods, the audiomagnetotellurics (AMT) method is appropriated for deep groundwater and geothermal studies (Tezkan, 1999; Meju, 2002) given the depth of its penetration and its tensorial nature. Over recent years, different studies of deep aquifers (Ritz et al., 1997; Wannamaker, 1997; Meju et al., 1999; Krivochieva and Chouteau, 2003) and geothermal problems (Ogawa et al., 1998; Volpi et al., 2003) have been reported.
480
C. Arango et al. / Journal of Applied Geophysics 68 (2008) 479–488
Fig. 1. Temperature distribution in producing water wells of the Llucmajor zone.
New methodological aspects have recently been developed in magnetotellurics, in particular: the use of Wavelet Transform for time series processing (Zhang and Paulson, 1997; Trad and Travassos, 2000; Arango, 2005); and the use of the invariants to analyze the geoelectrical dimensionality (Weaver et al., 2000; Martí et al., 2004). This work incorporates them into the study of this aquifer. The main objectives are: 1) to determine the 3D geometry of the aquifer system and 2) to identify the zone where the aquitard could be fractured. We first describe the geology of the area. We then describe the main aspects of the magnetotelluric method, emphasizing the methodological elements such as geoelectrical dimensionality and the use of the determinant of the impedance tensor in the interpretation. Later, the model is correlated with data from the exploratory well, drilled following the MT fieldwork. 2. Geological setting Majorca is the biggest island of the Balearic archipelago and is considered the NE continuation of the Betic chain, located in southern Spain. The main structures of the island are identified by three mountain ranges with a NE–SW direction (Tramuntana, Randa and Levante systems) and separated by two main depressions (Gelabert, 1998). Its current configuration is conditioned by Neogene distension, which formed small basins. Three of them (Palma, Inca and Sa Pobla) are aligned in the western area, limited by NE–SW faults and separated from each other by NW–SE structures. The fourth (Llucmajor–Campos), in the eastern area, shows the same fault pattern (Fig. 2). The Llucmajor platform is part of one of these depressions. It is mainly composed of Quaternary materials (red silts, conglomerates and aeolianites), and Upper Miocene materials (limestones). Reefal calcarenites and discontinuous Quaternary silts crop out in the northern zone (Pomar, 1991). Several regional geophysical surveys have been carried out in Majorca using gravity, magnetic, seismic and MT methods (Benedicto et al., 1993; Ayala et al., 1994; Pous et al., 1995). They confirm that the main regional structures observed on surface continue beneath the island.
3. MT data survey 3.1. The MT method The magnetotelluric method is a geophysical technique that determines ground electrical resistivity distribution from the simultaneous measurement of the fluctuations of the natural electromagnetic field. Under the plane wave assumption, the linear relationship between the horizontal components of electric and magnetic fields at a given frequency, f, are expressed as follow: Eð f Þ = Z ð f ÞHð f ÞZ
Ex Ex
=
Zxx Zyx
Zxy Zyy
Hx Hx
ð1Þ
where Z is the impedance tensor. The impedance is a complex magnitude, from which it is customary to define the apparent resistivity and the phase for each component of the tensor as: ρaij ð f Þ =
1 2 jZ ð f Þj 2πf μ ij
! − 1 Im Zij ð f Þ uij ð f Þ = arg Zij ð f Þ = tan Re Zij ð f Þ
ð2aÞ
ð2bÞ
where µ is the magnetic permeability, and i, j denote any horizontal component. In this work we have also used the apparent resistivity and phase calculated from the determinant of the impedance tensor, introduced by Berdichevsky and Dmitriev (1976) as effective impedance: ρa ð f Þ =
uð f Þ =
1 jdet Z ð f Þj 2πf μ 1 argðdet Z ð f ÞÞ 2
ð3aÞ
ð3bÞ
These responses have interesting properties for modeling: they are rotationally invariant, the phase is not affected by galvanic distortions,
C. Arango et al. / Journal of Applied Geophysics 68 (2009) 479–488
Fig. 2. Geological setting of the Majorca Island (Spain) modified from Gelabert (1998).
Fig. 3. Location of AMT and MT sites and interpreted profiles.
481
482
C. Arango et al. / Journal of Applied Geophysics 68 (2008) 479–488
and both the apparent resistivity and the phase contain information from the four complex elements of the impedance tensor. Its use simplifies the data employed in the 3D modeling and integrates information of all components of the impedance tensor. 3.2. Magnetotelluric data The fieldwork was designed following two main perpendicular profiles, the inner section being filled out using a grid to satisfactorily cover the area of interest (Fig. 3). The average distance between sites was 250 m. Data acquisition had been carried out in summer 2002 at 42 different sites with two different Metronix systems. 36 measurements were taken with the GSM06 system, whose frequency ranges from 8 · 103 to 2.5 · 10− 3 Hz.10 measurements were taken with the MMS03E system, whose frequency ranges from 2.5 · 102 to 2.5 · 10− 3 Hz. Four sites were recorded with both instruments to test the agreement between results. The magnetotelluric horizontal components had been measured along the N–S and E–W directions. The recording time of most of the AMT stations were only a few hours, so only responses for higher frequencies (N0.1–1 Hz) were computed, and the MT stations covered lower frequencies. Although simultaneous recordings between both systems were performed, no remote reference was possible, since synchronization was unfeasible.
The vertical magnetic component presented electronic problems, producing very scattered tipper responses, and was not considered. For an estimation of impedance from the time series we employed the Wavelet Transform algorithm as it is described in Zhang and Paulson (1997), and implemented and modified by Arango (2005). This technique allows for a detailed analysis in the time-frequency domain, as well as a description of frequency content along time, which is not possible using the classic method based on the Fourier analysis. The mother wavelet used in this work to compute the Wavelet Transform of the electromagnetic fields was the function proposed by Morlet et al. (1982). It was chosen for its capacity to concentrate energy in a narrow frequency band, as well as its simple lineal relationship between frequency and scale. The most important feature of this processing method is the selection of intervals with high energy content for each scale (frequency) as suggested in Zhang and Paulson (1997). These are used to calculate the impedance tensor elements to ensure a better signal in the dead band. After several attempts and for this work, the threshold level between high and low energy content was taken as the mean energy. This offers a proper tradeoff for our data quality, as shown in Fig. 4, where a comparison between the Wavelet Transform algorithm and a standard process is presented. Both processes present similar results, although the Wavelet Transform-based appears as less scattered.
Fig. 4. Apparent resistivity and phase curves computed using standard scheme(Mapros-software, Friedrichs, 2001) and the Wavelet Transform-based scheme.
C. Arango et al. / Journal of Applied Geophysics 68 (2009) 479–488
483
Fig. 5. (a) Directionality analysis for different frequency intervals using the algorithm by Martí et al. (2004). (b) Strike direction computed between 1 and 100 Hz.
The data is restricted to frequencies higher than 0.1 Hz, focusing on the shallow part of the geoelectrical structure. The sea is located 10– 15 km away. This is one to two times the skin depth, if both the frequency range, and a non-resistive basement in this region (Pous et al., 1995), are considered. Therefore there is no coast effect in the studied zone for the chosen range of frequencies.
An important step prior to interpretation is dimensionality analysis, which allows identification of the main directions of geoelectrical structures. The analysis of dimensionality used in this work is based on the invariants of the magnetotelluric impedance tensor, as described in Weaver et al. (2000), and implemented following the algorithm of Martí et al. (2004). This analysis showed
484
C. Arango et al. / Journal of Applied Geophysics 68 (2008) 479–488
that the behavior of the geoelectrical structure between 0.1 and 1000 Hz is three-dimensional, and supported this description with a 3D model (Fig. 5a). It can be also seen in this figure that several sites exhibit 2D behavior between 1 and 100 Hz, and this is employed to obtain a 2D model, used as preparatory model for the further 3D interpretation. According to Ledo (2005), in some cases it is possible to extract useful information from 3D data using a 2D algorithm, but this approach must be taken carefully however, since it depends on both geoelectrical configuration and the particular features of the media. Furthermore, a determination of the geoelectrical strike was also done to all sites using the Strike program (McNeice and Jones, 2001) based on the Groom–Bailey decomposition (Groom and Bailey, 1989). In this frequency interval, this analysis indicated a strike of N50°E (Fig. 5b). Static shift was corrected, assuming that there were no important variations on the very shallow structures, as suggested from the geological information. It was thus expected that the near surface resistivity would be essentially the same at all the sites. This assumption is also supported by the results of the vertical electrical soundings in this area (IGME, 1985). Most of the apparent resistivities coincided with the same value at high frequencies, taken as the average, which is taken as the near surface resistivity. Three perpendicular profiles along the NW–SE direction and identified as profiles 1, 2 and 3 (Fig. 3) were considered. The 2D resistivity models were computed using the REBOCC inversion algorithm (Siripurnvaraporn and Egbert, 2000), modified by Pedersen and Engels (2005) to use the determinant of the impedance tensor. This type of data ensures a better data fitting when a strike is not fully determined, as is the situation with 2D interpretation of 3D structures. The final 2D models are shown in Fig. 6 and the integrated rms value for each profile was 2.97, 3.14 and 8.61 respectively, with a 5% error floor in the impedance components. The higher misfit in profile 3 is probably caused by bad data quality, and one of the sites presents a strong 3D behavior (Fig. 5a).
In the 2D models we can identify three layers with different electrical properties along all profiles. The intermediate horizon has a low resistivity value, unlike the other two layers, and its thickness varies suddenly beneath the thermal zone. These results were used for creating a 3D model grid, as well as for building the starting 3D model. We assigned resistivity values obtained from 2D models to create a 3D model grid with information based on previous results. 3.3. 3D model The acquisition pattern designed allowed the construction of a 3D model to characterize electrical resistivity distribution. Such a model consisted of 35 × 35 cells in NS–EW between 115 and 660 m and 20 layers of variable thickness between 20 and 250 m. This configuration ensured that grid is refined over zones containing more information. To account for the boundary conditions, the grid has been laterally extended twice the size of the study area. The procedure for obtaining the model was established through a trial and error strategy with a varying resistivity value of the cells using the determinant of the impedance tensor as a model response for data comparison. Once such model is found, it is tested for adjusting the responses of the antidiagonal elements of the impedance tensor. Forward modeling was computed using the algorithm developed by Mackie et al. (1993), and modified by Booker and Mackie (personal communication, 1999). The data frequencies ranged from 1000 Hz to 0.1 Hz. The 3D geoelectrical model that we obtained shows the main subsurface structures up to 600 m in depth (Fig. 7). The fit between data and model responses was satisfactory, and the rms was 6.07 with an error floor of 5% in the impedance tensor components. This model would give a misfit of 5.52 for the sites located in profile 3, improving the misfit obtained by the previous 2D inversion. The fit of the model to the determinant responses, as well as to the antidiagonal components are shown for comparison in Fig. 8. To ensure the stability of the model
Fig. 6. 2D inversion of the determinant of the impedance tensor for profiles 1, 2 and 3.
C. Arango et al. / Journal of Applied Geophysics 68 (2009) 479–488
485
Fig. 7. Horizontal slices showing apparent resistivity distribution at different depths obtained from 3D model. Only the part of the model close to the sites is shown.
4. Discussion
(180 Ω m) can be associated with the upper aquifer, with an averaged thickness of 140 m, it is related to the Quaternary group and reef complex of the Upper Miocene. 2) The conductive structure (5–20 Ω m) with variable thickness may have a correlation with the confining unit or aquitard. The minimum thickness of this layer (around 200 m) is located below the thermal anomaly zone. 3) The other resistive media (N80 Ω m) are associated with the deeper aquifer unit, which supplies warm water. The main feature suggested by the 3D model is the presence of a zone where the conductor is thinner with a NE–SW preferential orientation, coincident with the thermal anomaly and with the general orientation of the regional structures. This thinner layer makes up the basement to the north and down towards the south, and may be related to a weakness or faulting zone. This model was tested with the results obtained from an exploratory well that was drilled at the end of 2002, based on preliminary AMT results and reaching 700 m depth. We compared the geological and log data information obtained from the well with the geoelectrical image determined by the 3D MT model (Fig. 9). Some correlations were derived from the log data recorded into the borehole:
The results from the 3D modeling allow us to quantify the suggested conceptual model: 1) The high, shallow resistive layer
1) The boundary between the shallow resistive zone (R1) and the less resistive domain (C) coincides with the limit between shelly
with regard to discretization of the mesh, we compared responses obtained from the model with those obtained for a more refined model (50 × 50 × 30) not observing significant differences. This model shows a resistivity distribution composed by three main layers: 1) A high resistive top layer, identified by R1, and evident in the 0 m to 100 m horizons. 2) A more conductive middle layer, between 200 m and 500 m, corresponding to domain C1. 3) A relatively resistive layer in the bottom (domain R2). The main characteristic of the geoelectrical model is the thickness variation of the conductive layer. This abrupt change is evident between 200 m and 300 m when this horizon disappears, showing the thinnest zone along NE–SW orientation, beneath the stations 6, 65 and 63. This behavior is evident until 600 m depth, where the structure becomes more resistive. In the 2D model, this abruptness is also detected in all profiles concurrent with the anomaly distribution. Nevertheless, the 3D model requires that the thickness of the conductive horizon be recovered towards to west to adjust the site 10.
486
C. Arango et al. / Journal of Applied Geophysics 68 (2008) 479–488
C. Arango et al. / Journal of Applied Geophysics 68 (2009) 479–488
487
Fig. 9. Comparison between lithology and log data observed into the exploration well drilled by IGME and 3D resistivity model.
limestone (reef complex) and the confining unit or aquitard composed meanly of marls. The low permeability of this type of rock is evident in the gamma ray log which reflects the shale content of formations, and enhances lithological interpretation. 2) The level R1 is not identified in normal resistivity logs, because its recording started at 140 m. However, its thickness coincides with information provided by lithology and the gamma ray log. 3) The resistivity value of the 3D model for domain C around 15 Ω m is similar to the value observed in the resistivity log measured into the borehole. 4) The resistivity logs show an increment of an order of magnitude at around 260 m. This variation could be related to the change from domain C to R2.
5) Below this conductive layer a resistant domain, R2, with a resistivity value of 60 Ω m can be observed. Inside this layer there are small variations that cannot be detected in the 3D model because the size of the model cells does not allow enough resolution. Data integration allows for detection of a 10 m thick layer around 245 m in depth. According to lithological analysis, it is a tectonic calcareous breccia which could be associated with the structure that allows warm water to rise. This variation cannot be detected by magnetotelluric sounding because of resolution, but by checking the 3D model results in the borehole zone, a sudden change in resistivity value is detected at this depth. This layer could be a plausible way for water circulation towards the upper aquifer, which is then exploited on surface.
Fig. 8. Curves of apparent resistivity and phase (points) and fit of 3D model (solid lines) for four representative sites. Left: data from the determinant responses; right: data from the antidiagonal components of the impedance tensor.
488
C. Arango et al. / Journal of Applied Geophysics 68 (2008) 479–488
5. Conclusions The 3D geoelectrical model confirms the conceptual one suggested for this zone, and allows us to identify the main features of the hydrological scheme. It is composed of a shallow unconfined aquifer, associated with reefal materials of Upper Miocene, as well as a deeper confined aquifer, related to Paleogene materials. They are separated by a level with poor permeability and low resistivity as well as a minimum thickness of around 200 m. The presence of hot water on the surface is related to this thinner zone of the aquitard. The integration of the MT model and well data suggest that this thinner part is associated with a fault or weakness zone located beneath the thermal manifestation, favoring vertical water movements. It has a preferential NE–SW alignment, coinciding with the regional structural direction. Acknowledgements This work was supported by the Ministerio de Ciencia y Tecnología, Spain, project REN 2002-4538-C0201 and by contract FBG-301994 between the Fundación Bosch Gimpera (U.B) and the Instituto Geológico y Minero de España. The stay of C. A. at the University of Barcelona was supported by the Consejo Nacional de Ciencia y Tecnología, (CONACYT), México. We would like to thank A. Martí, B. Gelabert and J. L. Plata, who provided helpful comments on this work, and R. M. Mateos, who provided with the temperature data and well information. We would also like to thank F. Santos and E. Almeida for equipment support. Lastly we thank our guest editors Gad El-Qady and Toivo Korja, as well as the two anonymous reviewers, for their comments for improvement of the manuscript. References Arango, C., 2005. Estudio Magnetotelúrico de la zona de Llucmajor (Mallorca): Avances en el proceso de datos y modelo 3D. Ph. D. Thesis. Departament de Geodinàmica i Geofísica. Universitat de Barcelona. Ayala, C., Pous, J., Sàbat, F., Casas, A., Rivero, L., Gelabert, B., 1994. Modelización gravimétrica de la isla de Mallorca. Rev. Soc. Geol. Esp. 7, 3–4. Benedicto, A., Ramos-Guerrero, E., Casas, A., Sàbat, F., Barón, A., 1993. Evolución tectonosedimentaria de la cubeta neógena de Inca (Mallorca). Rev. Soc. Geol. Esp. 6 (1–2), 167–176. Berdichevsky, M.N., Dmitriev, V.I., 1976. Distortion of magnetic and electric fields by near surface lateral inhomogeneities. Acta Geod. Geophys. Montan. Acad. Sci. Hung. 11, 447–483. Fernàndez, M., Cabal, J., 1992. Heat-flow data and shallow thermal regime on Mallorca and Menorca (western Mediterranean). Tectonophysics 203, 133–143. Friedrichs, B., 2001. MT Software-MAPROS. CooperPower Tools, Branschweiger, Germany.
Gelabert, B., 1998. La estructura geológica de la mitad occidental de la Isla de Mallorca. Instituto Tecnológico Geominero de España, Colección de memorias. Groom, R.W., Bailey, R.C., 1989. Decomposition of the magnetotelluric impedance tensor in the presence of local three-dimensional galvanic distortion. J. Geophys. Res. 94, 1913–1925. IGME, 1985. Investigación eléctrica en Lluchmayor para apoyo a estudios geotérmicos (Mallorca). Instituto Geológico y Minero de España, IGME. Document # 40270 (Dataset can be downloaded from www.igme.es/internet/sigeof/inicio_eng.html). Krivochieva, S., Chouteau, M., 2003. Integrating TDEM and MT methods for characterization and delineation of the Santa Catarina aquifer (Chalco Sub-Basin, Mexico). J. Appl. Geophys. 52, 23–43. Ledo, J., 2005. 2-D versus 3-D magnetotelluric data interpretation. Surv. Geophys. 26 (5), 511–543. Mackie, R.L., Madden, T.R., Wannamaker, P.E., 1993. Three-dimensional magnetotelluric modeling using difference equations — theory and comparisons to integral equation solutions. Geophysics 58 (2), 215–226. Martí, A., Queralt, P., Roca, E., 2004. Geoelectric dimensionality in complex geological areas: application to the Spanish Betic Chain. Geophys. J. Int. 157, 961–974. McNeice, G.W., Jones, A.G., 2001. Multisite, multifrequency tensor decomposition of magnetotelluric data. Geophysics 66 (1), 158–173. Meju, M.A., 2002. Geoelectromagnetic exploration for natural resources: models case studies and challenges. Surv. Geophys. 23, 133–205. Meju, M.A., Fontes, S.L., Oliveira, M.F.B., Lima, J.P.R., Ulugergerli, E.U., Carrasquilla, A.A., 1999. Regional aquifer mapping using combined VES-TEM-AMT/ EMAP methods in the semiarid eastern margin of Parnaiba Basin, Brazil. Geophysics 64, 337–356. Morlet, J., Arens, G., Fourgeau, E., Giard, D., 1982. Wave propagation and sampling theory — part 2: sampling theory and complex waves. Geophysics 47 (2), 222–236. Ogawa, Y., Matsushima, N., Oshima, H., Takakura, S., Utsugi, M., Hirano, K., Igarashi, M., Doi, T., 1998. A resistivity cross-section of Usu volcano, Hokkaido, Japan, by audiomagnetotelluric soundings. Earth Planets Space 50, 339–346. Pedersen, L.B., Engels, M., 2005. Routine 2D inversion of magnetotelluric data using the determinant of the impedance tensor. Geophysics 70 (2), G33–G41. Pomar, L., 1991. Reef geometries, erosion surfaces and high-frequency sea-level changes, upper Miocene reef complex, Mallorca, Spain. Sedimentology 38, 243–269. Pous, J., Ayala, C., Ledo, J., Marcuello, A., Sàbat, F., 1995. 3D modeling of magnetotelluric and gravity data of Mallorca island (Western Mediterranean). Geophys. Res. Lett. 22 (6), 735–738. Ritz, M., Descloitres, M., Robineau, B., Courteaud, M., 1997. Audiomagnetotelluric prospecting for groundwater in the Baril coastal area, Piton de la Fournaise Volcano, Reunion Island. Geophysics 62, 758–762. Siripurnvaraporn, W., Egbert, G., 2000. An efficient data subspace inversion method for 2-D magnetotelluric data. Geophysics 65, 791–803. Tezkan, B., 1999. A review of environmental applications of quasi-stationary electromagnetic techniques. Surv. Geophys. 20 (3–4), 279–308. Trad, D.O., Travassos, J.M., 2000. Wavelet filtering of magnetotelluric data. Geophysics 65, 482–491. Volpi, G., Manzella, A., Fiordelisi, A., 2003. Investigation of geothermal structures by magnetotellurics (MT): an example from the Mt. Amiata area, Italy. Geothermics 32, 131–145. Wannamaker, P.E., 1997. Tensor CSAMT survey over the Sulphur Springs thermal area, Valles Caldera, New Mexico, U.S.A.; part II, implications for CSAMT methodology. Geophysics 62, 466–476. Weaver, J.T., Agarwal, A.K., Lilley, F.E.M., 2000. Characterisation of the magnetotelluric tensor in terms of its invariants. Geophys. J. Int. 141, 321–336. Zhang, Y., Paulson, K.V., 1997. Enhancement of signal-to-noise ratio in natural-source transient magnetotelluric data with wavelet transform. Pure Appl. Geophys. 149, 405–419.
Contents lists available at ScienceDirect
Journal of Applied Geophysics j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / j a p p g e o
3D magnetotelluric characterization of the geothermal anomaly in the Llucmajor aquifer system (Majorca, Spain) C. Arango ⁎, A. Marcuello, J. Ledo, P. Queralt Departament de Geodinàmica i Geofísica, Facultat de Geología, Universitat de Barcelona, C/Martí i Franquès s/n, 08028 Barcelona, Spain
A R T I C L E
I N F O
Article history: Received 15 April 2008 Accepted 23 May 2008 Keywords: Llucmajor aquifer Magnetotellurics Geothermal exploration Geoelectrical structure 3D modeling
A B S T R A C T In the Llucmajor aquifer system (Majorca Island, Spain) some geothermal evidences have appeared. This phenomenon is not isolated to Majorca and it is present in other areas, where it can be associated with structural conditions, especially to the extensional event suffered by the island after the Alpine Orogeny. However, the origin of this anomaly in Llucmajor is not well known, and there is no surface geological evidence of these structural conditions. With the aim of delineating the geoelectrical structure of the zone and identifying the geological structure that allows the presence of this anomaly, an audiomagnetotelluric (AMT) survey was carried out. The AMT data was processed using a Wavelet Transform-based scheme. Dimensionality analysis indicates that the geoelectrical structure is mainly 3D. The 3D model was obtained by trial and error forward modeling, taking accounting of the responses from the determinant of the impedance tensor. The model shows a vertical resistivity distribution with three horizons associated with different units: on the top, a shallow high resistive media related to an unconfined shallow aquifer; in the middle, a conductive layer related to the aquitard, and below it, another resistive media related to the confined deeper aquifer. The intermediate horizon shows a sudden thinning beneath the thermal anomalous zone that can be identified as a weakness zone (fault or fracture) connecting both aquifers. An exploratory well was drilled after the AMT survey and reached almost 700 m in depth. This allowed correlating the resistivity distribution of the 3D model with data logging and lithology obtained from the well, showing a proper agreement between them. © 2008 Published by Elsevier B.V.
1. Introduction The Llucmajor aquifer system is located in the southern part of Majorca Island (Spain). In this area water is extracted from the unconfined aquifer of the system, and used for agricultural purposes. From the exploitation wells the inferred thickness of the aquifer is around 150 m. A geothermal anomaly has been detected in this area, manifested by anomalous high water temperatures in some wells aligned along a NS–EW direction, and reaching a maximum temperature of 50 °C (Fig. 1). Because of the lack of volcanic sources in this region, these high values must be associated to vertical movements of water. Taking the reported thermal gradient in the zone as 46 ± 8 mK/m (Fernàndez and Cabal, 1992), a minimum depth of 700 m is needed to obtain such a temperature. Understanding the origin of this hot water, as well as its possible economical exploitation was the aim of a research program conducted by the IGME (Spanish Geological Survey) and supported by the Balearic Government. To explain this behavior a simple conceptual model with two aquifer units was proposed: one unconfined aquifer
⁎ Corresponding author. Departamento de Geomagnetismo y Exploración, Instituto de Geofísica, UNAM, 04510 México, D.F., Mexico. E-mail address: claudiar@geofisica.unam.mx (C. Arango). 0926-9851/$ – see front matter © 2008 Published by Elsevier B.V. doi:10.1016/j.jappgeo.2008.05.006
on top, and another confined aquifer to supply the hot water. The aquitard between both aquifers must present a discontinuity (fault or fracture) to favor vertical fluid movements. This area is composed of Quaternary silts, and there is no surface geomorphologic evidence of structures that allows the rising of deep groundwater. In the studied zone, two geophysical surveys were reported. The first one employed vertical electrical soundings (IGME, 1985), but a high resistive shallower layer did not allow good penetration of the electrical current, and limited the range of the investigation depth to the upper aquifer. The second geophysical survey is related to the seismic reflection campaign carried out in Majorca between 1984 and 1990 (Benedicto et al., 1993). In this campaign, one of the profiles went through Llucmajor zone, but did not unfortunately have enough resolution to study the origin of the thermal anomaly. Taking this information into account, the magnetotelluric exploration was considered the proper technique for generating a model to explain the origin of the geothermal anomaly. Among the electromagnetic methods, the audiomagnetotellurics (AMT) method is appropriated for deep groundwater and geothermal studies (Tezkan, 1999; Meju, 2002) given the depth of its penetration and its tensorial nature. Over recent years, different studies of deep aquifers (Ritz et al., 1997; Wannamaker, 1997; Meju et al., 1999; Krivochieva and Chouteau, 2003) and geothermal problems (Ogawa et al., 1998; Volpi et al., 2003) have been reported.
480
C. Arango et al. / Journal of Applied Geophysics 68 (2008) 479–488
Fig. 1. Temperature distribution in producing water wells of the Llucmajor zone.
New methodological aspects have recently been developed in magnetotellurics, in particular: the use of Wavelet Transform for time series processing (Zhang and Paulson, 1997; Trad and Travassos, 2000; Arango, 2005); and the use of the invariants to analyze the geoelectrical dimensionality (Weaver et al., 2000; Martí et al., 2004). This work incorporates them into the study of this aquifer. The main objectives are: 1) to determine the 3D geometry of the aquifer system and 2) to identify the zone where the aquitard could be fractured. We first describe the geology of the area. We then describe the main aspects of the magnetotelluric method, emphasizing the methodological elements such as geoelectrical dimensionality and the use of the determinant of the impedance tensor in the interpretation. Later, the model is correlated with data from the exploratory well, drilled following the MT fieldwork. 2. Geological setting Majorca is the biggest island of the Balearic archipelago and is considered the NE continuation of the Betic chain, located in southern Spain. The main structures of the island are identified by three mountain ranges with a NE–SW direction (Tramuntana, Randa and Levante systems) and separated by two main depressions (Gelabert, 1998). Its current configuration is conditioned by Neogene distension, which formed small basins. Three of them (Palma, Inca and Sa Pobla) are aligned in the western area, limited by NE–SW faults and separated from each other by NW–SE structures. The fourth (Llucmajor–Campos), in the eastern area, shows the same fault pattern (Fig. 2). The Llucmajor platform is part of one of these depressions. It is mainly composed of Quaternary materials (red silts, conglomerates and aeolianites), and Upper Miocene materials (limestones). Reefal calcarenites and discontinuous Quaternary silts crop out in the northern zone (Pomar, 1991). Several regional geophysical surveys have been carried out in Majorca using gravity, magnetic, seismic and MT methods (Benedicto et al., 1993; Ayala et al., 1994; Pous et al., 1995). They confirm that the main regional structures observed on surface continue beneath the island.
3. MT data survey 3.1. The MT method The magnetotelluric method is a geophysical technique that determines ground electrical resistivity distribution from the simultaneous measurement of the fluctuations of the natural electromagnetic field. Under the plane wave assumption, the linear relationship between the horizontal components of electric and magnetic fields at a given frequency, f, are expressed as follow: Eð f Þ = Z ð f ÞHð f ÞZ
Ex Ex
=
Zxx Zyx
Zxy Zyy
Hx Hx
ð1Þ
where Z is the impedance tensor. The impedance is a complex magnitude, from which it is customary to define the apparent resistivity and the phase for each component of the tensor as: ρaij ð f Þ =
1 2 jZ ð f Þj 2πf μ ij
! − 1 Im Zij ð f Þ uij ð f Þ = arg Zij ð f Þ = tan Re Zij ð f Þ
ð2aÞ
ð2bÞ
where µ is the magnetic permeability, and i, j denote any horizontal component. In this work we have also used the apparent resistivity and phase calculated from the determinant of the impedance tensor, introduced by Berdichevsky and Dmitriev (1976) as effective impedance: ρa ð f Þ =
uð f Þ =
1 jdet Z ð f Þj 2πf μ 1 argðdet Z ð f ÞÞ 2
ð3aÞ
ð3bÞ
These responses have interesting properties for modeling: they are rotationally invariant, the phase is not affected by galvanic distortions,
C. Arango et al. / Journal of Applied Geophysics 68 (2009) 479–488
Fig. 2. Geological setting of the Majorca Island (Spain) modified from Gelabert (1998).
Fig. 3. Location of AMT and MT sites and interpreted profiles.
481
482
C. Arango et al. / Journal of Applied Geophysics 68 (2008) 479–488
and both the apparent resistivity and the phase contain information from the four complex elements of the impedance tensor. Its use simplifies the data employed in the 3D modeling and integrates information of all components of the impedance tensor. 3.2. Magnetotelluric data The fieldwork was designed following two main perpendicular profiles, the inner section being filled out using a grid to satisfactorily cover the area of interest (Fig. 3). The average distance between sites was 250 m. Data acquisition had been carried out in summer 2002 at 42 different sites with two different Metronix systems. 36 measurements were taken with the GSM06 system, whose frequency ranges from 8 · 103 to 2.5 · 10− 3 Hz.10 measurements were taken with the MMS03E system, whose frequency ranges from 2.5 · 102 to 2.5 · 10− 3 Hz. Four sites were recorded with both instruments to test the agreement between results. The magnetotelluric horizontal components had been measured along the N–S and E–W directions. The recording time of most of the AMT stations were only a few hours, so only responses for higher frequencies (N0.1–1 Hz) were computed, and the MT stations covered lower frequencies. Although simultaneous recordings between both systems were performed, no remote reference was possible, since synchronization was unfeasible.
The vertical magnetic component presented electronic problems, producing very scattered tipper responses, and was not considered. For an estimation of impedance from the time series we employed the Wavelet Transform algorithm as it is described in Zhang and Paulson (1997), and implemented and modified by Arango (2005). This technique allows for a detailed analysis in the time-frequency domain, as well as a description of frequency content along time, which is not possible using the classic method based on the Fourier analysis. The mother wavelet used in this work to compute the Wavelet Transform of the electromagnetic fields was the function proposed by Morlet et al. (1982). It was chosen for its capacity to concentrate energy in a narrow frequency band, as well as its simple lineal relationship between frequency and scale. The most important feature of this processing method is the selection of intervals with high energy content for each scale (frequency) as suggested in Zhang and Paulson (1997). These are used to calculate the impedance tensor elements to ensure a better signal in the dead band. After several attempts and for this work, the threshold level between high and low energy content was taken as the mean energy. This offers a proper tradeoff for our data quality, as shown in Fig. 4, where a comparison between the Wavelet Transform algorithm and a standard process is presented. Both processes present similar results, although the Wavelet Transform-based appears as less scattered.
Fig. 4. Apparent resistivity and phase curves computed using standard scheme(Mapros-software, Friedrichs, 2001) and the Wavelet Transform-based scheme.
C. Arango et al. / Journal of Applied Geophysics 68 (2009) 479–488
483
Fig. 5. (a) Directionality analysis for different frequency intervals using the algorithm by Martí et al. (2004). (b) Strike direction computed between 1 and 100 Hz.
The data is restricted to frequencies higher than 0.1 Hz, focusing on the shallow part of the geoelectrical structure. The sea is located 10– 15 km away. This is one to two times the skin depth, if both the frequency range, and a non-resistive basement in this region (Pous et al., 1995), are considered. Therefore there is no coast effect in the studied zone for the chosen range of frequencies.
An important step prior to interpretation is dimensionality analysis, which allows identification of the main directions of geoelectrical structures. The analysis of dimensionality used in this work is based on the invariants of the magnetotelluric impedance tensor, as described in Weaver et al. (2000), and implemented following the algorithm of Martí et al. (2004). This analysis showed
484
C. Arango et al. / Journal of Applied Geophysics 68 (2008) 479–488
that the behavior of the geoelectrical structure between 0.1 and 1000 Hz is three-dimensional, and supported this description with a 3D model (Fig. 5a). It can be also seen in this figure that several sites exhibit 2D behavior between 1 and 100 Hz, and this is employed to obtain a 2D model, used as preparatory model for the further 3D interpretation. According to Ledo (2005), in some cases it is possible to extract useful information from 3D data using a 2D algorithm, but this approach must be taken carefully however, since it depends on both geoelectrical configuration and the particular features of the media. Furthermore, a determination of the geoelectrical strike was also done to all sites using the Strike program (McNeice and Jones, 2001) based on the Groom–Bailey decomposition (Groom and Bailey, 1989). In this frequency interval, this analysis indicated a strike of N50°E (Fig. 5b). Static shift was corrected, assuming that there were no important variations on the very shallow structures, as suggested from the geological information. It was thus expected that the near surface resistivity would be essentially the same at all the sites. This assumption is also supported by the results of the vertical electrical soundings in this area (IGME, 1985). Most of the apparent resistivities coincided with the same value at high frequencies, taken as the average, which is taken as the near surface resistivity. Three perpendicular profiles along the NW–SE direction and identified as profiles 1, 2 and 3 (Fig. 3) were considered. The 2D resistivity models were computed using the REBOCC inversion algorithm (Siripurnvaraporn and Egbert, 2000), modified by Pedersen and Engels (2005) to use the determinant of the impedance tensor. This type of data ensures a better data fitting when a strike is not fully determined, as is the situation with 2D interpretation of 3D structures. The final 2D models are shown in Fig. 6 and the integrated rms value for each profile was 2.97, 3.14 and 8.61 respectively, with a 5% error floor in the impedance components. The higher misfit in profile 3 is probably caused by bad data quality, and one of the sites presents a strong 3D behavior (Fig. 5a).
In the 2D models we can identify three layers with different electrical properties along all profiles. The intermediate horizon has a low resistivity value, unlike the other two layers, and its thickness varies suddenly beneath the thermal zone. These results were used for creating a 3D model grid, as well as for building the starting 3D model. We assigned resistivity values obtained from 2D models to create a 3D model grid with information based on previous results. 3.3. 3D model The acquisition pattern designed allowed the construction of a 3D model to characterize electrical resistivity distribution. Such a model consisted of 35 × 35 cells in NS–EW between 115 and 660 m and 20 layers of variable thickness between 20 and 250 m. This configuration ensured that grid is refined over zones containing more information. To account for the boundary conditions, the grid has been laterally extended twice the size of the study area. The procedure for obtaining the model was established through a trial and error strategy with a varying resistivity value of the cells using the determinant of the impedance tensor as a model response for data comparison. Once such model is found, it is tested for adjusting the responses of the antidiagonal elements of the impedance tensor. Forward modeling was computed using the algorithm developed by Mackie et al. (1993), and modified by Booker and Mackie (personal communication, 1999). The data frequencies ranged from 1000 Hz to 0.1 Hz. The 3D geoelectrical model that we obtained shows the main subsurface structures up to 600 m in depth (Fig. 7). The fit between data and model responses was satisfactory, and the rms was 6.07 with an error floor of 5% in the impedance tensor components. This model would give a misfit of 5.52 for the sites located in profile 3, improving the misfit obtained by the previous 2D inversion. The fit of the model to the determinant responses, as well as to the antidiagonal components are shown for comparison in Fig. 8. To ensure the stability of the model
Fig. 6. 2D inversion of the determinant of the impedance tensor for profiles 1, 2 and 3.
C. Arango et al. / Journal of Applied Geophysics 68 (2009) 479–488
485
Fig. 7. Horizontal slices showing apparent resistivity distribution at different depths obtained from 3D model. Only the part of the model close to the sites is shown.
4. Discussion
(180 Ω m) can be associated with the upper aquifer, with an averaged thickness of 140 m, it is related to the Quaternary group and reef complex of the Upper Miocene. 2) The conductive structure (5–20 Ω m) with variable thickness may have a correlation with the confining unit or aquitard. The minimum thickness of this layer (around 200 m) is located below the thermal anomaly zone. 3) The other resistive media (N80 Ω m) are associated with the deeper aquifer unit, which supplies warm water. The main feature suggested by the 3D model is the presence of a zone where the conductor is thinner with a NE–SW preferential orientation, coincident with the thermal anomaly and with the general orientation of the regional structures. This thinner layer makes up the basement to the north and down towards the south, and may be related to a weakness or faulting zone. This model was tested with the results obtained from an exploratory well that was drilled at the end of 2002, based on preliminary AMT results and reaching 700 m depth. We compared the geological and log data information obtained from the well with the geoelectrical image determined by the 3D MT model (Fig. 9). Some correlations were derived from the log data recorded into the borehole:
The results from the 3D modeling allow us to quantify the suggested conceptual model: 1) The high, shallow resistive layer
1) The boundary between the shallow resistive zone (R1) and the less resistive domain (C) coincides with the limit between shelly
with regard to discretization of the mesh, we compared responses obtained from the model with those obtained for a more refined model (50 × 50 × 30) not observing significant differences. This model shows a resistivity distribution composed by three main layers: 1) A high resistive top layer, identified by R1, and evident in the 0 m to 100 m horizons. 2) A more conductive middle layer, between 200 m and 500 m, corresponding to domain C1. 3) A relatively resistive layer in the bottom (domain R2). The main characteristic of the geoelectrical model is the thickness variation of the conductive layer. This abrupt change is evident between 200 m and 300 m when this horizon disappears, showing the thinnest zone along NE–SW orientation, beneath the stations 6, 65 and 63. This behavior is evident until 600 m depth, where the structure becomes more resistive. In the 2D model, this abruptness is also detected in all profiles concurrent with the anomaly distribution. Nevertheless, the 3D model requires that the thickness of the conductive horizon be recovered towards to west to adjust the site 10.
486
C. Arango et al. / Journal of Applied Geophysics 68 (2008) 479–488
C. Arango et al. / Journal of Applied Geophysics 68 (2009) 479–488
487
Fig. 9. Comparison between lithology and log data observed into the exploration well drilled by IGME and 3D resistivity model.
limestone (reef complex) and the confining unit or aquitard composed meanly of marls. The low permeability of this type of rock is evident in the gamma ray log which reflects the shale content of formations, and enhances lithological interpretation. 2) The level R1 is not identified in normal resistivity logs, because its recording started at 140 m. However, its thickness coincides with information provided by lithology and the gamma ray log. 3) The resistivity value of the 3D model for domain C around 15 Ω m is similar to the value observed in the resistivity log measured into the borehole. 4) The resistivity logs show an increment of an order of magnitude at around 260 m. This variation could be related to the change from domain C to R2.
5) Below this conductive layer a resistant domain, R2, with a resistivity value of 60 Ω m can be observed. Inside this layer there are small variations that cannot be detected in the 3D model because the size of the model cells does not allow enough resolution. Data integration allows for detection of a 10 m thick layer around 245 m in depth. According to lithological analysis, it is a tectonic calcareous breccia which could be associated with the structure that allows warm water to rise. This variation cannot be detected by magnetotelluric sounding because of resolution, but by checking the 3D model results in the borehole zone, a sudden change in resistivity value is detected at this depth. This layer could be a plausible way for water circulation towards the upper aquifer, which is then exploited on surface.
Fig. 8. Curves of apparent resistivity and phase (points) and fit of 3D model (solid lines) for four representative sites. Left: data from the determinant responses; right: data from the antidiagonal components of the impedance tensor.
488
C. Arango et al. / Journal of Applied Geophysics 68 (2008) 479–488
5. Conclusions The 3D geoelectrical model confirms the conceptual one suggested for this zone, and allows us to identify the main features of the hydrological scheme. It is composed of a shallow unconfined aquifer, associated with reefal materials of Upper Miocene, as well as a deeper confined aquifer, related to Paleogene materials. They are separated by a level with poor permeability and low resistivity as well as a minimum thickness of around 200 m. The presence of hot water on the surface is related to this thinner zone of the aquitard. The integration of the MT model and well data suggest that this thinner part is associated with a fault or weakness zone located beneath the thermal manifestation, favoring vertical water movements. It has a preferential NE–SW alignment, coinciding with the regional structural direction. Acknowledgements This work was supported by the Ministerio de Ciencia y Tecnología, Spain, project REN 2002-4538-C0201 and by contract FBG-301994 between the Fundación Bosch Gimpera (U.B) and the Instituto Geológico y Minero de España. The stay of C. A. at the University of Barcelona was supported by the Consejo Nacional de Ciencia y Tecnología, (CONACYT), México. We would like to thank A. Martí, B. Gelabert and J. L. Plata, who provided helpful comments on this work, and R. M. Mateos, who provided with the temperature data and well information. We would also like to thank F. Santos and E. Almeida for equipment support. Lastly we thank our guest editors Gad El-Qady and Toivo Korja, as well as the two anonymous reviewers, for their comments for improvement of the manuscript. References Arango, C., 2005. Estudio Magnetotelúrico de la zona de Llucmajor (Mallorca): Avances en el proceso de datos y modelo 3D. Ph. D. Thesis. Departament de Geodinàmica i Geofísica. Universitat de Barcelona. Ayala, C., Pous, J., Sàbat, F., Casas, A., Rivero, L., Gelabert, B., 1994. Modelización gravimétrica de la isla de Mallorca. Rev. Soc. Geol. Esp. 7, 3–4. Benedicto, A., Ramos-Guerrero, E., Casas, A., Sàbat, F., Barón, A., 1993. Evolución tectonosedimentaria de la cubeta neógena de Inca (Mallorca). Rev. Soc. Geol. Esp. 6 (1–2), 167–176. Berdichevsky, M.N., Dmitriev, V.I., 1976. Distortion of magnetic and electric fields by near surface lateral inhomogeneities. Acta Geod. Geophys. Montan. Acad. Sci. Hung. 11, 447–483. Fernàndez, M., Cabal, J., 1992. Heat-flow data and shallow thermal regime on Mallorca and Menorca (western Mediterranean). Tectonophysics 203, 133–143. Friedrichs, B., 2001. MT Software-MAPROS. CooperPower Tools, Branschweiger, Germany.
Gelabert, B., 1998. La estructura geológica de la mitad occidental de la Isla de Mallorca. Instituto Tecnológico Geominero de España, Colección de memorias. Groom, R.W., Bailey, R.C., 1989. Decomposition of the magnetotelluric impedance tensor in the presence of local three-dimensional galvanic distortion. J. Geophys. Res. 94, 1913–1925. IGME, 1985. Investigación eléctrica en Lluchmayor para apoyo a estudios geotérmicos (Mallorca). Instituto Geológico y Minero de España, IGME. Document # 40270 (Dataset can be downloaded from www.igme.es/internet/sigeof/inicio_eng.html). Krivochieva, S., Chouteau, M., 2003. Integrating TDEM and MT methods for characterization and delineation of the Santa Catarina aquifer (Chalco Sub-Basin, Mexico). J. Appl. Geophys. 52, 23–43. Ledo, J., 2005. 2-D versus 3-D magnetotelluric data interpretation. Surv. Geophys. 26 (5), 511–543. Mackie, R.L., Madden, T.R., Wannamaker, P.E., 1993. Three-dimensional magnetotelluric modeling using difference equations — theory and comparisons to integral equation solutions. Geophysics 58 (2), 215–226. Martí, A., Queralt, P., Roca, E., 2004. Geoelectric dimensionality in complex geological areas: application to the Spanish Betic Chain. Geophys. J. Int. 157, 961–974. McNeice, G.W., Jones, A.G., 2001. Multisite, multifrequency tensor decomposition of magnetotelluric data. Geophysics 66 (1), 158–173. Meju, M.A., 2002. Geoelectromagnetic exploration for natural resources: models case studies and challenges. Surv. Geophys. 23, 133–205. Meju, M.A., Fontes, S.L., Oliveira, M.F.B., Lima, J.P.R., Ulugergerli, E.U., Carrasquilla, A.A., 1999. Regional aquifer mapping using combined VES-TEM-AMT/ EMAP methods in the semiarid eastern margin of Parnaiba Basin, Brazil. Geophysics 64, 337–356. Morlet, J., Arens, G., Fourgeau, E., Giard, D., 1982. Wave propagation and sampling theory — part 2: sampling theory and complex waves. Geophysics 47 (2), 222–236. Ogawa, Y., Matsushima, N., Oshima, H., Takakura, S., Utsugi, M., Hirano, K., Igarashi, M., Doi, T., 1998. A resistivity cross-section of Usu volcano, Hokkaido, Japan, by audiomagnetotelluric soundings. Earth Planets Space 50, 339–346. Pedersen, L.B., Engels, M., 2005. Routine 2D inversion of magnetotelluric data using the determinant of the impedance tensor. Geophysics 70 (2), G33–G41. Pomar, L., 1991. Reef geometries, erosion surfaces and high-frequency sea-level changes, upper Miocene reef complex, Mallorca, Spain. Sedimentology 38, 243–269. Pous, J., Ayala, C., Ledo, J., Marcuello, A., Sàbat, F., 1995. 3D modeling of magnetotelluric and gravity data of Mallorca island (Western Mediterranean). Geophys. Res. Lett. 22 (6), 735–738. Ritz, M., Descloitres, M., Robineau, B., Courteaud, M., 1997. Audiomagnetotelluric prospecting for groundwater in the Baril coastal area, Piton de la Fournaise Volcano, Reunion Island. Geophysics 62, 758–762. Siripurnvaraporn, W., Egbert, G., 2000. An efficient data subspace inversion method for 2-D magnetotelluric data. Geophysics 65, 791–803. Tezkan, B., 1999. A review of environmental applications of quasi-stationary electromagnetic techniques. Surv. Geophys. 20 (3–4), 279–308. Trad, D.O., Travassos, J.M., 2000. Wavelet filtering of magnetotelluric data. Geophysics 65, 482–491. Volpi, G., Manzella, A., Fiordelisi, A., 2003. Investigation of geothermal structures by magnetotellurics (MT): an example from the Mt. Amiata area, Italy. Geothermics 32, 131–145. Wannamaker, P.E., 1997. Tensor CSAMT survey over the Sulphur Springs thermal area, Valles Caldera, New Mexico, U.S.A.; part II, implications for CSAMT methodology. Geophysics 62, 466–476. Weaver, J.T., Agarwal, A.K., Lilley, F.E.M., 2000. Characterisation of the magnetotelluric tensor in terms of its invariants. Geophys. J. Int. 141, 321–336. Zhang, Y., Paulson, K.V., 1997. Enhancement of signal-to-noise ratio in natural-source transient magnetotelluric data with wavelet transform. Pure Appl. Geophys. 149, 405–419.