Imaging Thermal Anomalies in Hot Dry Rock Geothermal Systems from Near-Surface Geophysical Modelling

Abstract

Convective hydrothermal systems have been extensively studied using electrical and electromagnetic methods given the strong correlation between low conductivity anomalies associated with hydrothermal brines and high temperature areas. However, studies addressing the application of similar geophysical methods to hot dry rock geothermal systems are very limited in the literature. The Timanfaya volcanic area, located on Lanzarote Island (Canary Islands), comprises one of these hot dry rock systems, where ground temperatures ranging from 250 to 605 °C have been recorded in pyroclastic deposits at shallow (<70 m) depths. With the aim of characterizing the geophysical signature of the high ground temperature areas, three different geophysical techniques (ground penetrating radar, electromagnetic induction and magnetic prospecting) were applied in a well-known geothermal area located inside Timanfaya National Park. The area with the highest ground temperatures was correlated with the location that exhibited strong ground penetrating radar reflections, high resistivity values and low magnetic anomalies. Moreover, the high ground temperature imaging results depicted a shallow, bowl-shaped body that narrowed and deepened vertically to a depth greater than 45 m. The ground penetrating radar survey was repeated three years later and exhibited subtle variations of the signal reflection patterns, or signatures, suggesting a certain temporal variation of the ground temperature. By identifying similar areas with the same geophysical signature, up to four additional geothermal areas were revealed. We conclude that the combined use of ground penetrating radar, electromagnetic induction and magnetic methods constitutes a valuable tool to locate and study both the geometry at depth and seasonal variability of geothermal areas associated with hot dry rock systems.

1. Introduction

One of the most remarkable features of active volcanic areas is the presence of high-enthalpy geothermal systems that exhibit high temperature zones (>150–200 °C) at ground level (e.g., [1]). The most common are convective hydrothermal systems that consist of a magmatic body acting as a heat source, a groundwater system to transport the heat towards the surface and a confining, shallow, impermeable structure (e.g., [2]). Less common are hot dry rock (HDR) geothermal systems that consist of subsurface zones with very low fluid content (e.g., [3]). Most of the literature on the exploration of geothermal resources is related to hydrothermal systems and is primarily based on electrical and electromagnetic methods due to the strong correlation between the low conductivity anomalies associated with hydrothermal brines and the high temperature areas. Among the geophysical techniques more commonly used are magnetotellurics, transient electromagnetics, electrical resistivity, magnetics, self-potential and seismic, with investigation depths ranging from few to hundreds of meters (e.g., [1] and references herein). However, the use of such techniques for studying HDR systems is very limited because HDR systems are not as common as convective hydrothermal ones, and because the presence of vapour as dominant phase instead of hot water makes these kinds of systems more difficult to interpret using electromagnetic methods.

The Timanfaya volcanic area (Figure 1) is located on Lanzarote Island (Canary Archipelago, Spain). The Canary Islands constitute an intraplate oceanic volcanic archipelago consisting of seven major islands and four islets located in the eastern Central Atlantic, approximately 125 km from northwest Africa. The islands of Lanzarote and Fuerteventura developed in the earliest stage of the Canary magmatic activity and constitute the emergent part of the Eastern Canary Ridge, which is located northeast of the archipelago [4]. In Lanzarote, thermal anomalies associated with a HDR geothermal system are still present in different areas as a consequence of the historical eruptive activity that occurred between 1730 and 1736 (the Timanfaya eruption).

Over 30 volcanic vents formed during that eruption, mostly characterized by fissural eruptions and low-explosive basaltic magmas. The last eruption on Lanzarote Island occurred in 1824, with the formation of three new volcanoes. Typically, all of those vents are aligned along a system of fractures oriented to N70°E, of approximately 14 km in length, following the path of a central structural rift-type zone [5,6]. As a result of the 1730–1736 eruption, hot magma exists at approximately 4 km depth [7]. Although the magma body undergoes a cooling process, it still has a very high temperature, which explains the presence of superficial thermal anomalies in the area. The anomalies are isolated and limited to some fracture-related alignments and along craters’ rims [7,8,9].

The studied geothermal field is located on Islote de Hilario, within Timanfaya National Park, northwest of the Timanfaya volcanic system (see Figure 1 and Figure 2). The shallow geological structure of Islote de Hilario, determined from the previously available information [10] (surface geological units and several monitoring wells are used for temperature and gas composition determination), is very simple, consisting of a thick (>70 m) unit of pyroclastic materials that originated during the historical eruption of Timanfaya (1730–1736) overlying older basaltic lava flows. The measured temperature in this area is approximately 605 °C inside a 13-m deep well [10]. Other temperature measurements made in shallower wells range from 300 to 380 °C, 450 °C in a natural furnace and approximately 250 °C on the surface [11]. This geothermal field has been defined as a HDR system [8,11,12] that is connected to the base of the Timanfaya volcano through very shallow inclined fractures as revealed by previous surveys [13]. This work presents a study on the geothermal anomalies located at the Timanfaya volcanic area via an analysis and joint interpretation of ground penetrating radar (GPR), magnetic prospecting and electromagnetic induction (EMI) data obtained for areas not previously surveyed. First, we combined three shallow geophysical methods applied on a well-known geothermal anomaly (Figure 2) to fully characterize the geophysical signature and the relationship between the resistivity, magnetic and electromagnetic anomalies with temperature. Then, we have looked for the same geophysical signature allowing us to locate similar unknown geothermal anomalies in the surrounding area. The results help assess the reliability of the combination of the proposed methods to locate the geothermal zones at the surface and to obtain information regarding their distribution and geometry at depth.

2. Geophysical Techniques

2.1. Ground Penetrating Radar (GPR)

GPR is a geophysical method based on the propagation of electromagnetic pulses (1–20 ns) in the frequency band of 10 MHz–2.5 GHz. A transmitting antenna emits an electromagnetic signal into the ground, which is partly reflected when it encounters media with different dielectric properties and partly transmitted into deeper layers. Then, the reflections produced are recorded from the receiving antenna, which is either in a separate antenna box or in the same antenna box as the transmitter. The strength (amplitude) of the reflected fields is proportional to the magnitude of the dielectric constant change, which is characterized by the reflection coefficient (RC). The RC increases as the dielectric contrast increases. As the antenna is moved along the ground surface, a two-dimensional image (radargram) is obtained, which is an XZ graphic representation of the detected reflections. The x-axis represents the antenna displacement along the survey line, and the z-axis represents the two-way travel time of the pulse emitted (in terms of nanoseconds). If the time required to propagate to a reflector and back is measured, and the velocity of this pulse in medium is known, then the depth of the reflector can be determined [14,15,16].

GPR Data Acquisition and Processing

The same GPR transect was recorded in a time span of three years during two different field campaigns, May 2012 and April 2015 (Figure 2b). This was aimed to investigate possible variations on the subsurface temperature along a well-known geothermal anomaly (TA1) developed over pyroclastic deposits.

The GPR surveys were conducted using a RAMAC/GPR system from MALÅ Geosciences (Mala, Sweden). First, a common mid-point mode (CMP) test was performed to determine the velocity of propagation of the GPR waves in the subsurface medium, which is composed of lapilli. In the CMP configuration, the receiver (Rx) and transmitter (Tx) are moved one from the other at a constant step along the survey line, maintaining a common fixed centre point. The antenna spacing is varied at a fixed location and the change in the two-way travel time of the electromagnetic wave from the reflectors is measured [14]. The CMP data acquisition was carried out with a 200-MHz unshielded antenna. The geometry consisted of an 11-m antenna separation with a step size of 0.1 m (0.2 m antenna separation per trace), with reflections collected using a time window of 234 ns, defined by 512 samples per trace. A root-mean-squared velocity of 0.08 m/ns was obtained for the pyroclastic deposit (Figure 3).

Next, both survey seasons were conducted using a common-offset (CO) mode with a 100-MHz unshielded antenna. In the CO configuration, the receiver and transmitter are moved together along the radar line with a constant distance between them. In this work, point-to-point data were acquired with an offset between transmitter and receiver of 1 m, with reflections collected within a time window of 467 ns defined by 512 samples per trace and trace-distance intervals of 50 and 20 cm (in 2012 and 2015 campaigns, respectively). To measure the profile lengths, and to control the distance between traces, a tape measure was used as a guide (Figure 2b). The 100-MHz unshielded antenna was selected since this frequency provides sufficient penetration (up to 18 m depth) [14]. Moreover, unshielded antennas produce more powerful electromagnetic energy into the ground reaching higher depth penetration. However, they are very sensitive to background noise (e.g., electrical wires and metal) and should be used in cleaned spaces. For our time window, and considering a velocity of propagation of 0.08 m/ns, the maximum depth of penetration is about 18 m.

The reflection profiles recorded during the two different field campaigns were identically filtered using the processing sequence described in

Table 1. Data processing applied to the 100-MHz GPR data.

Step	Filter	Parameters
1	Time-zero correction	---
2	Subtract-mean (dewow)	Time window: 10 ns
3	Energy decay	Scaling value: 0.5
4	Subtracting average	Average traces: 100
5	Band-pass (butterworth)	Low pass: 40 MHz; high pass: 150 MHz
6	Topographic correction	Velocity: 0.08 m/ns

with the ReflexW software (version 8.5.7, Sandmeier geophysical research, Karlsruhe, Germany) [17]. The velocity of propagation obtained from the CMP (using the algorithm semblance analysis [18]) was used for topographic correction as well as to convert the time window of the radargrams into depth values (Figure 4).

2.2. Electromagnetic Survey

Electromagnetic induction (EMI) is a shallow geophysical method that obtains information about the apparent conductivity (or, its inverse, apparent resistivity) distribution at depth without requiring the direct contact of electrodes with the ground, as in the electrical resistivity method. EMI has been widely used to detect arsenic concentrations in groundwater [19], map saline plumes [20], locate faults [21], depict impact structures [22] and identify sinkholes [23], among other uses. Because temperature modifies the electromagnetic properties of the materials, the anomalies mapped by the EMI method might be related to variations in temperature, represented as differences in resistivity contrast.

Typical EMI equipment consists of two coils (a transmitter and receiver) separated by a fixed distance. A time-varying electric current at a fixed frequency generates a primary magnetic field at the transmitting coil. The primary field creates eddy currents in the ground that induce a secondary magnetic field. Both fields are recorded by the receiving coil [24]. From these data, the apparent resistivity of the subsoil can be obtained as it is linearly related to the quadrature component (i.e., the ratio of the secondary to the primary magnetic field) [24].

The depth of investigation depends on both the orientation of the coils (vertical or horizontal coplanar dipoles) and the intercoil spacing. We used a Geonics EM34-3 conductivity meter with three different intercoil spacings of 10 m, 20 m and 40 m, corresponding to transmitter frequencies of 6.4 kHz, 1.6 kHz and 0.4 kHz, respectively [24]. For each intercoil distance, we used two different coil configurations, vertical and horizontal coplanar dipoles (Figure 2c), selecting a station spacing of 1 m to attain a suitable horizontal resolution of the profiles. When the length of the profile is ≥40 m, a total of six measurements at each station can be obtained. The intercoil distances of 10 m, 20 m and 40 m correspond to investigation depths of 15 m, 30 m and 60 m, respectively, in the vertical dipole configuration and 7.5 m, 15 m and 30 m in the horizontal dipole configuration. Considering the availability of space within the study area, three profiles of different lengths and orientation were investigated (Figure 2a and Figure 5) in May 2013. The longest profile (EMI1), with a length of 49 m, runs NW-SE, laterally crossing a known geothermal anomaly area at the SE end (TA1) where the National Park’s visitor exhibitions are located (Figure 2c). The readings totalled 160 (80 for each dipole configuration). A second profile (EMI2) trending NE-SW, transverse to EMI1, also crosses the geothermal anomaly area (TA1) at its central part but is only 21 m long because it is limited by a parking area and other man-made walls. Because of space restrictions, an intercoil spacing of only 10 m was used with 26 readings. A third profile (EMI3) running N-S was obtained to image the resistivity values outside but close to the geothermal anomaly area (TA1). The length of this profile is 30 m, allowing the use of both 10 m and 20 m intercoil spacings, with readings totalling 64.

The field data were inverted using EM1DFM [25] software to produce a true resistivity depth image. EM1DFM is a 2D inversion programme that fits a modelled response to the data by minimizing an objective function [26]. The data required are the quadrature response of the secondary magnetic field, the orientation and intercoil separations, transmitter frequencies, height of the coils, and data error. The quadrature response data were inverted to a model consisting of 25 layers with thickness ranging from 0.5 to 2 m, increasing at depth. An L-curve inversion type, with a maximum number of iterations of 15 and a tolerance value of 0.4, was determined to be optimal after several trial inversions. The resulting inverted model presents the resistivity variations in both lateral and vertical directions along the profile. Because the data measured correspond to the midpoint of the different intercoil spacings, the lengths of the inverted models are less than the lengths of the corresponding profiles. Thus, the EMI1 model is 40 m long, the EMI2 model is 13 m long and the EMI3 model is 21 m long.

2.3. Magnetic Survey

Magnetic anomalies represent the contribution to the Earth’s magnetic field originating from geological structures due to the presence of small amounts of magnetic minerals in crustal rocks. These anomalies can be studied at very different scales, from a lithospheric perspective when they are measured from satellites to high-resolution local studies based on land surveys aimed at imaging the shallow subsurface structure. In particular, volcanic areas typically display intense magnetic anomalies due to the high concentration of magnetic minerals of volcanic rocks. Most surveys are conducted using scalar magnetometers that measure the magnitude of the Earth’s magnetic field vector. Under these conditions, the magnetic anomalies created by the subsurface structures are dipolar in shape, each of them consisting of a high and low (unless the anomalies are measured at the magnetic poles, where the Earth’s magnetic field is vertical).

Magnetization (J), which is the physical property responsible for magnetic anomalies, is a vector quantity resulting from the sum of a remanent component acquired in the past (J_rem) and an induced component (J_ind), which is proportional to the magnetic susceptibility (k) and is parallel to the Earth’s present magnetic field (H). Thus, the following relationship may be written: J = J_rem + J_ind = J_rem + k H. In volcanic rocks, the remanent magnetization is primarily acquired when the rocks are cooled in the presence of the Earth’s magnetic field (thermoremanence) and usually exceeds the induced magnetization. Therefore, in volcanic environments, the ratio J_rem/J_ind, known as the Köenigsberger ratio (Q), is greater than 1 in most cases. Magnetic anomalies can be related to magnetization contrasts among subsurface structures, which are temperature-dependent, since magnetic minerals lose their properties at temperatures above their Curie point. The primary mineral responsible for magnetic anomalies in volcanic environments is titanomagnetite (Fe_3-xTi_xO₄, where 0 ≤ x ≤ 1 and x represents the Ti content), which has a Curie temperature that increases as the Ti content decreases, reaching 575 °C for pure magnetite (x = 0) [27]. In active volcanic areas with associated thermal anomalies and geothermal systems, magnetic studies revealed that some of the detected anomalies are linked to the decrease of the rocks’ magnetization, which can be explained not only by the loss of magnetic properties at high temperatures but also by the chemical alteration affecting magnetic minerals and the circulation of hydrothermal fluids within the rocks [28,29,30,31,32].

In November 2012, a land magnetic survey was performed in the Timanfaya volcanic area. In the Islote de Hilario zone, the total field magnetic data were acquired using an Overhauser magnetometer (GSM-19W, GEM Systems Inc., Toronto, ON, Canada with 0.01 nT resolution and 0.2 nT absolute accuracy), following a 75 m-long profile that coincided in the first 50 m with the GPR profile (Figure 2a,d). Three additional profiles (two of them parallel and one orthogonal to the first) completed the survey. Data sampling along the profiles was performed every 80–100 cm. The magnetic sensor was located approximately 2 m above the ground.

Data in Islote de Hilario area were measured on November 26 (2012) between 4:41 p.m. and 5:12 p.m. local time (GMT-1). Magnetic data acquired at the Güímar magnetic observatory in Tenerife on that date show that the variation in the total field intensity (F) within the same time interval was less than 2 nT. This variation range was so minimal compared with the variations of the crustal field (several hundred nT) that correction of the temporal magnetic variations of external origin was unnecessary.

Magnetic Data Processing and Rock Magnetic Properties

Processing the magnetic data consisted of: (a) subtraction of the main Earth magnetic field modelled using the International Geomagnetic Reference Field [33]; (b) data interpolation into a regular grid with a cell size of 50 cm (Figure 6a); and (c) reduction to the pole assuming a magnetization of induced origin (parallel to the Earth’s magnetic field, with declination of –4.7° and inclination of 38.2°) (Figure 6b). This transformation removes the dipolar character of magnetic anomalies and makes their interpretation easier because the maximum of the reduced-to-the-pole anomaly is located roughly over the center of the causative source [34]. In a reduced-to-the-pole map, magnetic lows are associated with negative magnetization/susceptibility contrasts, whereas magnetic highs are related to positive magnetization/susceptibility contrasts in the subsurface beneath the anomaly.

A paleomagnetic study was conducted on Lanzarote in 2013 [35]. In the Timanfaya volcanic area, three representative basaltic lava flows were sampled: one pertaining to the 1824 Chinero eruption (TM1) and three corresponding to the 1730–1736 eruptions (TM2 to TM4). The magnetic properties relevant to the interpretation of magnetic anomalies are summarized in

Table 2. Summary of the rock magnetic properties of the lava flows in the Timanfaya area.

Site and Lithology	Magnetic Susceptibility (SI)	Induced Magnetization (A/m)	Remanent Magnetization, NRM				Q
			Intensity (A/m)	Declination (°)	Inclination (°)	α95
TM1 (lava flow 1)	0.011 ± 0.007	0.35 ± 0.22	22.5 ± 14.8	331.8	63.4	24.3	64
TM2 (lava flow 2)	0.019 ± 0.018	0.58 ± 0.57	10.4 ± 4.1	346.0	52.9	11.9	18
TM3 (lava flow 3)	0.020 ± 0.014	0.62 ± 0.42	15.2 ± 9.3	2.8	59.2	8.9	24
TM4 (lava flow 4)	0.025 ± 0.013	0.78 ± 0.40	12.9 ± 1.5	321.4	62.3	22.4	17

(personal communication). Laboratory experiments revealed that most of the magnetic signal is carried by titanomagnetite with low to medium Ti content and Curie temperatures ranging between approximately 420 °C and 575 °C [35].

Lava flows in the Timanfaya volcanic area reveal very intense natural remanent magnetizations (NRMs). Induced magnetization was calculated as the product of the magnetic susceptibility and magnetic field in Lanzarote (estimated as 31 A/m), and resulted in being much weaker. Both magnetizations present large variability within the same lava flow (see the high standard deviations). The Köenigsberger ratio (Q) reveals that the NRM is the most relevant component of the total magnetization vector in lava.

In the Islote de Hilario area, pyroclastic fall deposits have the same composition as the Timanfaya lava flows (basalts). Therefore, the expected differences between lava and pyroclastic deposits with regard to their magnetic response are: (1) the magnetization of the pyroclastic layer is only of an induced origin because the remanent component (acquired during the rapid cooling of pyroclasts while falling) is chaotically distributed within the layer, so the summation of the remanent magnetization vectors is null; and (2) the induced magnetization of the pyroclastic layer may exceed the induced magnetization of the sampled lava flows, which were quite vesicular and consequently less dense than the pyroclastic deposits.

3. Results

3.1. GPR Results

Figure 4 presents the reflection profiles obtained in May 2012 (A) and April 2015 (B). The depth axis in radargrams was obtained using the velocity of 0.08 m/ns provided by the CMP-velocity analysis (Figure 3). A closer look at the radargram obtained in 2012 shows that no continuous subhorizontal reflections are displayed up to approximately 8 m depth, and the signatures of the signal reflections patterns vary both laterally and vertically. The areas showing maximum signal reflections (or scattering) in depth agree with the location of the higher temperatures at the surface. Moreover, the reflections of high amplitude display a convex-upwards geometry that is clearly seen, especially below the well-known geothermal anomaly located at the beginning of the profile. The strength of the GPR signal reflections diminishes upwards. Thus, a subsurface area with higher temperatures can be outlined at the beginning of the radargram (from 0 to 35 m), indicating that the heat source is shallower at the beginning of the profile and extends laterally and progressively deeper towards the end. With respect to the radargram obtained in 2015, a similar signature is observed and the subsurface area with reflections of higher strength extends from 4 to 35 m, although the greatest signal strength is concentrated at 4 to 15 m. These subtle variations of the GPR signal reflection patterns suggest a certain temporal variation of the ground temperature, which will be discussed later.

3.2. EMI Results

The EMI1 profile (Figure 5, top panel) shows four high resistivity areas at the surface, which are located at the beginning of the profile (0 to 3 m, TA1), at the middle part (11 to 15 m, TA2, and 17 to 22 m, TA3) and at the end (30 to 40 m, TA4). These high resistivity anomaly zones continue vertically at a depth of, at least, 40 m. The first high-resistive zone (TA1) is located at a well-known circular geothermal anomaly area, whereas the last zone (TA4) is located close to some old monitoring wells used to measure ground temperature, reporting temperature values as high as 388 °C at 7.5 m depth [10]. Considering this, the intermediate high-resistive areas must also be related to geothermal anomalies (TA2 and TA3). Thus, the four high-resistive anomaly areas are associated with high ground temperature zones at the surface that extend vertically downwards through the pyroclastic deposits. The lowest resistivity anomaly values (WA1 and WA2 in Figure 5), placed at both the beginning of the profile and at 28–30 m, are located ~15 m below the surface, but they do not continue at depth. These low-resistive anomaly values agree with the location of man-made walls (Figure 2b,d) built using local basaltic lava rocks as bricks. Considering all this, we suggest that the occurrence of fractures in the underlying volcanic rocks act as conduits for the emission of gases at high temperature, and the surface geothermal anomaly areas are located immediately above these fractures. The high temperature values would modify the electromagnetic properties of the pyroclastic rocks, resulting in an increase in their resistivity values. The significance of the localized low resistivity zones WA1 and WA2 is that they are impacted by the effect the man-made walls have on the coils when they are very close to the walls; therefore, in this case, they can be considered an artefact of anthropogenic origin. Moreover, the intermediate resistivity values located at the surface between TA1-TA2 and TA3-TA4 represent pyroclastic deposits with a ground temperature that is lower than the geothermal anomaly areas.

The EMI2 profile (Figure 5, lower left panel) has a clearer image of the geothermal anomaly area TA1. A well-defined high-resistive area extends from the beginning up to 9 m, which is in good agreement with the dimensions of the geothermal anomaly area at the surface (Figure 2a–c). The high-resistive area that corresponds to TA1 in EMI1 is bowl-shaped and extends downwards to ~15 m depth. Close to its central part (2–3 m from the beginning of the profile), a medium resistivity area can be seen extending vertically up to, at least, 45 m depth. A wide area of lower resistivity values surrounds the high-resistive area both horizontally and vertically. As for the EMI1 profile, good correlation exists between the zone of high-resistive values and the location of the geothermal anomaly area (TA1), confirming the relationship between a high ground temperature and the increased resistivity. The narrow vertical area of medium resistivity values corresponds to the location of the vertically ascending hot gases column emanating from a fault located below the pyroclastic deposits.

The EMI3 profile (Figure 5, lower right panel) runs tangentially to the geothermal anomaly area. For reference, the location of the nearby geothermal anomaly zone has been projected. Low to medium resistivity values dominate most of the shallower parts of the profile (up to ~5–7 m depth); however, for profiles EMI1 and EMI2, high-resistive areas are imaged at greater depths extending downwards to ~45 m depth. The first area is located 1–2 m from the beginning of the profile and corresponds to TA2 in the EMI1 profile; the second area at 8–10 m is immediately beneath the projection of the centre of the geothermal anomaly (TA1 in the EMI1 and EMI2 profiles). The third area (TA5) is located at 17–19 m. Considering the previous profiles, the high-resistive areas correspond to the location of hot gases emanating from the faults buried by pyroclastic deposits, which would also extend laterally close to the surface.

3.3. Magnetic Modelling Results

As previously mentioned, the shallow subsurface structure of the studied zone consists of a pyroclastic deposit that is several tens of metres (>70 m) thick overlying the older basaltic lava flows. From a magnetic point of view, it is expected that the chaotic distribution of lapilli fragments cancels the effect of the thermoremanent magnetization acquired during cooling. Consequently, we can assume the magnetic signal in this survey (caused by the shallowest deposit) is only due to the induced magnetization of the pyroclastic blanket, which is proportional to the magnetic susceptibility and parallel to Earth’s present magnetic field. Within this context, magnetic anomalies can be associated with variations in the magnetic susceptibility produced by high temperatures and the eventual alteration caused by the circulation of high temperature gases within the rocks. In fact, the reduced-to-the-pole anomaly map displays some magnetic lows (Figure 6b), which we interpret as the effect of high temperatures at shallow depths causing a partial or total demagnetization of the subsurface rocks.

Based on this hypothesis, we performed a 2.75D forward modelling (Figure 7) with Geosoft GM-SYS software (version 8.1, Geosoft Inc., Toronto, Canada) along the two orthogonal profiles shown in Figure 6: profile ab, which coincides with both the GPR and EMI1 profiles; and profile cd, which runs parallel to the electromagnetic profile EMI2 (see location in Figure 2a). Profile cd was selected for modelling considering that it cuts the “magnetic low 1” symmetrically. Both profile data (ab and cd) were obtained by “slicing” the reduced-to-the-pole magnetic grid (Figure 6b).

We assumed that the pyroclastic layer was characterized by a uniform induced magnetization of 2 A/m. The three magnetic lows identified in the reduced-to-the-pole magnetic anomaly map (Figure 6b) were modelled as totally (source of low 1) or partially (sources of lows 2 and 3) demagnetized vertical bodies (lower magnetization than the surrounding pyroclastic deposits) located beneath each magnetic low. We hypothesized that the source of magnetic low 1, which is the most intense and is located inside the circular ground thermal anomaly area, is completely demagnetized. The vertical extension of these demagnetized volumes cannot be established due to the size of the survey (magnetization variations below a depth greater than approximately 20 m are not detected), but we can assume that they have a deep source. The main magnetic high was modelled, including the presence of massive surface lava blocks used to build a wall located very close to the profile (WA2). The magnetic highs located at the beginning of both profiles can also be related to the presence of another man-made wall (WA1) (Figure 2b).

4. Discussion

4.1. Relation between GPR Signal and Ground Temperature

The response of a material to the propagation of electromagnetic waves, and consequently the signal reflection and its strength depends on three factors (e.g., [16,36]): the relative permittivity (ε) of the materials through which the electromagnetic waves pass, the conductivity (σ) and the magnetic permeability (µ). The relative permittivity depends on several factors, primarily the nature of the material and the water content. For most of the GPR studies, the greater the change of material and/or water content, the stronger the reflections obtained. However, there is also a thermal effect resulting in a dependence of the relative permittivity upon the temperature [15,16]. In general, ε tends to increase with decreasing temperature. Although for most practical applications of GPR, this temperature dependence is negligible; it can be important when strong variations of temperature are present in the materials surveyed, and this is the case for the geothermal areas in this study. Similarly, σ primarily depends on the nature of the material and the water content (essentially the charges of dissolved anions and cations), whereas µ depends on the ferromagnetic minerals present in the materials (primarily magnetite, maghemite and hematite). Most of the materials do not contain ferromagnetic minerals; thus, their magnetic effect has little effect on the propagation of the GPR signal [37], which is not the case for volcanic materials because their ferromagnetic mineral content is relevant and has a considerable effect on the GPR wave velocity and signal attenuation [16]. As previously mentioned, the magnetic properties of volcanic minerals depend on temperature, exhibiting a decrease in magnetization at high temperatures. This temperature dependence on the value of µ is very significant for the GPR signal response in geothermal areas.

Because the material composition of Islote de Hilario zone is uniform up to at least 70 m depth, as demonstrated by the loose pyroclastic rocks with similar characteristics found from the drilling of several wells, the differences in the strength of the GPR reflections can only be due to differences in water content and/or temperature. However, the water content is extremely low, both at the surface (the landscape is typical of an arid environment, with a mean annual accumulated rainfall of ~150 mm (source: AEMET, Spanish Meteorological Agency)) and at great depths (0.0003–0.002% at the bottom of the monitoring wells [38]). Therefore, temperature is the only factor that can explain the differences in the strength of the signal reflections observed in the radargrams.

By considering the temperature dependence of both ε and µ, the strongest GPR reflections will correspond to the areas where a strong difference in ground temperature exists. This relationship is confirmed at the beginning of the radargram (Figure 4a,b), with the strongest reflections located immediately below the well-known geothermal anomaly area. The convex-upwards geometry of the reflections is interpreted as the geometry of the overheated gases column emanating from the underlying fractured lava flows and propagating towards the surface. Strong GPR reflections, although of lower strength, can also be seen at the central part of the radargram, indicating that the presence of gases at high temperatures is also present below the surface but with a lower contrast of temperature. The geometry of the reflections, decreasing slightly and shallower towards the end of the profile, seems to be the consequence of a combination of an upwards and lateral movement as the gases extend from the emission point, with a consequent decrease in temperature due to the adiabatic expansion. Thus, in geothermal areas of homogeneous materials and low water content, GPR results are a fast and accurate method for locating the distribution at depth of the higher temperature areas inside the geothermal site.

4.2. Seasonal Temperature Variations from the GPR Signal

Once it is determined that variations in the GPR signatures are related to ground temperature changes at the Islote de Hilario zone, with the strongest reflections associated with the greater temperature contrasts, it becomes possible to study if seasonal variations of temperature exist in the area by repeating the GPR survey at different epochs. Ref. [39] has linked the variations in ground temperature and deformation for the Timanfaya volcanic area with changes in surface thermal anomalies produced by the movement of hydrothermal fluids, and changes due to fluctuations in the water table, which is altered by rainfall periods. Thus, it would be interesting to determine if GPR can detect such seasonal variations of ground temperature. Consequently, the GPR survey performed on May 2012 (Figure 4a,b) at the Islote de Hilario zone was repeated in April 2015 (Figure 4c,d). The same processing sequence was selected to ensure that differences in the GPR signatures (if any) were due to changes on the material properties and not artefacts derived from a different data processing. Both radargrams look very similar at a first glance, although some differences are present: the zone showing greater GPR signal reflections, or scattering, is wider for the May 2012 profile, extending along the first 35 m of the profile and from ~2 m depth to the bottom of the profile, becoming deeper towards the end of the radargram (from ~10 m to the bottom). Looking at the April 2015 profile, instead of a single wide zone of high amplitude signatures, three narrower ones can be seen: one at the beginning of the profile and extending from ~2–4 m depth to the bottom; a second one at the middle of the radargram (25–30 m) and a third one near the end of the profile at ~12 m depth. The three areas can be interpreted as three different zones of overheated gas emissions that could potentially be located over underlying fractures affecting the lava flows located at a depth greater than 70 m. A greater volume of overheated gas emissions (and, consequently, a higher ground temperature) would permit the connection of the three temperature anomaly areas laterally converging in a single wider zone of high GPR signal strength, as observed in May 2012. If the volume of gas emissions decreases, the three areas would laterally disconnect, as seen in the April 2015 profile.

Therefore, a certain seasonal variation of the ground temperature can be confirmed from the GPR surveys, which constitutes a novel result considering that the application of GPR to active geothermal areas is very limited in the literature, primarily researching the location of vents, fractures and sinter deposits [40,41,42] but not to the study of the ground temperature distribution and its temporal variations.

4.3. Relation between Resistivity Anomalies and Ground Temperature

Electromagnetic methods have been extensively used in the exploration of geothermal systems ([1,43,44,45] among others). Low resistivity anomalies associated with ascending hydrothermal brines can easily be detected with these methods, and a clear relationship between low resistivity anomalies and high temperature zones can be established ([44], and references herein) which is the case for different hydrothermal areas that have been extensively studied, such as Iceland [43,46], New Zealand [47,48,49], Guadeloupe Island [50], and Italy [51]. In the Canary Islands, a geophysical and geochemical study conducted at Tenerife Island revealed a prominent low-resistivity structure, interpreted as a clay alteration cap related to the mechanism of the upward motion of deep-seated gases from a volcano-geothermal system [52]. However, the application of electromagnetic systems to HDR geothermal systems (such as the Timanfaya volcanic zone) is very limited in the literature.

Our study has revealed a positive correlation between high-resistive anomaly areas determined from EMI measurements and zones of high ground temperatures. This correlation is very clear at the beginning of profiles EMI1 and EMI2 (Figure 5), where bowl-shaped high (>5000 ohm·m) resistivity anomaly areas are located on a well-known geothermal anomaly of the Islote de Hilario zone, extending downwards to, at least, 45 m depth. This correlation is the opposite of that observed in classical geothermal systems associated with hydrothermal convection. That is, the upper part of a HDR geothermal system is characterized by very low liquid water content, with vapour being the dominant phase. As previously mentioned, earlier studies on the composition of emanating gases in the Islote de Hilario zone have revealed that the composition is primarily N2 (95–98%) with very low values of CO2 (0.5–1%) and water vapour (0.0003–0.0020%). Whereas the presence of liquid water greatly decreases the resistivity value of rocks, the occurrence of overheated gases (the temperature values obtained for the Islote de Hilario zone range from 250 to 605 °C) in the absence of a liquid phase dramatically increases the resistivity, and thus a positive correlation between resistivity and ground temperature is found. This relationship has been described in previous works. For example, Ref. [51] in a geophysical study of the geothermal area of Travale (Italy) states that ‘depending on the predominant source of the fluids (liquid or vapour), the fluids produce low- or high-resistivity anomalies.’ Similarly, Ref. [53], in a study conducted at the geothermal site of The Geysers (California), concludes that, if the temperature is above boiling in a borehole, a local zone of high-resistive values could indicate the presence of a steam-filled fracture.

Considering this, the high-resistive anomalies displayed in the EMI profiles correspond to high temperature areas due to the emission of overheated gases from underlying fractures. Along the longest profile, EMI1, up to four different zones of gas emission can be obtained (TA1 to TA4) corresponding to, at least, the location of four fractures affecting the underlying lava flows located below the maximum depth of investigation. Profile EMI2 provides a more detailed image of the TA1 anomaly with a narrow (~2 m) high-resistive subvertical zone that corresponds to the vertical ascent of the overheated gases through the pyroclastic deposit. Additionally, a shallower but wide (~10 m) high-resistive area, which is in good agreement with the extension of the geothermal anomaly observed at the surface, constitutes the lateral expansion of the gases reaching the surface. It must be considered that the resolution of the resistivity model obtained from the inversion of the EMI data decreases with depth, and thus the high temperature areas are much better defined near the ground surface. The two low-resistive anomalies WA1 and WA2 observed in the profile EMI1 correspond to the effect of two man-made walls built from basaltic rock, and these have very different resistivity and magnetic properties than those of the pyroclastic deposits. Accordingly, these anomalies are not related to the temperature variations of the geothermal anomaly areas, which constitutes one limitation of the EMI method because, when anthropogenic structures are very close to the transmitter or receiver coil, the effect can be remarkable, generating resistivity anomalies that are not related to any geological feature. Moreover, due to the specific data acquisition and inversion procedures of the method, these anomaly values horizontally displace half the intercoil spacing and at a depth much deeper than its actual position, which is evident for anomalies WA1 and WA2 that are placed 5 m (i.e., half the minimum intercoil spacing of 10 m) to the right of their actual position and at a mean depth of 15 m (the depth related to the inversion method to the EMI data obtained for both horizontal and vertical coplanar dipole modes for the 10 m intercoil spacing).

All of the above suggest that the EMI method constitutes a useful tool to obtain reliable information about the location of fractures where the emission of overheated gases occurs, as well as the lateral and vertical extension of the rocks affected by the diffusion of the gases and the associated high temperature areas. Although it is very difficult to find examples that apply this method for the exploration of HDR geothermal systems, we consider it a very valuable geophysical prospection technique, as it provides clear images of the temperature distribution up to several tens of metres depth.

4.4. Relation between Magnetic Anomalies and Ground Temperature

The results from the magnetic modelling (Figure 7) have allowed the locations of some narrow areas of demagnetization associated with the presence of magnetic lows to be identified. Up to three different totally or partially demagnetized volumes, ranging in width from 1.5 to 4 m, are located at 6 m (source of the “magnetic low 1”), 13 m (source of the “magnetic low 2”) and 19 m (source of the “magnetic low 3”) of the horizontal distance in profile ab. The source of the “magnetic low 1” was also modelled in the orthogonal direction along profile cd, which can be compared with the TA1 anomaly in the resistivity profile EMI2 (Figure 5). The first demagnetized area (located at 6 m) agrees well with the location of the geothermal anomaly TA1 at the surface. Because magnetic susceptibility is dependent on temperature, it is straightforward to interpret the demagnetized volumes as being discrete high temperature geothermal anomalies in the subsoil. Temperatures measured in the different wells of the Islote de Hilario zone ranged from 250 to 605 °C [10,11]. Considering that Curie temperatures of the Timanfaya basalts vary from approximately 420 to 575°C [35], it becomes evident that the magnetic lows may be associated with the presence of volumes affected by the demagnetization of thermal origin within the pyroclastic deposits. This interpretation is found in previous works where magnetic anomalies were used as exploration tools in geothermal areas i.e., [28,29,30,31,32,54,55]. For example, Ref. [54] linked the negative magnetic anomalies in New Zealand’s high temperature geothermal systems with the demagnetization of the host volcanic rocks. In a similar manner but using aeromagnetic data, Ref. [55] interpreted a negative magnetic anomaly located at the Ngatamariki geothermal field (New Zealand) as demagnetized rocks extending to >1 km depth. A similar interpretation of the negative magnetic anomalies can be found in a more regional study of the high temperature reservoirs of the Taupo volcanic zone (New Zealand), where the authors [30] linked the demagnetization of the volcanic rocks with the alteration of magnetite, caused by geothermal systems expressing both liquid and vapour dominant phases.

The positive correlation between high temperature areas and rocks’ demagnetization constitutes a valuable and simple tool to locate discrete geothermal anomaly areas because of the fast acquisition of the magnetic data. However, in volcanic areas with a more complex geology, the interpretation of the magnetic lows is not so simple and additional information coming from independent geophysical techniques and/or geological data would be needed to provide reliable results.

4.5. Joint Interpretation of GPR, EMI and Magnetic Data for the Detection of Shallow Geothermal Anomaly Areas

Our results have provided the geophysical signature of shallow, high temperature areas inside a HDR geothermal area using three independent geophysical methods. A positive correlation with the temperature exits between the strength of the GPR signal and the high-resistive areas, whereas, for the magnetic anomalies, the correlation is negative (Figure 8). Up to five geothermal anomaly areas have been outlined, and three of them (TA1, TA2 and TA3) have been detected using the three methods. Only TA4, which is associated with a high GPR signal strength and high-resistive values, cannot be correlated to a magnetic low due to the absence of magnetic data at that zone. The magnetic profile deviates by a few metres from the GPR and EMI profiles, precisely where TA4 is located (see Figure 2a). It must also be noted that the data acquisition of the three geophysical methods was conducted at different times (2012, 2013 and 2015). As a seasonal variation of the subsurface, ground temperature has been revealed by means of GPR data, and differences between the locations of the geothermal anomalies revealed by the three techniques may exist.

Due to the ambiguity inherent to geophysical methods, by using only one single technique, some alternative interpretations might be possible. Special care must be taken in areas of complex geology because the different electromagnetic, resistivity and magnetic anomalies might be associated with causes other than temperature variations. For example, Ref. [56], in an area of the Timanfaya volcanic zone very close to the Islote de Hilario zone, has detected hidden lava tubes by using similar GPR, EMI and gravity anomalies that consequently are not related to changes in ground temperature. Moreover, some anomalies have an anthropogenic origin, such as those corresponding to man-made walls that are present in the EMI and magnetic profiles of the present study and that could be misinterpreted as areas of lower ground temperature if additional information is not available. Thus, the best way to reduce any ambiguity is to combine as many independent geophysical methods as possible. Accordingly, the most accurate interpretation when looking for discrete areas of high temperature inside HDR geothermal systems is by referencing the combined occurrence of a high strength GPR signal, a high resistivity anomaly and magnetic low.

Using a combination of different geophysical methods for the exploration of volcanic areas is frequent in the literature. For example, the hydrothermal system of Solfatara (Italy) has been extensively surveyed using resistivity and self-potential methods [57] as well as resistivity and gravity [58]. A combination of magnetotelluric-time domain electromagnetics surveys with resistivity, gravity and magnetics has been used in the geothermal fields of New Zealand [59]. Recently, Ref. [60] surveyed the Piton de la Fournaise volcano (La Réunion Island, France) using self-potential, resistivity, GPR and microgravimetric methods to unravel the air convection and associated thermal field inside the volcano. These works are only a sample among many others that can be found in the literature that primarily focus on the characterization of geothermal and volcanic systems up to several hundreds or even thousands of metres depth. However, we propose here a novel combination of shallow, high-resolution geophysical methods that provide reliable and accurate information for the shallowest part of the geothermal areas. The combination of GPR, EMI and magnetics allows a good compromise between resolution and depth for the investigation of the first several tens of metres below the surface.

5. Conclusions

The shallower part of the HDR geothermal system of the Islote de Hilario, inside the Timanfaya volcanic zone, has been studied using a combination of three independent geophysical techniques: GPR, EMI and magnetic prospecting. The survey of a discrete, well-known geothermal anomaly in the area has been used to depict the geophysical signature of the high-temperature zones, which was then used to infer the location of similar unknown zones distributed around the area. The results from the GPR revealed that, when the material is homogeneous, the signature of the reflections is more intense in the areas with high temperature values, due to the temperature dependence of the physical parameters that determine the response of the rocks to the electromagnetic waves. Similarly, a variation in the subsurface distribution of the thermal anomaly has been detected for the first time through the analysis of GPR at different periods. The resistivity models obtained from the inversion of EMI data demonstrated that high-resistive areas are associated with high temperature zones, due to the increase in resistivity experienced by the volcanic materials when overheated gases propagate through them towards the surface. Regarding the magnetic method, the zones with high temperature values are associated with magnetic lows due to the demagnetization of the volcanic materials when they are heated to temperatures close to or higher than the Curie point of the involved magnetic minerals.

We propose here that a combination of the GPR, EMI and magnetic prospecting methods provides the best results for identifying the location of high ground temperature areas in HDR geothermal systems. Because data acquisition is fast and the methods are non-destructive, they constitute a very useful tool in areas that require special protection, as in the case of Timanfaya National Park. If these techniques can be performed at the same place but at different periods of the year, they can also provide valuable information regarding the seasonality of temperature variations, which can be related to changes in the volume of the gas emissions that cause geothermal anomalies.