3-D interpretation of short-period magnetotelluric data at Furnas Volcano, Azores Islands

3-D interpretation of short-period magnetotelluric data at Furnas Volcano, Azores Islands Summary Accurate geophysical imaging of shallow subsurface features provides crucial constraints on understanding the dynamics of volcanic systems. At Furnas Volcano (Azores), intense circulation of volcanic fluids at depth leading to high CO2 outgassing and flank destabilization poses considerable threat to the local population. Presented is a novel 3-D electrical resistivity model developed from 39 magnetotelluric soundings that images the hydrothermal system of the Furnas Volcano to a depth of 1 km. The resistivity model images two conductive zones, one at 100 m and another at 500 m depth, separated by a resistive layer. The shallow conductor has conductivity less than 1 S m−1, which can be explained by clay mineral surface conduction with a mass fraction of at least 20 per cent smectite. The deeper conductor extends across the majority of the survey area. This deeper conductor is located at depths where smectite is generally replaced by chlorite and we interpret it as aqueous fluids near the boiling point and infer temperatures of at least 240 °C. The less conductive layer found between these conductors is probably steam-dominated, and coincides within the mixed-clay zone found in many volcanic hydrothermal systems. Electrical properties, Hydrothermal systems, Magnetotellurics, Numerical modeling 1 INTRODUCTION The Azores archipelago (Portugal) is formed by nine volcanic islands located in the North Atlantic Ocean where the American, Eurasian and African plates meet at a triple junction (Buforn & Udías 2010). São Miguel Island is the largest (≈65 × 15 km2) of the archipelago and is located over the Terceira Rift, which extends from the island of Santa Maria towards the Mid-Atlantic Ridge following approximately a NW–SE trend (Machado 1959). São Miguel Island comprises three trachytic polygenetic volcanoes: Sete Cidades, Fogo (Água de Pau) and Furnas (Fig. 1). Figure 1. View largeDownload slide The Azores islands are located at the triple junction between the North American, Eurasian and African plates. São Miguel Island has three active strato-volcanoes, Sete Cidades, Fogo (Água de Pau) and Furnas. The red frame outlines the survey area of this work around the inner caldera at Furnas. Modified from Guest et al. (1999) and Wallenstein et al. (2007). Figure 1. View largeDownload slide The Azores islands are located at the triple junction between the North American, Eurasian and African plates. São Miguel Island has three active strato-volcanoes, Sete Cidades, Fogo (Água de Pau) and Furnas. The red frame outlines the survey area of this work around the inner caldera at Furnas. Modified from Guest et al. (1999) and Wallenstein et al. (2007). Furnas, a volcanic complex, is the eastern-most of the three trachytic active central volcanoes of São Miguel Island. It does not have a well-developed edifice and is a superposition of two calderas. The external and largest caldera was formed about 30 ka BP, followed by a period of caldera filling eruptions. The inner caldera is smaller and was formed around 10–12 ka BP, which is responsible for the steep topography of more than 500 m in our study area. It contains several younger eruptive centres, a shallow caldera lake and several zones of intensive fumarolic activity within the inhabited caldera (Guest et al. 1999, 2015). Phlegrean fields volcano (Naples, Italy) may be comparable to Furnas volcano in several aspects, such as the existence of a caldera and the composition of the hydrothermal fumarolic emissions. The last magmatic eruption in this volcanic system took place in 1630 AD (Cole et al. 1995) causing about 200 fatalities. Present-day activity comprises secondary manifestations of volcanism that are characterized by boiling temperature fumaroles, steaming ground, thermal springs, cold CO2 rich springs and areas of intense diffuse degassing especially in the northern part of the inner caldera, north of Furnas Lake, where detailed maps of surface CO2 and 222Rn emissions were recently made available (Viveiros et al. 2010; Silva et al. 2015). CO2 emissions can pose serious threats to populations that reside in volcanic regions. This is for instance the case at Furnas, where 60-90 per cent of dwellings are built over hydrothermal CO2 emission zones (Viveiros et al. 2010). In addition, hydrothermal explosions were reported in one of the fumaroles located close to Furnas village (Pedone et al. 2015, and references therein). Along with this risk, intense circulation of volcanic fluids at depth, coupled with alteration can lead to associated mechanical weakening of the volcanic edifice and may lead to caldera wall instabilities. This process may be triggered and/or coupled with elevated rainfall. The evaluation and forecasting of volcanic risks require the understanding of the geometry of the hydrothermal circulation system and quantifying the relationship between the surface degassing features and deeper processes driving the heat and fluid flows. The complex tectonic structure of Furnas is affected by fault systems trending WNW–ESE, NE–SW, N–S and NNW–SSE, expressed by volcanic alignments, linear valleys and caldera outlines (Carmo et al. 2015). The caldera walls have orientations that coincide with the main fault trends (WNW–ESE and NE–SW), suggesting that caldera collapses were controlled by structural weaknesses as proposed by Guest et al. (1999) and Carmo et al. (2015). The tectonic regime of the area is thus complex; because of the low amount of local seismicity, the recent state of tectonic stress is not well constrained. Even less is known on the existence, location, and geometry of the magmatic source, though the consensus is that there is a magmatic reservoir region, constrained to shallow crustal depths of 3 to 6 km beneath the Furnas Caldera by tectonics (Machado 1959), GPS (Sigmundsson et al. 1995), magnetic data (Blanco et al. 1997), gravity data (Camacho et al. 1997; Montesinos et al. 1999), seismic tomography (Zandomeneghi et al. 2008) and geochemistry (Jeffery 2016). The exact location and structure of the magmatic reservoir is somewhat unknown, with indications of a position south of the caldera centre (Zandomeneghi et al. 2008, and references therein). High electrical conductivity values are typically associated with volcanic-hydrothermal systems and the elevated conductivity structure can provide constraints on these volcanic and hydrothermal processes (Moore et al. 2008; Revil et al. 2010). These settings are a good target for electrical geophysical methods that are sensitive to conductivity. The resistivity (the reciprocal of conductivity) of a geological unit is a function of the amount and interconnectivity of a conductive fluid contained within the matrix of the lithology and also the surface conductivity. In general, bulk resistivity decreases with increases in pore water connectivity, saturation, salinity, temperature, and clay content due to rock alteration (e.g. Revil et al. 1998; Ussher et al. 2000; Revil et al. 2002; Muñoz 2014). Studies to date have found that conductive areas correspond to the pathway and phase states of hydrothermal fluids, and to the fluid bearing structure in volcanic and geothermal areas (e.g. Ogawa et al. 1998; Aizawa et al. 2005, 2009; Revil et al. 2011). Thus, resistivity anomalies can provide important constraints leading to an improved understanding on structure and dynamics of hydrothermal systems. The first magnetotelluric (MT) measurements in the Furnas Caldera date from 1983, when the United States Geological Survey (USGS) carried out geophysical surveys for assessing geothermal potential on the three volcanoes on São Miguel Island (Hoover et al. 1983, 1984). Though using only scalar Audio-MagnetoTelluric (AMT) soundings, they were able to identify zones of high geothermal potential, in particular, the geothermal reservoir at the NW flank of the Fogo (Água de Pau) volcano, which now hosts two geothermal power plants (Franco 2015) providing more than a third of the São Miguel's energy demand. Aiming at an improved imaging of subsurface electrical conductivity, and a better understanding of the shallow, fumarole-related features within the Furnas hydrothermal system, 39 AMT including 15 Broad-Band MagnetoTelluric (BBMT) systems were deployed in 2015 and 2016. In this study, we focus on the AMT portion of the data in order to identify the main features of the Furnas hydrothermal system. The use of the remaining long-period data requires accurate knowledge of the bathymetry around the island and will be presented in a follow-up paper. 2 PREVIOUS STUDIES AT FURNAS VOLCANO Furnas volcano is situated on the eastern part of the island and its earliest lithology is dated from about 100 000 yr BP (Moore 1990). It comprises a multicyclic forming caldera which last erupted in 1630 AD. Furnas exhibits a 5 × 8 km summit depression formed by juxtaposed calderas controlled by NW–SE and NE–SW faults. Conjugate faulting, parallel to the Terceira Rift regional fault system are observed and exhibit a NW–SE trend (Guest et al. 1999). The present caldera system is flanked by steep walls as high as 500 m, that formed following multi-collapse episodes and subsequent caldera-filling eruptions. The younger caldera contains several craters, and the western part is covered by a lake that occupies an area of 1.9 km2 (Cruz & França 2006). Recent volcanic activity at Furnas includes two historic events since the island became inhabited in the 15th century. These two intracaldera volcanic eruptions have occurred at Furnas volcano in 1439–1443 (which led to the creation of Gaspar Dome also referred to as Pico do Gaspar) and in 1630 (referred to as the ‘1630 Dome’) that led to the formation of two tuff rings with central domes that occupy the central part of the caldera floor as shown in Fig. 2. A detailed discussion on the composition, stratigraphy and eruptive history of Furnas can be found in Guest et al. (1999, 2015, and references therein). At present, volcanic activity of Furnas volcano is characterized by the secondary manifestations of volcanism mostly located in the northern region of the caldera, within Furnas village and north of the Furnas Lake, and are characterized by boiling temperature fumaroles, steaming ground, thermal springs, cold CO2-rich springs and areas of intense diffuse degassing, namely CO2 and radon (Viveiros et al. 2010; Caliro et al. 2015; Silva et al. 2015). Figure 2. View largeDownload slide Locations of MT stations around the Furnas volcano. The red inverted triangles represent AMT only sites, while the blue ones represent combined AMT and BBMT sites. The structural features are also mapped: solid line (fault), dashed line (inferred fault) and dotted line (volcanic alignment), courtesy of Carmo et al. (2015). The main fumarole area at Furnas Lake is marked by the yellow circle. The pink circle marks the location of site FUR007A (Fig. 3). Elevation data provided by the IVAR GIS Group at the University of the Azores. Figure 2. View largeDownload slide Locations of MT stations around the Furnas volcano. The red inverted triangles represent AMT only sites, while the blue ones represent combined AMT and BBMT sites. The structural features are also mapped: solid line (fault), dashed line (inferred fault) and dotted line (volcanic alignment), courtesy of Carmo et al. (2015). The main fumarole area at Furnas Lake is marked by the yellow circle. The pink circle marks the location of site FUR007A (Fig. 3). Elevation data provided by the IVAR GIS Group at the University of the Azores. The subsurface structure of the Furnas volcanic-hydrothermal system is not well constrained. In 1985, two aeromagnetic surveys were made covering Faial and São Miguel Islands. The results, first published by Miranda et al. (1991), yielded high amplitude positive magnetic anomalies, coherent with the existence of highly magnetized and young eruptive rocks covering the flanks of most of the volcanic systems on the island, however two calderas, Sete Cidades and Furnas match relative magnetic lows, whereas Fogo (Água de Pau) corresponds to a large negative magnetic anomaly. Miranda et al. (2015) inferred the decrease in magnetization at Furnas to be the result of hydrothermal alteration, similar to the conclusion derived from Blanco et al. (1997) for the decrease in the magnetic field anomaly at Furnas, as sensed by land based magnetic surveys. The intense degassing and fumarolic features at Furnas (Viveiros et al. 2010; Silva et al. 2015) suggest the presence of such a mature hydrothermal system beneath the central volcano and supports an equivalent interpretation. As mentioned, Blanco et al. (1997) disclosed the existence of magnetic anomalies within the Furnas area. Forjaz (1996) documents that this alteration would affect the shallowest materials more intensively within the caldera that are derived from different structures in and around the volcano. The inner part of the caldera, around Furnas Lake and the historic volcanic domes which are most relevant to our study, coincides with a broad negative magnetic anomaly. The authors suggest that the cause may be related to the oxidation of the magnetic minerals as a result of hydrothermal alteration from the surface to a depth of 500 m, where the hydrothermal system seems most active, however its effect probably exists at much greater depths. The magnetic anomaly map derived from this work is only sensitive to structures to a depth of approximately 3 km, thus the proposed magma chamber, presumed to be at a depth of at least 5 km, could not be imaged. Geodetic studies, between 1991 and 1994, across central São Miguel revealed a slight inflation in the area of Furnas, with the centre of this suggested inflation in the northern part of the caldera. Sigmundsson et al. (1995) infers this phenomenon being caused by growing pressures of fluids within a deep hydrothermal system beneath the caldera which is sealed by an impermeable zone. It is worth noting that the authors suggest that the centre of deformation is located within the zone of the negative anomaly detected by Blanco et al. (1997). Camacho et al. (1997) and Montesinos et al. (1999) analysed and modelled gravity data collected around Furnas, yielding gravity anomaly maps that isolated low-density zones that may be the result of several distinct caldera collapses and infilling. One such anomaly, a principal gravity minimum is located below Pico do Gaspar, just east of Furnas Lake, and is quite distinct at a depth of 1000 m from the work of Montesinos et al. (1999). This structure is probably based on the formation of the inner caldera that cuts the old caldera followed by another phase of infilling. Similar low-gravity anomalies related to low-density products of explosive volcanic activity and caldera collapse processes were reported by Sherburn et al. (2003) for Taupo zone calderas and Vanorio et al. (2005) for Campi Flegrei. A seismic tomography study covering the centre of São Miguel between the Fogo (Água de Pau) and Furnas volcanoes was conducted by Zandomeneghi et al. (2008), identifying a low velocity region beneath Furnas. It is characterized by low P-wave velocities and low Vp/Vs ratios, which decrease with depth until they reach a minimum of 2 km below sea level. The authors’ results coincide spatially with the gravity inversion low-density anomaly. The authors infer the low Vp/Vs ratios changes being related to changes in fluid conditions within an extensive steam-dominated geothermal system. Geothermal exploration and development commenced on São Miguel as early as 1975, however, this exploration was confined to the area on the northern slope of Fogo (Água de Pau). In 1982 and 1983, the USGS conducted reconnaissance AMT surveys of the three quiescent central volcanoes on the island, to evaluate the potential for geothermal systems (Hoover et al. 1983, 1984). While their survey yielded reasonable results, it was strongly limited by the narrow frequency range and use of a scalar system, that is, only measuring the electric (E) and corresponding orthogonal magnetic (H) field in an assumed fixed coordinate system (NW–SE and NE–SW). Thus only two of the four elements of the impedance tensor (eq. 1, Section 3) were determined. Nevertheless, their results at Furnas indicated an east-west trending, low resistivity zone through the centre of the caldera. This zone correlates with the known fumaroles and hot springs at Furnas Lake and at Furnas village. The observed apparent resistivities at this zone were estimated to be 10 Ωm from their apparent resistivity maps at 7.5 Hz (this was the lowest frequency they recorded at which signal strengths were adequate to give good data quality) in the near surface and drop to approximately 1.5 Ωm at a depth of 50 m from 1-D estimates. Similar low apparent resistivities were observed to the south and east of the caldera near the coast, however the authors suggest that sea water intrusion on these data cannot be ruled out because of the proximity of the Atlantic Ocean, however it should also be noted that the studies of Viveiros et al. (2010) and Silva et al. (2015) also identify fumarolic fields and soil diffuse degassing anomalous areas (CO2 and 222Rn) at the coast south of the Furnas Caldera near Ribeira Quente, concluding that localized geothermal systems are present outside the caldera. 3 MT METHOD The MT method is a passive geophysical exploration technique that yields information on the 3-D spatial distribution of electrical conductivity of the Earth's subsurface by measuring and relating the natural time-varying horizontal electric (E) and magnetic (H) fields at its surface (Tikhonov 1950; Cagniard 1953). The relationship between horizontal and mutually perpendicular field components provides amplitude responses (apparent resistivity) and phase responses at various frequencies, known as MT response functions, for each site recorded. Amplitudes of the fields decrease exponentially with increasing depth in a uniform conductor, known as the skin-depth phenomenon. The depth of penetration of the MT signal into the earth depends on both the conductivity of the rocks and the frequency of the variation of the MT signal; with low-frequency (or long-period) variations sensing deeper. The conductivity structure of the subsurface can be determined from the amplitude and phase relationships between the surface electric and magnetic field vectors as a function of period and sounding location. These relationships are expressed mathematically by the impedance tensor Z, where Z (Vozoff 1972), is a complex 2 × 2 matrix as   \begin{equation} (E_x\ E_y) = \left( {\begin{array}{*{20}{c}}{Z_{xx}}&{{Z_{yx}}} \\ {{Z_{xy}}}&{{Z_{yy}}} \end{array}} \right)\left( {\begin{array}{*{20}{c}}{H_x} \\ {{H_y}} \end{array}} \right). \end{equation} (1) From the elements of the impedance tensor, the more intuitive apparent resistivities and phases can be calculated as   \begin{equation} {\rho _{a,ij}} = \frac{{{\mu _0}}}{\omega }{\left| {{Z_{ij}}(\omega )} \right|^2},{{\rm }}{\phi _{ij}} = \arctan \frac{{\Im ({Z_{ij}})}}{{\Re ({Z_{ij}})}} \end{equation} (2)where μ0 is the free-space magnetic permeability. The angular frequency ω is defined 2πf, while ℜ and ℑ represent the real and imaginary parts, respectively. A further quantity calculated from the MT data is the local magnetic transfer function (also complex) between the horizontal and vertical magnetic field components, T, commonly called induction vector. For this study the sign of this vector was chosen so that the induction vector points in the direction of increasing conductivity (Parkinson 1962):   \begin{equation} \left({\begin{array}{*{20}{c}}{H_z} \end{array}} \right) = \left( {\begin{array}{*{20}{c}}{T_{x}}&{{T_{y}}} \end{array}} \right)\left( {\begin{array}{*{20}{c}}{H_x} \\ {{H_y}} \end{array}} \right). \end{equation} (3) All of these quantities provide information on the subsurface distribution of resistivity. In particular, the structure of the impedance tensor and its rotational invariants can be used to determine whether data may be interpreted in a 1-D or 2-D model, or require 3-D modelling. Further details can be found in recent comprehensive monographs (such as Simpson & Bahr 2005; Chave & Jones 2012). 4 FIELD DATA AND PROCESSING Following a pilot survey of 13 AMT soundings and Electrical Resistivity Tomography (ERT) data collected along two profiles in the eastern part of Furnas Caldera in 2015, a second campaign was completed in June 2016, yielding a total of 39 soundings including 15 BBMT soundings, using the Phoenix Geophysics MTU-5A data-logger along with MTC-30 and MTC-50 induction coils for AMT and BBMT magnetic field data, respectively. Due to the steep topography within the caldera, obtaining a site distribution in a regular grid proved impossible, however good coverage was obtained within and outside the caldera (Fig. 2). The survey area was free of major noise sources. At each site, the two horizontal, perpendicular magnetic field components Hx and Hy were recorded as well as the vertical magnetic field component Hz. The two horizontal, perpendicular electric field components Ex and Ey were measured using non-polarizing Phoenix Geophysics Pb–PbCl (lead–lead chloride) electrodes (Petiau & Dupis 1980), laid out in a cross with a dipole length of typically 80 m. The time-series were recorded overnight and have a duration of approximately 12 hr. Time-series processing techniques applied to the data have yielded apparent resistivity and phase response curves for each site. The acquired MT data were processed from the raw time-series via commercial remote-reference software from Phoenix Geophysics (based on Jones & Jödicke 1984). Remote reference processing was computed with concurrent sites and also with the remote site, however processed data from the remote site did not always improve the final response significantly for the MT sites. The best response (i.e. the smoothest transfer function with the best responses in the low signal frequency band between 1 and 5 kHz) from the far and local remote site was selected as the final response of the individual site. Both magnetic field and electric field references were processed, however the magnetic field data were used as it gave superior responses. Following processing, the physical consistency of the resistivity and phase estimates at each station were checked using the D+ algorithm (Parker 1980; Parker & Whaler 1981; Beamish & Travassos 1992; Parker & Booker 1996), using the WinGlink® software package. Outliers and inconsistent data points were removed from the data set, prior to data analysis and inversion. Site locations with both AMT and BBMT data had their responses merged to yield one response for the site. An example of a high-frequency MT response can be seen in Fig. 3. With the full five-component measurements now available, it was possible to analyse the dimensionality and directionality of the observations with different methods used within the MT community. Figure 3. View largeDownload slide Response curves of an example AMT site. Apparent resistivity (ρa, top) curves, phase curves (ϕ, middle), and induction arrows (bottom) in Parkinson convention are shown. The site location is marked with a pink circle on Fig. 2. Figure 3. View largeDownload slide Response curves of an example AMT site. Apparent resistivity (ρa, top) curves, phase curves (ϕ, middle), and induction arrows (bottom) in Parkinson convention are shown. The site location is marked with a pink circle on Fig. 2. 5 DATA ANALYSIS Prior to analysing and modelling of the acquired MT data, the presence of the Atlantic Ocean around the survey area was investigated to quantify the sharp electrical contrast between the ocean and land. The skin-depth of a typical BBMT sounding can easily reach great depths (several tens of km) depending on the actual conductivity structure, thus the ocean has a significant impact on the observed responses when the separation distance between the coast and the survey area is smaller than the depth of interest. This makes it difficult to determine reliable geoelectrical structures at deeper regions of the subsurface. In order to quantify the impact of the Atlantic Ocean, forward modelling was used to create synthetic responses from a pre-determined 3-D resistivity model. A simple 3-D model of São Miguel Island with a seawater resistivity value of 0.3 Ωm extending to a depth of 6000 m and a conservative uniform land resistivity of 100 Ωm was designed. The results from this forward modelling, indicate that the ocean does not impact on data with periods less than 1 s, therefore for the purpose of this study, only data up to a period of 1 s are used in the analysis and modelling. Analysis of the geoelectric dimensionality of MT data has acquired relevance due to the advancement of numerical codes that have made it possible, with confidence, to model and invert the data using 1-D, 2-D or 3-D approaches. Simultaneously, one can gain insights into the direction of the geoelectrical strike at various depths, which in turn can be used to qualitatively correlate with different geological processes and subsurface geological structures. A dimensionality study focussing on the rotational invariants of the MT tensor was presented by Weaver et al. (2000), however, in real and noisy data, the invariants are rarely zero, yielding geoelectric dimensionality of noisy responses to be interpreted as 3-D. Weaver et al. (2000) overcame this difficulty by implementing a threshold value, τ, beneath which the invariants are set to zero. The work of Martí et al. (2005) and Martí et al. (2010), that led to the development of WALDIM (Martí et al. 2009), a FORTRAN application to perform dimensionality analysis of MT data, found that a value of τ between 0.1 and 0.2 determines the dimensionality very well, and that the geoelectric strike direction can be estimated, when 2-D structure is inferred through a Monte Carlo approach with the addition of Gaussian noise. Along with classifying the dimensionality into bands of periods for each site, the WALDIM application has been extended to identify hints of anisotropy through the assessment of the invariants and strike directions from a pattern of sites and frequencies. According to Martí et al. (2010), these indicators of anisotropy occur in the inferred 2-D cases when there is an inconsistency in strike directions and when relationships between the tensors do not imply isotropic media. In this study, a τ of 0.15 was used and the dimensionality of each site was estimated within each period decade band. The results, along with the estimates of geoelectrical strike direction whenever quasi-2-D cases are found can be seen in Fig. 4. Figure 4. View largeDownload slide Results of dimensionality analysis and geoelectric strike estimations (purple coloured arrows) for four bands using the WALDIM code (Martí et al. 2009), superimposed on the main structural lineaments in the Furnas region as derived from Carmo et al. (2015). Note the dominant change in dimensionality of the geoelectrical subsurface from 1-D (purple squares) at the highest frequency range (10 000 to 1000 Hz) to a more 3-D (green squares) environment at the deepest regions (10 to 1 Hz). Hints of anisotropy are also estimated and can be seen to be coincident with steep topography and adjacent to main structural features. Legend (following the definitions of Martí et al. (2009): 1-D, one dimensional; 2-D, two dimensional; 3-D/2-Dtw, two -dimensional data affected by distortion (only twist); 3-D/2-D, general case of galvanic distortion over a 2-D structure; 3-D, affected or not by galvanic distortion; 3-D/2-Ddrt, 3-D or 2-D with diagonal regional tensor; 3-D/2-D or 3-D/1-D, galvanic distortion over a 1-D or 2-D structure with a non-recoverable strike direction; Anis1, homogeneous anisotropic medium; Anis 2, anisotropic body within a 2-D medium. Figure 4. View largeDownload slide Results of dimensionality analysis and geoelectric strike estimations (purple coloured arrows) for four bands using the WALDIM code (Martí et al. 2009), superimposed on the main structural lineaments in the Furnas region as derived from Carmo et al. (2015). Note the dominant change in dimensionality of the geoelectrical subsurface from 1-D (purple squares) at the highest frequency range (10 000 to 1000 Hz) to a more 3-D (green squares) environment at the deepest regions (10 to 1 Hz). Hints of anisotropy are also estimated and can be seen to be coincident with steep topography and adjacent to main structural features. Legend (following the definitions of Martí et al. (2009): 1-D, one dimensional; 2-D, two dimensional; 3-D/2-Dtw, two -dimensional data affected by distortion (only twist); 3-D/2-D, general case of galvanic distortion over a 2-D structure; 3-D, affected or not by galvanic distortion; 3-D/2-Ddrt, 3-D or 2-D with diagonal regional tensor; 3-D/2-D or 3-D/1-D, galvanic distortion over a 1-D or 2-D structure with a non-recoverable strike direction; Anis1, homogeneous anisotropic medium; Anis 2, anisotropic body within a 2-D medium. The results of the dimensionality analysis indicate that the Furnas Caldera region is quite complex, with dimensionality ranging from 1-D to 3-D for the four frequency bands. At the highest frequency band (10000 Hz – 1000 Hz), 1-D structure is dominant across the whole survey area, complementing the overlapping XY and YX modes in the MT responses. As the longer-period data probes deeper into the subsurface, the geoelectrical structure becomes more 3-D and this can be clearly seen in the 10 Hz – 1 Hz panel. Hints of anisotropy can be seen throughout the data set in all frequency bands, however the majority of these sites appear to be located in regions of steep topography (particularly north of Furnas Lake), and also in structurally controlled regions. The WALDIM code was also used to quantify the geoelectrical strike direction within the data set where the dimensionality showed 2-D structure and where hints of anisotropy were sensed. The general trend of this strike is NE–SW, with the average strike direction for each decade ranging from 42° to 55°. The overall complexity of the subsurface structure beneath Furnas can be further illustrated by plotting maps of the phase tensor ellipses from all sites as a function of period. The MT phase response becomes more sensitive to structure at greater depth and less sensitive to structure at shallow depth as the period increases. In contrast, the magnitude or apparent resistivity response will be distorted by the galvanic effect of near-surface conductivity heterogeneities at all periods. However, galvanic effects do not significantly affect the phase tensor (Caldwell et al. 2004) and phase tensor ellipse maps (Fig. 5) provide a method of visualizing spatial gradients in the conductivity structure directly. The principal axes of the phase tensor ellipse, ϕmax and ϕmin, indicate the horizontal directions of the maximum and minimum induction current, which reflects lateral variations in the resistivity structure. The variation of the orientation of the major axis of the ellipse can help identify changes in the dimensionality of the subsurface structure. In the case of 1-D structure, ϕmax = ϕmin and therefore the ellipse will be a circle. In the 2-D case, ϕmax ≠ϕmin and the phase tensor is an ellipse but also symmetric, however when a 3-D subsurface is probed, the ellipse is not symmetric. Four maps at frequencies that are approximately halfway within each decade band along with induction arrows are plotted in Fig. 5 in which the ellipses are scaled so that the long axis (representing ϕmax) is of uniform length, and are colour coded according to the value of ϕmin. At the highest frequency ranges (centred at 5263 Hz and 532 Hz), a significant amount of the ellipses in the survey area are quite circular, inferring the near surface conductivity structure to be essentially 1-D. This compliments the results from the WALDIM dimensionality study. The only exceptions to these are a number of sites that show a small degree of ellipticity which may be a response to the steep topography. At these regions, the magnitude of ϕmin is relatively low (colder colours) inferring a geoelectrical subsurface that is resistive. The magnitude of ϕmin shows that much higher phase values—an indicator of increasing conductivity (warmer colours) with depth—occur over the whole survey area at longer periods which can be seen in the frequency ranges (centred at 49 and 5.6 Hz). A remarkable conductive area is clearly shown just beneath the northeastern edge of Furnas Lake (Figs 5c and d). Induction arrows can be used to illustrate the relationship between the measured vertical and horizontal magnetic field components. Adopting the convention proposed by Parkinson (1962), the real part of such arrows point towards conductive bodies or concentrations of induced currents. Analysis of the induction arrows confirm that there is a conductive body beneath the region of Furnas Lake (Figs 5b and c), however, at the longest periods (Fig. 5d) the induction arrows provide a strong indication of the presence of a larger-scale conductivity structure located beyond the eastern and the northeastern part of Furnas Lake. Figure 5. View largeDownload slide Maps of phase tensor ellipses and real part of the induction arrows at four frequencies. The colour of the phase tensor ellipses (axes normalized by ϕmax) indicates ϕmin. The induction arrows are plotted using the Parkinson convention. Figure 5. View largeDownload slide Maps of phase tensor ellipses and real part of the induction arrows at four frequencies. The colour of the phase tensor ellipses (axes normalized by ϕmax) indicates ϕmin. The induction arrows are plotted using the Parkinson convention. Following this analysis, it is clear that the conductivity structure associated with the hydrothermal system beneath Furnas is 3-D and the observed MT data need to be inverted accordingly. Furthermore, the steep topography associated with the caldera can only be accounted for by true 3-D inverse modelling. 6 3-D MODELLING As discussed in Section 5, the potential influence of the Atlantic Ocean on the data set is absent for periods ≤1 s. The 3-D inversions thus were performed including only high-resolution topography and Furnas Lake bathymetry data employing the parallel version of the Modular system for Electromagnetic inversion code ModEM (Meqbel 2009; Egbert & Kelbert 2012; Kelbert et al. 2014). The 3-D mesh consisted of 119 × 119 × 149 (plus 10 air layers) in the north–south, east–west and vertical directions, respectively. The centre of the mesh comprised a uniform mesh of 89 × 89 × 149 cells, with a horizontal cell size of 50 m × 50 m. The lateral extents of the padding cells, 15 planes in each of the north, south, east and west directions, increased by a factor of 1.5. In the vertical direction, the cell thickness was of 7.5 m for 698.264 m < z ≤ 817.5 m; subsequent layer thicknesses were increased by a factor of 1.2. The input data were decimated to 17 periods per site with regular distribution for periods from 0.0001 to 1 s. In total, 39 sites were used in the inversion process. The 3-D inversions were run with full impedance tensor elements (Zxx, Zxy, Zyx and Zyy) and vertical magnetic transfer functions. Error floors were set to an absolute value of 3 per cent of (∣Zxy × Zyx∣)1/2. For the vertical transfer functions, a constant error of 0.015 was used. The starting model resistivity was set to 100 Ωm. The high-resolution bathymetry of the lake was included as a priori information and kept fixed during the inversion. A resistivity value of 63 Ωm was used for the lake (Andrade et al. 2016). Model regularization employed a smoothing parameter of 0.4 for both horizontal directions, and 0.3 for the vertical direction. A commonly used metric for the fit of the modelled to the observed data is the normalized root mean square nRMS, defined as   \begin{equation} {\rm nRMS} = \sqrt{\left( {{N_d} - 1} \right)_{}^{ - 1} \sum \nolimits _{j = 1}^{N_d} \left( \frac{d_{{\rm obs},j} - d_{{\rm calc},j}}{\sigma _j} \right)^2 } {\, }, \end{equation} (4)where dobs,j and dcalc,j are observed and calculated data, respectively, and σj represents the data errors for all N data points. A value of this metric close to 1 is commonly interpreted as a data fit within the range of observational error, that is, in the case of this study the error floors, which somewhat reduces the statistical significance. For the final inversion model, 112 iterations were required to fit the full impedance tensor and the tipper data, with an nRMS of 1.8. As shown in Fig. 6, this value is strongly influenced by very few large-residual sites, as is characteristic for a non-robust metric of the least-squares type. The fit to the majority of the observations becomes much clearer from the data fit curves and detailed maps of nRMS, which were calculated for each data components and for each frequency band at each station, are presented in Supporting Information Sections S1 and S2, respectively. The 3-D inversion results are summarized in Figs 7 and 8, and discussed in the following section. Figure 6. View large Download slide Map of nRMS misfit values for each station. A total nRMS of 1.8 was achieved. Maps of variation in nRMS for each frequency band and for all the elements of the impedance tensor are shown in Supporting Information Section S2. Figure 6. View large Download slide Map of nRMS misfit values for each station. A total nRMS of 1.8 was achieved. Maps of variation in nRMS for each frequency band and for all the elements of the impedance tensor are shown in Supporting Information Section S2. Figure 7. View largeDownload slide Horizontal slices of the final resistivity model at various depth. Panels (a) and (b) depict the depth extent of the shallow conductor beneath the fumarolic field on the northern edge of Furnas Lake. Panels (c) and (d) indicate the limits of the main conductor, C2, beneath the inner caldera at Furnas. Figure 7. View largeDownload slide Horizontal slices of the final resistivity model at various depth. Panels (a) and (b) depict the depth extent of the shallow conductor beneath the fumarolic field on the northern edge of Furnas Lake. Panels (c) and (d) indicate the limits of the main conductor, C2, beneath the inner caldera at Furnas. Figure 8. View largeDownload slide Extracted 2-D slices from our final 3-D resistivity model. (a) P1 is an NS profile passing through the main fumarolic region on the northern edge of Furnas Lake towards the most recent site of volcanism, the 1630 Dome, (b) an EW profile through the centre of the survey imaging beneath Furnas Lake and Pico do Gaspar and (c) a profile passing through the centres of the main domes in the inner caldera, Gaspar Dome and the 1630 Dome in the south of the profile. Figure 8. View largeDownload slide Extracted 2-D slices from our final 3-D resistivity model. (a) P1 is an NS profile passing through the main fumarolic region on the northern edge of Furnas Lake towards the most recent site of volcanism, the 1630 Dome, (b) an EW profile through the centre of the survey imaging beneath Furnas Lake and Pico do Gaspar and (c) a profile passing through the centres of the main domes in the inner caldera, Gaspar Dome and the 1630 Dome in the south of the profile. 7 RESULTS Depth slices through the 3-D resistivity inversion volume are shown in Fig. 7 and selected cross-sections through various regions of interest are shown in Fig. 8. The overall structure comprises a resistive surface layer (R1) of varying thickness with resistivities in the region of 100–1000 Ωm. It overlays zones of high conductivity (C1 and C2) beneath the main caldera, ranging from depths of 100 m to a general depth of at approximately 500 m. Almost all of the MT soundings recorded in the campaign are characterized by relatively high near-surface resistivity, in the region of 100–200 Ωm, with the exception of one site (namely FUR004, Supporting Information Section S1), adjacent to the main fumarolic fields, with a resistivity range in the order of tens of Ωm at the near surface. The shallow subsurface of the entire area within this survey (R1) is caldera infill, dominated by material derived from volcanic activity and subsequent erosion. The deposits include coarse lapilli, fine ashes and pumice blocks in the upper deposits, however due to the lack of any borehole data the full extent of this stratigraphy is unclear. The pore space is filled with water of meteoric origin, probably similar to the water within Furnas Lake (≈60 Ωm) (Andrade et al. 2016). A shallow, very localized, conductive zone (C1) is found in the northern region of Furnas Lake, at ≈100 m beneath the surface, adjacent to a mapped fault, that is steeply dipping (≈70°), known as Furnas Fault 2 (FF2, see Carmo 2013; Fig. 5). This conductive zone extends vertically ≈100 m, with a resistivity of 0.5–10 Ωm, which explains the behaviour of the high-frequency induction vectors, indicating a concentration of conductive material at this location and depth (Figs 5b and c). This zone may not be completely controlled by FF2 as the orientation of the lens like structure is not parallel to FF2 strike. Geoelectric strike direction estimates, derived from our WALDIM analysis, infer ranges of 40°–55° (Fig. 4), however, the strike of FF2 is ≈100° (Carmo et al. 2015). Note that the orientation of C1 is almost parallel to the steep caldera wall which may suggest a relation to the limits of caldera infill, or tectonic features related to the caldera wall. We interpret this conductive zone as a shallow, cap-like structure that partially seals the ascent of vapour feeding the surface hydrothermal manifestations directly above, possibly influenced/controlled by the fault, FF2. Fig. 9 compares the resistivity sections obtained by AMT and ERT methods approximately along an NE–SW, 1.5 km-long profile (P4) along with the soil CO2 degassing results from Viveiros et al. (2010). The corresponding ERT profile transects the FF2 structure on the northern edge of Furnas Lake. Data were acquired using a multielectrode ABEM system (Terrameter LS) using a roll-along Wenner-Schlumberger configuration (64 electrodes with a spacing of 10 m) with an injected current of 100–400 mA. The resulting apparent resistivity data including topography information were then inverted using the RES2DINV® software (Loke & Barker 1996; Loke 2016). Both geophysical methods clearly and coherently image the shallow conductive zone C1, that is found directly beneath the highest concentration of CO2 degassing in the region of 400–900 m on P4. This illustrates the consistency of using both electrical and EM geophysical methods in imaging subsurface features and isolating pathways beneath high CO2 degassing output in this part of the caldera. Due to the lack of coverage of data around Furnas village, which also exhibits localized zones of high CO2 degassing, we cannot conclude if there is a distinct network of shallow conductive zones distributed beneath the village feeding other zones of degassing. Figure 9. View largeDownload slide Comparison of AMT and ERT results along a profile P-4 (Fig. 8) across the fumarole area. Top: Satellite image of surface and topography, overlain by the CO2 flux distribution from Viveiros et al. (2010). The dashed white line marks the location of the main fault in the area, FF2. Centre: inversion results for an ERT profile starting at the NW rim of the lake following a roughly NNE direction. Bottom: slice through the 3-D AMT model approximately along the ERT profile above. Both methods agree well in their overlap. Figure 9. View largeDownload slide Comparison of AMT and ERT results along a profile P-4 (Fig. 8) across the fumarole area. Top: Satellite image of surface and topography, overlain by the CO2 flux distribution from Viveiros et al. (2010). The dashed white line marks the location of the main fault in the area, FF2. Centre: inversion results for an ERT profile starting at the NW rim of the lake following a roughly NNE direction. Bottom: slice through the 3-D AMT model approximately along the ERT profile above. Both methods agree well in their overlap. Beneath the localized conductor C1 we find a thin, more resistive region (R2) of approximately 100 m thickness that increases in resistivity with depth. It is approximately at sea level and separates C1 from the dominant conductive zone beneath the survey area, C2. This large conductive zone has an upper boundary at 150 m b.s.l., approximately 500–600 m beneath the surface. It has a very conductive upper layer in the region of 1–10 Ωm, and lies directly beneath C1, however the upper boundary of C2 varies spatially in our inversion, but this may result from poor data coverage especially in the region of Furnas Lake. The highest conductivity found within C2 lies directly beneath C1, which could indicate they are linked, however sensitivity analysis of our 3-D model in this region could not confirm that C1 and C2 are separate zones. Unfortunately, there are no borehole data in the Furnas region so is not possible to confirm or refute this hypothesis. The contrast between C2 and R2 can be clearly seen in 2-D slice, P1, that is orientated north-south along the eastern edge of Furnas Lake through our survey area (Fig. 8a). It shows C2 thinning towards the south and becoming more resistive but interestingly shows no conductive body beneath the 1630 Dome, the remnant of the most recent event of volcanism. C2 is clearly visible beneath Gaspar Dome and can be seen in both Profiles P2 (E–W across the survey area, Fig. 8b) and P3 (passing through the two volcanic domes, Fig. 8c). The top of this anomalous feature is found at a depth of ≈250 m b.s.l., that is, ≈600 m beneath the surface, with a thickness of at least 500 m. We cannot confirm the eastern boundary of C2, however C2 appears to connect the main fumarolic fields at Furnas with Gaspar Dome, but thins towards the south. The resolution of the resulting 3-D model is significantly reduced ≈1 km below the surface, underneath the C2 anomaly, though it cannot be excluded that this is due to the shielding effects of the high-conductivity zones above. The absence of C2 in the south of the survey implies that saline intrusion from the Atlantic Ocean is not being sensed by our high-frequency MT survey at these depths. Earlier modelling studies, prior to analysis and inversion, confirmed the ocean effect to be absent at these high-frequencies. 8 DISCUSSION The electrical conductivity of the subsurface is a crucial parameter in the characterization of volcanic and hydrogeothermal settings and has been exploited widely by MT surveys in various settings around the world in recent decades (Muñoz 2014, and references therein). Regions of high conductivity within the shallow Earth, can be attributed to partial melt related to magma reservoirs, the presence of graphite, metallic material, or conducting fluids. Magma reservoirs, in particular in view of the high melt fraction required, are highly improbable in the depth range found in this study as there is no volcanological evidence pointing to this explanation. Carbon deposits in the form of graphite conductors are equally improbable because the depth ranges of the conductive zones do not favour this explanation, in particular when the required large volume of connected conducting material is considered. Franco (2015) reports disseminated pyrite found in boreholes in a geothermal system at the northern flank of the near-by Fogo (Água de Pau) Volcano, ranging over the whole depth of the geothermal system, though no estimates of their quantity are given. The contribution of pyrite to the conduction process will be minor, but cannot be excluded. For hydrothermal systems, several possible explanations for the high conductivities exist: hydrothermal systems may appear conductive because of the fluids involved (which would need to be highly saline), or the alteration products caused by recent and ancient fluid movements, for example, a conductive clay cap overlaying the fluid system itself. Both may occur side by side in a given hydrothermal system. The resistivity structure determined by our 3-D inversion of high-frequency MT data, maps the structure(s) of the hydrothermal system beneath the Furnas Caldera and we will discuss this hypothesis in detail in the following. Our 3-D inversion of the upper hydrothermal system at Furnas yielded two highly conductive zones C1 and C2, respectively. Though a robust feature of our models, sensitivity studies of locally perturbed models show that we cannot be certain that they are fully separated. The data fit is not substantially changed by introducing a conductive connection between the two regions, however, their representative resistivity characteristics at different depths imply a complex system of pathways, which feeds the fumaroles at the surface above. For the purpose of this discussion we will treat them as separate zones. In volcanic areas, a conductor with resistivities in the region of 0.5–10 Ωm has been attributed to the presence of hydrothermal water or alteration minerals (Ussher et al. 2000; Muñoz 2014), however both involve contrasting permeability structures to be present. A layer of smectite clay typically exists around the reservoir in a geothermal system. Smectite and smectite-illite clays are the predominant products derived from rock alteration processes that form in the 50–200 °C zone over and adjacent to volcanic hosted systems (Browne 1978). These clay-rich zones are important components for geothermal hydrology, quite often partially sealing the reservoir, where high permeability is required for hydrothermal fluids to migrate within the system, usually through an extensive network of fractures. Typically, values of 1–10 Ωm are observed at clay-rich zones, whereas for the adjacent hydrothermal reservoir, conductivity ranges between 5 and 100 Ωm are believed to be typical (Wright et al. 1985; Pellerin et al. 1996). Thus, the hottest regions of a geothermal system are often characterized by higher resistivity than observed in zones of hydrothermal alteration (Ussher et al. 2000; Muñoz 2014). Geothermal manifestations (fumaroles) are found at the surface directly above the high conductivity zones C1 and C2 at the northern edge of Furnas Lake. Recent work by Woitischek et al. (2017), characterizing the origin of hydrothermal waters on São Miguel by their chemical and isotopic composition, confirm that the hydrothermal waters found at Furnas are generally of meteoric origin. At the surface they are slightly saline with resistivities near 10 Ωm, while the resistivity of the water within Furnas Lake is about 60 Ωm (Andrade et al. 2016). It is unknown how deep the meteoric waters migrate and which role they play in the full system. In the case of C1, the presence of clay material (smectite) is required along with the presence of water to account for the high conductivity values, sensed at these depths between 100 and 200 m beneath the surface. Petrophysical reasoning based on the work of Revil et al. (2002) quantifying the electrical properties of altered volcaniclastic materials can be used to estimate the percentage of smectite needed to account for a conductivity value, once an estimate of the temperature is known. Fig. 10 shows the range of smectite required to account for the conductivity of C1 for a broad range of temperatures and fluid conductivities. The petrophysical model assumed is described in Appendix. From the results we conclude that a mass fraction of approximately 20 per cent smectite is sufficient to explain the observed bulk conductivity of C1, where the remaining 80 per cent comprise fluid and the unaltered original volcaniclastic material. Figure 10. View largeDownload slide Calculation of smectite content based on the studies of Waxman & Smits (1968), Revil et al. (1998, 2002) and Flóvenz et al. (2005) for the temperature range 80–200 °C. The conductivity of the shallow conductor C1 (upper grey zone) structure cannot be explained by ionic conduction in the fluid assuming Archie's law (lower grey zone) even when comparatively high temperatures and salinities are used. Taking into account the surface conductivity of smectite allows us to explain the behaviour of the C1 zone assuming a moderate mass fraction of this clay mineral. An explanation of the underlying estimate is given in Appendix. Figure 10. View largeDownload slide Calculation of smectite content based on the studies of Waxman & Smits (1968), Revil et al. (1998, 2002) and Flóvenz et al. (2005) for the temperature range 80–200 °C. The conductivity of the shallow conductor C1 (upper grey zone) structure cannot be explained by ionic conduction in the fluid assuming Archie's law (lower grey zone) even when comparatively high temperatures and salinities are used. Taking into account the surface conductivity of smectite allows us to explain the behaviour of the C1 zone assuming a moderate mass fraction of this clay mineral. An explanation of the underlying estimate is given in Appendix. Unfortunately, there are no local borehole data to constrain our observations, however borehole data from geothermal exploration at the NW flank of Fogo Volcano (Ribeira Grande, approximately 20 km away from our study area) may provide some clues, even though the morphological position and hydrodynamical characteristic of that hydrothermal system differs from Furnas, the geochemical composition was found to be quite similar (Caliro et al. 2015). The work of Franco (2015) at Fogo (Água de Pau), identified two zones of hydrothermal alteration beneath a shallow zone of relatively unaltered rocks, based on borehole results. Evidence from this study documents a smectite layer at comparable depth to the conductive zone C1, that is observed in our study and the author assigns a temperature range 200 °C to 220 °C. This observation complements our petrophysically based estimates that confirm the requirement of smectite to yield the high conductivity of C1. We propose that the conductivity of C1 is controlled by the presence of smectite clay material, forming a partial seal of a localized vapour driven aquifer, at a temperature of up to 200 °C. Near Furnas Lake it is perturbed by the fault zone FF2. The horizontal resistivity slices in Fig. 7 show that at shallow depths, the strike of the anomalous zone deviates from FF2, with a significant clockwise rotation of the axis at larger depth when we reach C2. No low resistivity zones are found along the eastern continuation of this fault zone. This structure is consistent with the idea that the continuous formation of the low-permeability seal zone inhibits direct ascent, with preferential fluid or steam pathways changing with depth. Beneath C1 we find the thin, more resistive zone R2. Sensitivity analysis on the 3-D modelling cannot rule out that both C1 and the deeper, more extensive conductive zone C2 are connected. The resistivity change may be explained by a change in clay alteration mineralogy that can also act as a proxy for changes in temperature. The study of Arnason et al. (2008) documents that in the temperature range from 220 °C to 240–250 °C, low temperature zeolites gradually disappear and smectite is transformed into chlorite, yielding a mixed-clay zone, where smectite and chlorite coexist. Thus, R2 probably represents a mixed clay zone. This zone may exhibit a wide range of resistivities (Arnason et al. 2008; Hersir & Árnason 2015), from 10 Ωm to more than 1000 Ωm. The resistivity is strongly influenced by residual smectite content, less by porosity and pore/fracture characteristics for a low-salinity system as assumed in this study. Weisenberger et al. (2016) recently have shown that surface conduction, which critically depends on the cation exchange capacity (CEC, see Appendix), is considerably lower in chlorite-rich zones. The emergence of C2, the dominant conductive zone beneath Furnas may indicate the top of the chlorite-dominated zone, where all of the smectite has disappeared from the overlying mixed clay zone. Arnason et al. (2008) shows that the chlorite zone commences at 250 °C, and that at higher temperatures 260–270 °C, epidote becomes the dominant mineral in the so-called chlorite-epidote zone, where conductivities may greater than 100 Ωm (Ussher et al. 2000). However, these temperature estimates apply to freshwater systems, while brine fluids would shift the zoning to higher values. We conclude here that the more recent temperature estimates of ≈275 °C from Caliro et al. (2015) for the hydrothermal fluids feeding Furnas Lake fumaroles, are associated with the deeper reservoir C2. Ocean water (0.3 Ωm) can be eliminated as a source of the high conductivity within the system from the geochemical signature found by Woitischek et al. (2017). Furthermore, the observed conductive zone C2 does not extend towards the south, and we do not find a connection to the ocean there. One may speculate that the lack of southward connection originates from the recent explosive volcanism within the caldera, which has been migrating from the NE to SW since at least 1.9 ka ago (Guest et al. 2015) and may have destroyed earlier fluid pathways. It has to be noted, however, that ocean water is not the only potential source of salinity as the fluid associated with the C2 zone may derive from mixing with highly saline magma-derived fluids (Nicholson 1993; Yardley & Bodnar 2014). However, from our measurements it is not possible to further constrain fluid composition. If we assume that electrical conduction occurs mainly through the fluid phase filling the reservoir, we need the hydrothermal system to be in fluid phase. Fig. 11 shows vertical resistivity profiles extracted from the 3-D model for the sites nearest to the fumarole area, together with the boiling curves for different values of salinities, calculated using the method of Haas (1971) as modified by Canet et al. (2011). The position of C2 (centre) in the resistivity profiles at ≈500 m beneath the surface would imply a temperature of 250 °C, independent of the assumed salinity. This agrees well with the thermal structure found at the Ribeira Grande system (Rangel 2014; Franco 2015), where similar temperatures of about 240 °C were found. The recent geothermometric estimate of 275 °C (Caliro et al. 2015) is significantly higher than the value derived here, but still falls into the depth range of the C2 conductor. Figure 11. View largeDownload slide Boiling curves for different salinities between 1 and 8 wt per cent, calculated using the method described in Haas (1971) and Canet et al. (2011) for hydrostatic conditions and near-adiabatic ascent. Ocean water has an average salinity value of 3.5 wt per cent. The orange curves are vertical profiles at the sites nearest to the fumarole area, extracted from the final 3-D modelling results. The vertical blue line marks the recent temperature estimate of ≈270 °C from Caliro et al. (2015). Assuming that the conductor C2 is related to the fluid-steam interface, an estimate of reservoir temperature at about 250 °C could be derived from our models. This agrees notably well with the temperatures found in the geothermal wells in the Cachaços-Lombadas sector of the Ribeira Grande geothermal field at the NW flank of the Fogo (Água de Pau) Volcano. Note that temperatures close to the surface, where mixing with meteoric waters is important, may be considerably lower than the boiling curves suggest. Figure 11. View largeDownload slide Boiling curves for different salinities between 1 and 8 wt per cent, calculated using the method described in Haas (1971) and Canet et al. (2011) for hydrostatic conditions and near-adiabatic ascent. Ocean water has an average salinity value of 3.5 wt per cent. The orange curves are vertical profiles at the sites nearest to the fumarole area, extracted from the final 3-D modelling results. The vertical blue line marks the recent temperature estimate of ≈270 °C from Caliro et al. (2015). Assuming that the conductor C2 is related to the fluid-steam interface, an estimate of reservoir temperature at about 250 °C could be derived from our models. This agrees notably well with the temperatures found in the geothermal wells in the Cachaços-Lombadas sector of the Ribeira Grande geothermal field at the NW flank of the Fogo (Água de Pau) Volcano. Note that temperatures close to the surface, where mixing with meteoric waters is important, may be considerably lower than the boiling curves suggest. The complex conductivity image of the Furnas Caldera can possibly be explained by the following hypothetical process: an extended hydrothermal system within a large caldera develops, leading to the development of the typical associated alteration zones (clay minerals). Within this caldera, a succession of smaller eruptions occur (e.g. Gaspar and 1630 Domes), which perturb the pre-existing hydrothermal system locally, and add to the sedimentary filling of the original caldera. Subsequently, a new generation of hydrothermal systems develop, destroying or replacing the older ones. As pointed out by Guest et al. (2015), the area of the Furnas Caldera has a documented history of frequent eruptions since the formation of the Povoção Caldera at least 30 ka ago. 9 CONCLUSIONS This MT study has yielded, for the first time, a 3-D electrical resistivity model of the shallow hydrothermal system beneath Furnas Volcano (Azores, Portugal) and has been used to create a revised conceptual model of this dynamic area. The inversion of high frequency data identified highly conductive zones (0.5–10 Ωm) beneath the main central caldera, that feed the surface manifestations of volcano-hydrothermal activity that are characterized by boiling temperature fumaroles, steaming ground, thermal springs, and areas of intense diffuse degassing of CO2 and 222Rn. As no borehole data or temperature profile for this area is available, our petrophysical reasoning (see Appendix) led to the conclusion that the shallowest conductive zone (≈100 m beneath the surface) comprise at least 20 per cent of smectite. This high clay fraction could partially seal the deeper reservoir comprising meteoric water feeding the fumaroles on the surface above. New insights into the structural regime controlling this shallow conductive aquifer have been achieved. The geoelectrical strike of the shallow conductor is not aligned with the strike of the main fault in the area, which was previously interpreted to directly control the orientation of the feeding aquifers. The geoelectrical strike of the conductor is parallel to the steep caldera wall adjacent to it, which may infer that morphology of the caldera, especially beneath the infilling from various phases of caldera collapses and subsequent erosion is a factor in its structural character, however the main fault acts as a pathway from this aquifer towards the surface. The deeper, more dominant conductor imaged by MT outlines the broader extent of the Furnas hydrothermal system within the caldera. The extent of this conductive zone has not been fully mapped by our study due to access issues adjacent to Furnas village and the presence of a caldera lake. This conductor however does not extend towards the south, implying it is not associated with salt water intrusion from the Atlantic Ocean, confirming the recent geochemical analysis of the waters at Furnas to be mainly meteoric in origin. This conductive zone extends beneath Gaspar Dome but is not evident beneath the most recent centre of volcanism, the 1630 Dome, implying the current hydrothermal system is a result of various volcanic events. Borehole and thermal gradient data from geothermal exploration studies at Fogo (Água de Pau), an adjacent volcano, provide clues in characterizing the hydrothermal system in our study. We infer the deeper conductive zone beneath Furnas is dominated by chlorite, inferring a temperature of at least 250 °C. Imaging the proposed magma reservoir that fuels the volcano-hydrothermal system at Furnas, is not within the scope of this study. The inclusion of longer period MT data, combining topography and bathymetry data in the inversion, will be required to successfully image this deeper structure. ACKNOWLEDGEMENTS This multiyear study was primarily financed by DIAS. We appreciate the funding granted to AR (ISTerre) through the INSU project ‘Improving imaging of hydrothermal systems by joint inversion of ERT-AMT data. Case study of Furnas Caldera (Azores, Portugal)’. The fieldwork would not have been possible without the help of Ernesto Sousa, Victor Sousa, Bruno Medeiros, Lucía Rodríguez, Joana Eleutério, Annika Rödder, Damien Jougnot and Niels Grobbe. The Direção Regional do Ambiente, Serviços de Recursos Hídricos e Ordenamento do Território (Ponta Delgada) is thanked for providing the high-resolution bathymetry for the Furnas Lake. The tectonic maps were produced with the data from Azores Regional Government through a PhD grant to RC from Fundo Regional da Ciência e Tecnologia, and through Serviço Regional de Proteção Civil e Bombeiros dos Açores in the scope of the scientific and technical protocols to guarantee the Azores Seismovolcanic Surveillance and the Emergency Planning Studies of Centro de Informação e Vigilância Sismovulcânica dos Açores (CIVISA). The Irish Centre for High Performance Computing (ICHEC) is thanked for availing the Fionn cluster to carry out the numerical computations. Finally, G. Egbert, A. Kelbert, and N. Meqbel are gratefully thanked for making their ModEM code available to us. We thank Léa Levy for her comments on temperature stability range of secondary minerals, in particular smectite. Last, but not least, we thank the two reviewers, J. Hübert and J. Peacock, for their constructive and helpful comments. REFERENCES Aizawa K. et al.  , 2005. Hydrothermal system beneath Mt. Fuji volcano inferred from magnetotellurics and electric self-potential, Earth planet. Sci. Lett.  235( 1-2), 343– 355. https://doi.org/10.1016/j.epsl.2005.03.023 Google Scholar CrossRef Search ADS   Aizawa K., Ogawa Y., Ishido T., 2009. Groundwater flow and hydrothermal systems within volcanic edifices: Delineation by electric self-potential and magnetotellurics, J. geophys. Res.  114( B1), B01208. https://doi.org/10.1029/2008JB005910 Andrade C., Viveiros F., Cruz J.V., Coutinho R., Silva C., 2016. Estimation of the CO2 flux from Furnas volcanic lake (Sao Miguel, Azores), J. Volcanol. Geotherm. Res.  315 51– 64. https://doi.org/10.1016/j.jvolgeores.2016.02.005 Google Scholar CrossRef Search ADS   Archie G.E., 1942. The electrical resistivity log as an aid in determining some reservoir characteristics, Trans. Am. Inst. Min. Metall. Pet. Eng.  146 54– 62. Arnason K., Karlsdottir R., Eysteinsson H., Flovenz O.G., Gudlaugsson S.T., 2008. The resistivity structure of high-temperature geothermal systems in iceland, in Short Course III on Exploration for Geothermal Resources, Lake Naivasha, Kenya, October 24 - November 17  UNU-GTP, (http://www.unugtp.is, last accessed 15 September 2017). Beamish D., Travassos J., 1992. The use of the D+ solution in magnetotelluric interpretation, J. Appl. Geophys.  29( 1), 1– 19. https://doi.org/10.1016/0926-9851(92)90009-A Google Scholar CrossRef Search ADS   Blanco I., García A., Torta J.M., 1997. Magnetic study of the Furnas volcano (Azores), Ann. Geofis.  XL 341– 358. Browne P.R.L., 1978. Hydrothermal alteration in active geothermal fields, Annu. Rev. Earth Planet. Sci.  6 229– 250. https://doi.org/10.1146/annurev.ea.06.050178.001305 Google Scholar CrossRef Search ADS   Buforn E., Udías A., 2010. Azores-Tunisia, a tectonically complex plate boundary, in Advances in Geophysics  vol. 52, chap. 3, pp. 139– 182, ed. Dmowska R., Elsevier. Google Scholar CrossRef Search ADS   Cagniard L., 1953. Basic theory of the magneto-telluric method of geophysical prospecting, Geophysics  18 605– 635. https://doi.org/10.1190/1.1437915 Google Scholar CrossRef Search ADS   Caldwell T.G., Bibby H.M., Brown C., 2004. The magnetotelluric phase tensor, Geophys. J. Int.  158 457– 469. https://doi.org/10.1111/j.1365-246X.2004.02281.x Google Scholar CrossRef Search ADS   Caliro S., Viveiros F., Chiodini G., Ferreira T., 2015. Gas geochemistry of hydrothermal fluids of the S. Miguel and Terceira Islands, Azores, Geochim. Cosmochim. Acta  168 43– 57. https://doi.org/10.1016/j.gca.2015.07.009 Google Scholar CrossRef Search ADS   Camacho A.G., Montesinos F.G., Vieira R., 1997. A three-dimensional gravity inversion applied to Sao Miguel Island (Azores), J. geophys. Res.  102 7717– 7730. https://doi.org/10.1029/96JB03667 Google Scholar CrossRef Search ADS   Canet C., Franco S.I., Prol-Ledesma R.M., González-Partida E., Villanueva-Estrada R.E., 2011. A model of boiling for fluid inclusion studies: application to the Bolaños Ag-Au-Pb-Zn epithermal deposit, Western Mexico, J. Geochem. Explor.  110 118– 125. https://doi.org/10.1016/j.gexplo.2011.04.005 Google Scholar CrossRef Search ADS   Carmo R., 2013. Estudos de neotectónica na ilha de S. Miguel, uma contribuição para o estudo do risco sísmico no arquipélago dos Açores., PhD thesis  University of the Azores. Carmo R., Madeira J., Ferreira T., Queiroz G., Hipólito A., 2015. Volcano-tectonic structures of São Miguel Island, Azores, in Volcanic Geology of S. Miguel Island (Azores Archipelago)  chap. 6, pp. 65– 86, eds Gaspar J.L., Guest J.E., Duncan A.M., Barriga F., Chester D.K., Geological Society of London. Google Scholar CrossRef Search ADS   Carvalho M., Forjaz V., Almeida C., 2006. Chemical composition of deep hydrothermal fluids in the Ribeira Grande geothermal field (São Miguel, Azores), J. Volcanol. Geotherm. Res.  156( 1-2), 116– 134. https://doi.org/10.1016/j.jvolgeores.2006.03.015 Google Scholar CrossRef Search ADS   Chave A., Jones A., 2012. The Magnetotelluric Method: Theory and Practice  Cambridge Univ. Press. Google Scholar CrossRef Search ADS   Cole P., Queiroz G., Wallenstein N., Gaspar J., Duncan A., Guest J., 1995. An historic subplinian/phreatomagmatic eruption: the 1630 AD eruption of Furnas volcano, São Miguel, Azores, J. Volcanol. Geotherm. Res.  69( 1–2), 117– 135. https://doi.org/10.1016/0377-0273(95)00033-X Google Scholar CrossRef Search ADS   Cruz J.V., França Z., 2006. Hydrogeochemistry of thermal and mineral water springs of the Azores archipelago (Portugal), J. Volcanol. Geotherm. Res.  151( 4), 382– 398. https://doi.org/10.1016/j.jvolgeores.2005.09.001 Google Scholar CrossRef Search ADS   Egbert G.D., Kelbert A., 2012. Computational recipes for electromagnetic inverse problems, Geophys. J. Int.  189( 1), 251– 267. https://doi.org/10.1111/j.1365-246X.2011.05347.x Google Scholar CrossRef Search ADS   Flóvenz O.G., Spangenberg E., Kulenkampff J., Árnason K., Karlsdóttir R., Huenges E., 2005. The role of electrical interface conduction in geothermal exploration, in Proceedings World Geothermal Congress , Antalya, Turkey, 24– 29 April 2005. Forjaz V., 1996. Geothermal reservoirs of Furnas Laboratory Volcano, in 2nd Workshop on European Laboratory Volcanoes  Santorini, Greece. Franco A., 2015. Subsurface geology and hydrothermal alteration of Cachaços-Lombadas sector, Ribeira Grande geothermal field, (São Miguel island, Azores), Tech. Rep. UNU-GTP-2015-10, United Nations University, Geothermal Training Programme, (http://www.unugtp.is, last accessed 15 September 2017). Guest J., Pacheco J.M., Cole P.D., Duncan A.M., Wallenstein N., Queiroz G., Gaspar J.L., Ferreira T., 2015. The volcanic history of Furnas Volcano, São Miguel, Azores, in Volcanic Geology of São Miguel Island (Azores Archipelago)  chap. 6, pp. 125– 134, eds Gaspar J.L., Guest J.E., Duncan A.M., Barriga F., Chester D.K., Geological Society of London. Google Scholar CrossRef Search ADS   Guest J.E., Gaspar J.L., Cole P.D., Queiroz G., Duncan A.M., Wallenstein N., Ferreira T., Pacheco J.-M., 1999. Volcanic geology of Furnas Volcano, São Miguel, Azores, J. Volcanol. Geotherm. Res.  92( 1–2), 1– 29. https://doi.org/10.1016/S0377-0273(99)00064-5 Google Scholar CrossRef Search ADS   Haas J.L., 1971. The effect of salinity on the maximum thermal gradient of a hydrothermal system at hydrostatic pressure, Econ. Geol.  66 940– 946. https://doi.org/10.2113/gsecongeo.66.6.940 Google Scholar CrossRef Search ADS   Hersir G., Árnason K., 2015. Resistivity of rocks, in Short Course X on Exploration for Geothermal Resources  organized by UNU-GTP, GDC and KenGen, at Lake Bogoria and Lake Naivasha, Kenya, 2015 November 9–December 1, (http://www.unugtp.is, last accessed 15 September 2017). Hoover D., Amaral R., Broker M., 1983. Preliminary report on audio-magnetotelluric survey on São Miguel island, Azores, Portugal., Open-File Report 82-441, USGS, Menlo Park CA, USA. Hoover D., Rodrigues Da Silva A., Pierce H.A., Amaral R., 1984. Application of audio-magnetotelluric surveys on São Miguel island, Azores, Portugal, in Transactions of the Geothermal Research Council , no. 8, pp. 499– 503, Menlo Park CA, USA. Jeffery A.J., 2016. Petrogenesis and contrasting eruption styles of peralkaline silicic magmas from Terceira and São Miguel, Azores, PhD thesis , Keele University, UK. Jones A., Jödicke H., 1984. Magnetotelluric transfer function estimation improvement by a coherence-based rejection technique, in SEG Technical Program Expanded Abstracts 1984 , pp. 51– 55. Kelbert A., Meqbel N., Egbert G.D., Tandon K., 2014. ModEM: A modular system for inversion of electromagnetic geophysical data, Comput. Geosci.  66 440– 53. https://doi.org/10.1016/j.cageo.2014.01.010 Google Scholar CrossRef Search ADS   Loke M., 2016. ‘Tutorial: 2-D and 3-Delectrical imaging surveys’. Available at: http://www.geotomosoft.com, last accessed 15 September 2017. Loke M., Barker R.D., 1996. Rapid least-squares inversion of apparent resistivity pseudosections by a quasi-Newton method, Geophys. Prospect.  44( 1), 131– 152. https://doi.org/10.1111/j.1365-2478.1996.tb00142.x Google Scholar CrossRef Search ADS   Machado F., 1959. Submarine pits of the Azores plateau, Bull. Volcanol.  21 109– 116. https://doi.org/10.1007/BF02596510 Google Scholar CrossRef Search ADS   Martí A., Queralt P., Jones A.G., Ledo J., 2005. Improving Bahr's invariant parameters using the WAL approach, Geophys. J. Int.  163 38– 41. https://doi.org/10.1111/j.1365-246X.2005.02748.x Google Scholar CrossRef Search ADS   Martí A., Queralt P., Ledo J., 2009. WALDIM: A code for the dimensionality analysis of magnetotelluric data using the rotational invariants of the magnetotelluric tensor, Comput. Geosci.  35( 12), 2295– 2303. https://doi.org/10.1016/j.cageo.2009.03.004 Google Scholar CrossRef Search ADS   Martí A., Queralt P., Ledo J., Farquharson C., 2010. Dimensionality imprint of electrical anisotropy in magnetotelluric responses, Phys. Earth planet. Inter.  182( 3-4), 139– 151. https://doi.org/10.1016/j.pepi.2010.07.007 Google Scholar CrossRef Search ADS   Meqbel N.M., 2009. The electrical conductivity structure of the Dead Sea Basin derived from 2D and 3D inversion of magnetotelluric data, PhD thesis , Freie Universität Berlin, Germany. Miranda J.M., Luis J.F., Abreu I., Victor L.A.M., Galdeano A., Rossignol J.C., 1991. Tectonic framework of the Azores Triple Junction, Geophys. Res. Lett.  18( 8), 1421– 1424. https://doi.org/10.1029/91GL01607 Google Scholar CrossRef Search ADS   Miranda J.M., Luis J.F., Lourenço N., Fernandes R.M.S., 2015. The structure of the Azores triple junction: implications from São Miguel Island, in Volcanic Geology of São Miguel Island (Azores Archipelago), chap. 2 , pp. 5– 13, eds Gaspar J.L., Guest J.E., Duncan A.M., Barriga F., Chester D.K., Geological Society of London. Google Scholar CrossRef Search ADS   Montesinos F.G., Camacho A.G., Vieira R., 1999. Analysis of gravimetric anomalies in Furnas caldera (São Miguel, Azores), J. Volcanol. Geotherm. Res.  92( 1-2), 67– 81. https://doi.org/10.1016/S0377-0273(99)00068-2 Google Scholar CrossRef Search ADS   Moore J.N., Allis R.G., Necmčok M., Powell T.S., Bruton C.J., Wannamaker P.E., Raharjo I.B., Normanc D.I., 2008. The evolution of volcano-hosted geothermal systems based on deep wells from Karaha-Telaga Bodas, Indonesia, Am. J. Sci.  308 1– 48. https://doi.org/10.2475/01.2008.01 Google Scholar CrossRef Search ADS   Moore R.B., 1990. Volcanic geology and eruption frequency, São Miguel, Azores, Bull. Volcanol.  52 602– 614. https://doi.org/10.1007/BF00301211 Google Scholar CrossRef Search ADS   Muñoz G., 2014. Exploring for geothermal resources with electromagnetic methods, Surv. Geophys.  35( 1), 101– 122. https://doi.org/10.1007/s10712-013-9236-0 Google Scholar CrossRef Search ADS   Nicholson K., 1993. Geothermal Fluids  Springer-Verlag, Berlin Heidelberg, ISBN: 978-3-642-77846-9. Google Scholar CrossRef Search ADS   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. https://doi.org/10.1186/BF03352120 Google Scholar CrossRef Search ADS   Parker R.L., 1980. The inverse problem of electromagnetic induction: existence and construction of solutions based on incomplete data, J. geophys. Res.  85( 8), 4421– 4428. https://doi.org/10.1029/JB085iB08p04421 Google Scholar CrossRef Search ADS   Parker R.L., Booker J.R., 1996. Optimal one-dimensional inversion and bounding of magnetotelluric apparent resistivity and phase measurements, Phys. Earth planet. Inter.  98( 3-4), 269– 282. https://doi.org/10.1016/S0031-9201(96)03191-3 Google Scholar CrossRef Search ADS   Parker R.L., Whaler K., 1981. Numerical methods for establishing solutions to the inverse problem of electromagnetic induction, J. geophys. Res.  86 9574– 9584. https://doi.org/10.1029/JB086iB10p09574 Google Scholar CrossRef Search ADS   Parkinson W., 1962. The influence of continents and oceans on geomagnetic variations, Geophys. J. R. astr. Soc.  6 441– 449. https://doi.org/10.1111/j.1365-246X.1962.tb02992.x Google Scholar CrossRef Search ADS   Pedone M., Viveiros F., Aiuppa A., Giudice G., Grassa F., Gagliano A.L., Francofonte V., Ferreira T., 2015. Total (fumarolic + diffuse soil) CO2 output from Furnas volcano, Earth Planets Space  67 doi:10.1186/s40623-015-0345-5. https://doi.org/10.1186/s40623-015-0345-5 Google Scholar CrossRef Search ADS   Pellerin L., Johnston J.M., Hohmann G.W., 1996. A numerical evaluation of electromagnetic methods in geothermal exploration, Geophysics  61( 01), 121– 130. https://doi.org/10.1190/1.1443931 Google Scholar CrossRef Search ADS   Petiau G., Dupis A., 1980. Noise temperature coefficient and long-time stability of electrodes for telluric observations, Geophys. Prospect.  28( 05), 792– 804. https://doi.org/10.1111/j.1365-2478.1980.tb01261.x Google Scholar CrossRef Search ADS   Rangel G., 2014. Temperature model and tracer test analysis for the Ribeira Grande geothermal sield, (São Miguel island, Azores), Tech. Rep. UNU-GTP-2014-29, United Nations University, Geothermal Training Programme, (http://www.unugtp.is, last accessed 15 September 2017). Revil A., Cathles III L.M., Losh S., 1998. Electrical conductivity in shaly sands with geophysical applications, J. geophys. Res.  103( B10), 23 925–23 936. https://doi.org/10.1029/98JB02125 Revil A., Hermitte D., Spangenberg E., Cochemé J.J., 2002. Electrical properties of zeolitized volcaniclastic materials, J. geophys. Res.  107( B8), 2168, doi:10.1029/2001JB000599. https://doi.org/10.1029/2001JB000599 Google Scholar CrossRef Search ADS   Revil A., Johnson T.C., Finizola A., 2010. Three-dimensional resistivity tomography of Vulcan's forge, Vulcano Island, southern Italy, Geophys. Res. Lett.  37( 15), L15308, doi:10.1029/2010GL043983. Google Scholar CrossRef Search ADS   Revil A. et al.  , 2011. Hydrogeology of Stromboli volcano, Aeolian Islands (Italy) from the interpretation of resistivity tomograms, self-potential, soil temperature and soil CO2 concentration measurements, Geophys. J. Int.  186( 3), 1078– 1094. https://doi.org/10.1111/j.1365-246X.2011.05112.x Google Scholar CrossRef Search ADS   Sherburn S., Bannister S., Bibby H., 2003. Seismic velocity structure of the central Taupo Volcanic Zone, New Zealand, from local earthquake tomography, J. Volcanol. Geotherm. Res.  122( 1), 69– 88. https://doi.org/10.1016/S0377-0273(02)00470-5 Google Scholar CrossRef Search ADS   Sigmundsson F., Tryggvason E., Alves M.M., Alves J.L., Pálsson K., Ólafsson H., 1995. Slow inflation of the Furnas volcano, São Miguel, Azores, suggested from initial leveling and Global Positioning System measurements, Geophys. Res. Lett.  22( 13), 1681– 1684. https://doi.org/10.1029/95GL01604 Google Scholar CrossRef Search ADS   Silva C., Viveiros F., Ferreira T., Gaspar J.L., Allard P., 2015. Diffuse soil emanations of radon and hazard implications at Furnas Volcano, (São Miguel Island (Azores), in Volcanic Geology of São Miguel Island (Azores Archipelago)  chap. 15, pp. 5– 13, eds Gaspar J.L., Guest J.E., Duncan A.M., Barriga F., Chester D.K., Geological Society of London. Google Scholar CrossRef Search ADS   Simpson F., Bahr K., 2005. Practical Magnetotellurics  Cambridge University Press, Cambridge, UK. Google Scholar CrossRef Search ADS   Tikhonov A., 1950. The determination of the electrical properties of deep layers of the earth's crust, Dokl. Acad. Nauk. SSR  73 295– 297. Ussher G., Harvey C., Johnstone R., Anderson E., 2000. Understanding the resistivities observed in geothermal systems, in Proceedings World Geothermal Congress  Kyushu-Tohoku, Japan, 2000 May 28–June 10. Vanorio T., Virieux J., Capuano P., Russo G., 2005. Three-dimensional seismic tomography from P wave and S wave microearthquake travel times and rock physics characterization of the Campi Flegrei Caldera, J. geophys. Res.  110 B3, B03201, doi:10.1029/2004JB003102. https://doi.org/10.1029/2004JB003102 Google Scholar CrossRef Search ADS   Viveiros F., Cardellini C., Ferreira T., Caliro S., Chiodini G., Silva C., 2010. Soil CO2 emissions at Furnas volcano, São Miguel Island, Azores archipelago: Volcano monitoring perspectives, geomorphologic studies, and land use planning application, J. geophys. Res.  115( B12), B12208, doi:10.1029/2010JB007555. https://doi.org/10.1029/2010JB007555 Google Scholar CrossRef Search ADS   Vozoff K., 1972. The magnetotelluric method in the exploration of sedimentary basins, Geophysics  37 98– 141. https://doi.org/10.1190/1.1440255 Google Scholar CrossRef Search ADS   Wallenstein B., Duncan A., Chester D.K., Rui Marques D., 2007. Fogo volcano (São Miguel, Azores): a hazardous edifice, Géomorphologie : Relief, Processus, Environnement  13( 3), 259– 270. https://doi.org/10.4000/geomorphologie.2853 Google Scholar CrossRef Search ADS   Waxman M.H., Smits L.J.M., 1968. Electrical conductivities in oil-bearing shaly sands, Soc. Pet. Eng. J.  8 107– 122. https://doi.org/10.2118/1863-A Google Scholar CrossRef Search ADS   Weaver J.T., Agarwal A.K., Lilley F.E.M., 2000. Characterization of the magnetotelluric tensor in terms of its invariants, Geophys. J. Int.  141 321– 321. https://doi.org/10.1046/j.1365-246x.2000.00089.x Google Scholar CrossRef Search ADS   Weisenberger T.B., Ingimarsson H., Eyjólfsdóttir E.I., Lévy L., Hersir G.P., Ólafur G.Flóvenz, 2016. Validation of the influence of cation-exchange capacity on resistivity logs, in Proceedings European Geothermal Congress, Strasbourg, France, 19-24 Sept 2016 , (https://www.geothermal-energy.org/publications_and_services/conference_paper_database, last accessed 15 September 2017). Woitischek J., Dietzel M., Inguaggiato C., Böttcher M.E., Leis A., Cruz J.V., Gehre M., 2017. Characterisation and origin of hydrothermal waters at São Miguel (Azores) inferred by chemical and isotopic composition, J. Volcanol. Geotherm. Res.  346 104– 117. https://doi.org/10.1016/j.jvolgeores.2017.03.020 Google Scholar CrossRef Search ADS   Wright P.M., Ward S.H., Ross H.P., West R.C., 1985. State of the art geophysical exploration for geothermal resources, Geophysics  50( 12), 2666– 2696. https://doi.org/10.1190/1.1441889 Google Scholar CrossRef Search ADS   Yardley B.W.D., Bodnar R.J., 2014. Fluids in the continental crust, Geochem. Perspect.  3( 1), 1– 123. https://doi.org/10.7185/geochempersp.3.1 Google Scholar CrossRef Search ADS   Zandomeneghi D., Almendros J., Ibáñez J.M., Saccorotti G., 2008. Seismic tomography of Central São Miguel, Azores, Phys. Earth planet. Inter.  167( 1–2), 8– 18. https://doi.org/10.1016/j.pepi.2008.02.005 Google Scholar CrossRef Search ADS   SUPPORTING INFORMATION Supplementary data are available at GJI online. Figure S1. Apparent resistivity and phase responses plotted for all four components of the impedance tensor. While the solid circles represent the observed data, the solid lines represent the 3-D model responses. Figure S2. Apparent resistivity and phase responses plotted for all four components of the impedance tensor. While the solid circles represent the observed data, the solid lines represent the 3-D model responses. Figure S3. Apparent resistivity and phase responses plotted for all four components of the impedance tensor. While the solid circles represent the observed data, the solid lines represent the 3-D model responses. Figure S4. Apparent resistivity and phase responses plotted for all four components of the impedance tensor. While the solid circles represent the observed data, the solid lines represent the 3-D model responses. Figure S5. Maps of data fit (nRMS) for the high-frequency bands 10,000 1,000 Hz and 1,000 100 Hz. Figure S6. Maps of data fit (nRMS) for the low-frequency bands 100 10 Hz and 10 1 Hz. Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper. APPENDIX: CALCULATION OF SMECTITE MASS FRACTION Here, we summarize the parameter settings used for the calculation of smectite content following the studies by Waxman & Smits (1968), Revil et al. (1998, 2002) and Flóvenz et al. (2005). In Fig. 10, we took a broad interval of fluid conductivity, temperature, and porosity based on the borehole data and modelling results from São Miguel Island (Carvalho et al. 2006), the borehole data obtained by Flóvenz et al. (2005) from a geothermal system in Iceland that we consider as a comparable system to Furnas Volcano, and the laboratory results of Revil et al (2002) on altered samples (see Table A1). The purpose of this study is to demonstrate that the fluid conductivity alone cannot explain the observed value of bulk resistivity of the C1 body; a certain amount of clays characterized by high surface conductivity has to be assumed to explain our inversion results. We suggest the presence of smectite as it has the highest surface conductivity: the mass fraction of about 20 per cent of smectite is enough to explain the observed bulk conductivity of the conductive zone C1 and would still be coherent with high permeability suggested in the fault zone (FF2) from geochemical observations. An additional evidence is given by borehole data in Ribeira Grande (Franco 2015) where the smectite-layer was observed at a depth that is comparable to the depth of the C1 conductive zone. The smectite can be used as indicators of temperature; they are found at a certain temperature interval only and are gradually replaced by less conductive chlorites at a temperature around 200–230 °C (e.g. Franco 2015). The transition from the low to high resistivity below the C1 body (the R2 resistive body or alternatively, the gradual increase in resistivity in the C2 body) could signify this temperature transition. Table A1. Physical properties used for calculation of smectite mass fraction (Fig. 10). Density of rock (borehole data from Flóvenz et al. 2005)  ρ  (kg m−3)  Eq. (A3)  Fluid conductivity  σf  (S m−1)  0.02–3  Hittorf number for NaCl brines (Revil et al. 2002)  t(+)    0.38  Cation exchange capability (Revil et al. 2002)  CEC  (C kg−1)  87 000  Porosity ((borehole data from Flóvenz et al. 2005)  ϕ    0.02–0.25  Cementation exponent (Flóvenz et al. 2005)  m    2.75  Tortuosity (Flóvenz et al. 2005)  a    0.7  Temperature range    (°C)  80–200  Temperature coefficient (eq. A4)  υf    0.02  Temperature coefficient (eq. A4)  υs    0.04  Apparent mobility of counterions at 25 °C (Revil et al. 2002)  β(s)  (m2 s−1 V−1)  5.2 × 10−9  Density of rock (borehole data from Flóvenz et al. 2005)  ρ  (kg m−3)  Eq. (A3)  Fluid conductivity  σf  (S m−1)  0.02–3  Hittorf number for NaCl brines (Revil et al. 2002)  t(+)    0.38  Cation exchange capability (Revil et al. 2002)  CEC  (C kg−1)  87 000  Porosity ((borehole data from Flóvenz et al. 2005)  ϕ    0.02–0.25  Cementation exponent (Flóvenz et al. 2005)  m    2.75  Tortuosity (Flóvenz et al. 2005)  a    0.7  Temperature range    (°C)  80–200  Temperature coefficient (eq. A4)  υf    0.02  Temperature coefficient (eq. A4)  υs    0.04  Apparent mobility of counterions at 25 °C (Revil et al. 2002)  β(s)  (m2 s−1 V−1)  5.2 × 10−9  View Large A number of parameters characterizing the pore fluid or the rock are used in the equations listed below: the fluid conductivity σf, the surface conductivity σs; their ratio is a so-called Dukhin number ζ = σs/σf. The t(+) is the Hittorf number, the fraction of electrical current transported by cations in the electrolyte (see Table A1 for values). The F = aϕ−m is a modified intrinsic formation factor (Archie 1942), depends on the porosity of rock ϕ, tortuosity a, and the cementation exponent m. At low salinities (ζ ≥ t(+)), the bulk conductivity σ is calculated with the Waxman–Smits model (Waxman & Smits 1968; Revil et al. 2002) as a combination of the surface and fluid conductivity contributions   \begin{equation} \sigma = \sigma _s + \frac{\sigma _f}{F}\left( 1 - t_{(+)}\right). \end{equation} (A1) The surface conductivity σs is a function of the cation exchange capacity CEC,   \begin{equation} \sigma _s = \frac{2}{3} \rho \beta _s {\rm CEC} \end{equation} (A2)where the ρ is the density, while βs is the surface mobility of the counterions (Table A1). Smectite bulk density can be calculated following the regression determined in the work of Flóvenz et al. (2005) as   \begin{equation} {\rho _{\rm{smectite}}} = 2.95 - 0.032\phi. \end{equation} (A3) The dependence of σf and σs on temperature can be approximated linearly by   \begin{equation} \sigma _{f,s}(T) = \sigma _{f,s}(T_0)\left[ 1 + \upsilon _{f,s}\left( T - T_0\right) \right], \end{equation} (A4)where T is the temperature, T0 is its reference value (20 °C), and υf, s are the regression coefficients for fluid and surface conductivities, respectively. Following Revil et al. (1998), at higher salinities the bulk electrical conductivity can be expressed as the sum of conductivity contributions by cations (σ(+)) and anions (σ(−)) as   \begin{equation} \sigma = \sigma _{(+)} + \sigma _{(-)}. \end{equation} (A5)These contributions are given by   \begin{eqnarray} \sigma _{(+)} &=& \frac{\sigma _f}{F} \left[ F \zeta + \frac{1}{2} ( t_{(+)} - \zeta )\right. \nonumber\\ &&\left.\times \left( 1 - \frac{\zeta }{t_{(+)} } + \sqrt{\left( 1 - \frac{\zeta }{t_{(+)}} \right) ^2 + \frac{4F}{t_{(+)}} \zeta } \right) \right], \end{eqnarray} (A6)and   \begin{equation} \sigma _{(-)} = \frac{\sigma _{f}}{F} \left( 1 - t_{(+)}\right). \end{equation} (A7) © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

3-D interpretation of short-period magnetotelluric data at Furnas Volcano, Azores Islands

Loading next page...
 
/lp/ou_press/3-d-interpretation-of-short-period-magnetotelluric-data-at-furnas-w3EmU35Xx3
Publisher
The Royal Astronomical Society
Copyright
© The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx512
Publisher site
See Article on Publisher Site

Abstract

Summary Accurate geophysical imaging of shallow subsurface features provides crucial constraints on understanding the dynamics of volcanic systems. At Furnas Volcano (Azores), intense circulation of volcanic fluids at depth leading to high CO2 outgassing and flank destabilization poses considerable threat to the local population. Presented is a novel 3-D electrical resistivity model developed from 39 magnetotelluric soundings that images the hydrothermal system of the Furnas Volcano to a depth of 1 km. The resistivity model images two conductive zones, one at 100 m and another at 500 m depth, separated by a resistive layer. The shallow conductor has conductivity less than 1 S m−1, which can be explained by clay mineral surface conduction with a mass fraction of at least 20 per cent smectite. The deeper conductor extends across the majority of the survey area. This deeper conductor is located at depths where smectite is generally replaced by chlorite and we interpret it as aqueous fluids near the boiling point and infer temperatures of at least 240 °C. The less conductive layer found between these conductors is probably steam-dominated, and coincides within the mixed-clay zone found in many volcanic hydrothermal systems. Electrical properties, Hydrothermal systems, Magnetotellurics, Numerical modeling 1 INTRODUCTION The Azores archipelago (Portugal) is formed by nine volcanic islands located in the North Atlantic Ocean where the American, Eurasian and African plates meet at a triple junction (Buforn & Udías 2010). São Miguel Island is the largest (≈65 × 15 km2) of the archipelago and is located over the Terceira Rift, which extends from the island of Santa Maria towards the Mid-Atlantic Ridge following approximately a NW–SE trend (Machado 1959). São Miguel Island comprises three trachytic polygenetic volcanoes: Sete Cidades, Fogo (Água de Pau) and Furnas (Fig. 1). Figure 1. View largeDownload slide The Azores islands are located at the triple junction between the North American, Eurasian and African plates. São Miguel Island has three active strato-volcanoes, Sete Cidades, Fogo (Água de Pau) and Furnas. The red frame outlines the survey area of this work around the inner caldera at Furnas. Modified from Guest et al. (1999) and Wallenstein et al. (2007). Figure 1. View largeDownload slide The Azores islands are located at the triple junction between the North American, Eurasian and African plates. São Miguel Island has three active strato-volcanoes, Sete Cidades, Fogo (Água de Pau) and Furnas. The red frame outlines the survey area of this work around the inner caldera at Furnas. Modified from Guest et al. (1999) and Wallenstein et al. (2007). Furnas, a volcanic complex, is the eastern-most of the three trachytic active central volcanoes of São Miguel Island. It does not have a well-developed edifice and is a superposition of two calderas. The external and largest caldera was formed about 30 ka BP, followed by a period of caldera filling eruptions. The inner caldera is smaller and was formed around 10–12 ka BP, which is responsible for the steep topography of more than 500 m in our study area. It contains several younger eruptive centres, a shallow caldera lake and several zones of intensive fumarolic activity within the inhabited caldera (Guest et al. 1999, 2015). Phlegrean fields volcano (Naples, Italy) may be comparable to Furnas volcano in several aspects, such as the existence of a caldera and the composition of the hydrothermal fumarolic emissions. The last magmatic eruption in this volcanic system took place in 1630 AD (Cole et al. 1995) causing about 200 fatalities. Present-day activity comprises secondary manifestations of volcanism that are characterized by boiling temperature fumaroles, steaming ground, thermal springs, cold CO2 rich springs and areas of intense diffuse degassing especially in the northern part of the inner caldera, north of Furnas Lake, where detailed maps of surface CO2 and 222Rn emissions were recently made available (Viveiros et al. 2010; Silva et al. 2015). CO2 emissions can pose serious threats to populations that reside in volcanic regions. This is for instance the case at Furnas, where 60-90 per cent of dwellings are built over hydrothermal CO2 emission zones (Viveiros et al. 2010). In addition, hydrothermal explosions were reported in one of the fumaroles located close to Furnas village (Pedone et al. 2015, and references therein). Along with this risk, intense circulation of volcanic fluids at depth, coupled with alteration can lead to associated mechanical weakening of the volcanic edifice and may lead to caldera wall instabilities. This process may be triggered and/or coupled with elevated rainfall. The evaluation and forecasting of volcanic risks require the understanding of the geometry of the hydrothermal circulation system and quantifying the relationship between the surface degassing features and deeper processes driving the heat and fluid flows. The complex tectonic structure of Furnas is affected by fault systems trending WNW–ESE, NE–SW, N–S and NNW–SSE, expressed by volcanic alignments, linear valleys and caldera outlines (Carmo et al. 2015). The caldera walls have orientations that coincide with the main fault trends (WNW–ESE and NE–SW), suggesting that caldera collapses were controlled by structural weaknesses as proposed by Guest et al. (1999) and Carmo et al. (2015). The tectonic regime of the area is thus complex; because of the low amount of local seismicity, the recent state of tectonic stress is not well constrained. Even less is known on the existence, location, and geometry of the magmatic source, though the consensus is that there is a magmatic reservoir region, constrained to shallow crustal depths of 3 to 6 km beneath the Furnas Caldera by tectonics (Machado 1959), GPS (Sigmundsson et al. 1995), magnetic data (Blanco et al. 1997), gravity data (Camacho et al. 1997; Montesinos et al. 1999), seismic tomography (Zandomeneghi et al. 2008) and geochemistry (Jeffery 2016). The exact location and structure of the magmatic reservoir is somewhat unknown, with indications of a position south of the caldera centre (Zandomeneghi et al. 2008, and references therein). High electrical conductivity values are typically associated with volcanic-hydrothermal systems and the elevated conductivity structure can provide constraints on these volcanic and hydrothermal processes (Moore et al. 2008; Revil et al. 2010). These settings are a good target for electrical geophysical methods that are sensitive to conductivity. The resistivity (the reciprocal of conductivity) of a geological unit is a function of the amount and interconnectivity of a conductive fluid contained within the matrix of the lithology and also the surface conductivity. In general, bulk resistivity decreases with increases in pore water connectivity, saturation, salinity, temperature, and clay content due to rock alteration (e.g. Revil et al. 1998; Ussher et al. 2000; Revil et al. 2002; Muñoz 2014). Studies to date have found that conductive areas correspond to the pathway and phase states of hydrothermal fluids, and to the fluid bearing structure in volcanic and geothermal areas (e.g. Ogawa et al. 1998; Aizawa et al. 2005, 2009; Revil et al. 2011). Thus, resistivity anomalies can provide important constraints leading to an improved understanding on structure and dynamics of hydrothermal systems. The first magnetotelluric (MT) measurements in the Furnas Caldera date from 1983, when the United States Geological Survey (USGS) carried out geophysical surveys for assessing geothermal potential on the three volcanoes on São Miguel Island (Hoover et al. 1983, 1984). Though using only scalar Audio-MagnetoTelluric (AMT) soundings, they were able to identify zones of high geothermal potential, in particular, the geothermal reservoir at the NW flank of the Fogo (Água de Pau) volcano, which now hosts two geothermal power plants (Franco 2015) providing more than a third of the São Miguel's energy demand. Aiming at an improved imaging of subsurface electrical conductivity, and a better understanding of the shallow, fumarole-related features within the Furnas hydrothermal system, 39 AMT including 15 Broad-Band MagnetoTelluric (BBMT) systems were deployed in 2015 and 2016. In this study, we focus on the AMT portion of the data in order to identify the main features of the Furnas hydrothermal system. The use of the remaining long-period data requires accurate knowledge of the bathymetry around the island and will be presented in a follow-up paper. 2 PREVIOUS STUDIES AT FURNAS VOLCANO Furnas volcano is situated on the eastern part of the island and its earliest lithology is dated from about 100 000 yr BP (Moore 1990). It comprises a multicyclic forming caldera which last erupted in 1630 AD. Furnas exhibits a 5 × 8 km summit depression formed by juxtaposed calderas controlled by NW–SE and NE–SW faults. Conjugate faulting, parallel to the Terceira Rift regional fault system are observed and exhibit a NW–SE trend (Guest et al. 1999). The present caldera system is flanked by steep walls as high as 500 m, that formed following multi-collapse episodes and subsequent caldera-filling eruptions. The younger caldera contains several craters, and the western part is covered by a lake that occupies an area of 1.9 km2 (Cruz & França 2006). Recent volcanic activity at Furnas includes two historic events since the island became inhabited in the 15th century. These two intracaldera volcanic eruptions have occurred at Furnas volcano in 1439–1443 (which led to the creation of Gaspar Dome also referred to as Pico do Gaspar) and in 1630 (referred to as the ‘1630 Dome’) that led to the formation of two tuff rings with central domes that occupy the central part of the caldera floor as shown in Fig. 2. A detailed discussion on the composition, stratigraphy and eruptive history of Furnas can be found in Guest et al. (1999, 2015, and references therein). At present, volcanic activity of Furnas volcano is characterized by the secondary manifestations of volcanism mostly located in the northern region of the caldera, within Furnas village and north of the Furnas Lake, and are characterized by boiling temperature fumaroles, steaming ground, thermal springs, cold CO2-rich springs and areas of intense diffuse degassing, namely CO2 and radon (Viveiros et al. 2010; Caliro et al. 2015; Silva et al. 2015). Figure 2. View largeDownload slide Locations of MT stations around the Furnas volcano. The red inverted triangles represent AMT only sites, while the blue ones represent combined AMT and BBMT sites. The structural features are also mapped: solid line (fault), dashed line (inferred fault) and dotted line (volcanic alignment), courtesy of Carmo et al. (2015). The main fumarole area at Furnas Lake is marked by the yellow circle. The pink circle marks the location of site FUR007A (Fig. 3). Elevation data provided by the IVAR GIS Group at the University of the Azores. Figure 2. View largeDownload slide Locations of MT stations around the Furnas volcano. The red inverted triangles represent AMT only sites, while the blue ones represent combined AMT and BBMT sites. The structural features are also mapped: solid line (fault), dashed line (inferred fault) and dotted line (volcanic alignment), courtesy of Carmo et al. (2015). The main fumarole area at Furnas Lake is marked by the yellow circle. The pink circle marks the location of site FUR007A (Fig. 3). Elevation data provided by the IVAR GIS Group at the University of the Azores. The subsurface structure of the Furnas volcanic-hydrothermal system is not well constrained. In 1985, two aeromagnetic surveys were made covering Faial and São Miguel Islands. The results, first published by Miranda et al. (1991), yielded high amplitude positive magnetic anomalies, coherent with the existence of highly magnetized and young eruptive rocks covering the flanks of most of the volcanic systems on the island, however two calderas, Sete Cidades and Furnas match relative magnetic lows, whereas Fogo (Água de Pau) corresponds to a large negative magnetic anomaly. Miranda et al. (2015) inferred the decrease in magnetization at Furnas to be the result of hydrothermal alteration, similar to the conclusion derived from Blanco et al. (1997) for the decrease in the magnetic field anomaly at Furnas, as sensed by land based magnetic surveys. The intense degassing and fumarolic features at Furnas (Viveiros et al. 2010; Silva et al. 2015) suggest the presence of such a mature hydrothermal system beneath the central volcano and supports an equivalent interpretation. As mentioned, Blanco et al. (1997) disclosed the existence of magnetic anomalies within the Furnas area. Forjaz (1996) documents that this alteration would affect the shallowest materials more intensively within the caldera that are derived from different structures in and around the volcano. The inner part of the caldera, around Furnas Lake and the historic volcanic domes which are most relevant to our study, coincides with a broad negative magnetic anomaly. The authors suggest that the cause may be related to the oxidation of the magnetic minerals as a result of hydrothermal alteration from the surface to a depth of 500 m, where the hydrothermal system seems most active, however its effect probably exists at much greater depths. The magnetic anomaly map derived from this work is only sensitive to structures to a depth of approximately 3 km, thus the proposed magma chamber, presumed to be at a depth of at least 5 km, could not be imaged. Geodetic studies, between 1991 and 1994, across central São Miguel revealed a slight inflation in the area of Furnas, with the centre of this suggested inflation in the northern part of the caldera. Sigmundsson et al. (1995) infers this phenomenon being caused by growing pressures of fluids within a deep hydrothermal system beneath the caldera which is sealed by an impermeable zone. It is worth noting that the authors suggest that the centre of deformation is located within the zone of the negative anomaly detected by Blanco et al. (1997). Camacho et al. (1997) and Montesinos et al. (1999) analysed and modelled gravity data collected around Furnas, yielding gravity anomaly maps that isolated low-density zones that may be the result of several distinct caldera collapses and infilling. One such anomaly, a principal gravity minimum is located below Pico do Gaspar, just east of Furnas Lake, and is quite distinct at a depth of 1000 m from the work of Montesinos et al. (1999). This structure is probably based on the formation of the inner caldera that cuts the old caldera followed by another phase of infilling. Similar low-gravity anomalies related to low-density products of explosive volcanic activity and caldera collapse processes were reported by Sherburn et al. (2003) for Taupo zone calderas and Vanorio et al. (2005) for Campi Flegrei. A seismic tomography study covering the centre of São Miguel between the Fogo (Água de Pau) and Furnas volcanoes was conducted by Zandomeneghi et al. (2008), identifying a low velocity region beneath Furnas. It is characterized by low P-wave velocities and low Vp/Vs ratios, which decrease with depth until they reach a minimum of 2 km below sea level. The authors’ results coincide spatially with the gravity inversion low-density anomaly. The authors infer the low Vp/Vs ratios changes being related to changes in fluid conditions within an extensive steam-dominated geothermal system. Geothermal exploration and development commenced on São Miguel as early as 1975, however, this exploration was confined to the area on the northern slope of Fogo (Água de Pau). In 1982 and 1983, the USGS conducted reconnaissance AMT surveys of the three quiescent central volcanoes on the island, to evaluate the potential for geothermal systems (Hoover et al. 1983, 1984). While their survey yielded reasonable results, it was strongly limited by the narrow frequency range and use of a scalar system, that is, only measuring the electric (E) and corresponding orthogonal magnetic (H) field in an assumed fixed coordinate system (NW–SE and NE–SW). Thus only two of the four elements of the impedance tensor (eq. 1, Section 3) were determined. Nevertheless, their results at Furnas indicated an east-west trending, low resistivity zone through the centre of the caldera. This zone correlates with the known fumaroles and hot springs at Furnas Lake and at Furnas village. The observed apparent resistivities at this zone were estimated to be 10 Ωm from their apparent resistivity maps at 7.5 Hz (this was the lowest frequency they recorded at which signal strengths were adequate to give good data quality) in the near surface and drop to approximately 1.5 Ωm at a depth of 50 m from 1-D estimates. Similar low apparent resistivities were observed to the south and east of the caldera near the coast, however the authors suggest that sea water intrusion on these data cannot be ruled out because of the proximity of the Atlantic Ocean, however it should also be noted that the studies of Viveiros et al. (2010) and Silva et al. (2015) also identify fumarolic fields and soil diffuse degassing anomalous areas (CO2 and 222Rn) at the coast south of the Furnas Caldera near Ribeira Quente, concluding that localized geothermal systems are present outside the caldera. 3 MT METHOD The MT method is a passive geophysical exploration technique that yields information on the 3-D spatial distribution of electrical conductivity of the Earth's subsurface by measuring and relating the natural time-varying horizontal electric (E) and magnetic (H) fields at its surface (Tikhonov 1950; Cagniard 1953). The relationship between horizontal and mutually perpendicular field components provides amplitude responses (apparent resistivity) and phase responses at various frequencies, known as MT response functions, for each site recorded. Amplitudes of the fields decrease exponentially with increasing depth in a uniform conductor, known as the skin-depth phenomenon. The depth of penetration of the MT signal into the earth depends on both the conductivity of the rocks and the frequency of the variation of the MT signal; with low-frequency (or long-period) variations sensing deeper. The conductivity structure of the subsurface can be determined from the amplitude and phase relationships between the surface electric and magnetic field vectors as a function of period and sounding location. These relationships are expressed mathematically by the impedance tensor Z, where Z (Vozoff 1972), is a complex 2 × 2 matrix as   \begin{equation} (E_x\ E_y) = \left( {\begin{array}{*{20}{c}}{Z_{xx}}&{{Z_{yx}}} \\ {{Z_{xy}}}&{{Z_{yy}}} \end{array}} \right)\left( {\begin{array}{*{20}{c}}{H_x} \\ {{H_y}} \end{array}} \right). \end{equation} (1) From the elements of the impedance tensor, the more intuitive apparent resistivities and phases can be calculated as   \begin{equation} {\rho _{a,ij}} = \frac{{{\mu _0}}}{\omega }{\left| {{Z_{ij}}(\omega )} \right|^2},{{\rm }}{\phi _{ij}} = \arctan \frac{{\Im ({Z_{ij}})}}{{\Re ({Z_{ij}})}} \end{equation} (2)where μ0 is the free-space magnetic permeability. The angular frequency ω is defined 2πf, while ℜ and ℑ represent the real and imaginary parts, respectively. A further quantity calculated from the MT data is the local magnetic transfer function (also complex) between the horizontal and vertical magnetic field components, T, commonly called induction vector. For this study the sign of this vector was chosen so that the induction vector points in the direction of increasing conductivity (Parkinson 1962):   \begin{equation} \left({\begin{array}{*{20}{c}}{H_z} \end{array}} \right) = \left( {\begin{array}{*{20}{c}}{T_{x}}&{{T_{y}}} \end{array}} \right)\left( {\begin{array}{*{20}{c}}{H_x} \\ {{H_y}} \end{array}} \right). \end{equation} (3) All of these quantities provide information on the subsurface distribution of resistivity. In particular, the structure of the impedance tensor and its rotational invariants can be used to determine whether data may be interpreted in a 1-D or 2-D model, or require 3-D modelling. Further details can be found in recent comprehensive monographs (such as Simpson & Bahr 2005; Chave & Jones 2012). 4 FIELD DATA AND PROCESSING Following a pilot survey of 13 AMT soundings and Electrical Resistivity Tomography (ERT) data collected along two profiles in the eastern part of Furnas Caldera in 2015, a second campaign was completed in June 2016, yielding a total of 39 soundings including 15 BBMT soundings, using the Phoenix Geophysics MTU-5A data-logger along with MTC-30 and MTC-50 induction coils for AMT and BBMT magnetic field data, respectively. Due to the steep topography within the caldera, obtaining a site distribution in a regular grid proved impossible, however good coverage was obtained within and outside the caldera (Fig. 2). The survey area was free of major noise sources. At each site, the two horizontal, perpendicular magnetic field components Hx and Hy were recorded as well as the vertical magnetic field component Hz. The two horizontal, perpendicular electric field components Ex and Ey were measured using non-polarizing Phoenix Geophysics Pb–PbCl (lead–lead chloride) electrodes (Petiau & Dupis 1980), laid out in a cross with a dipole length of typically 80 m. The time-series were recorded overnight and have a duration of approximately 12 hr. Time-series processing techniques applied to the data have yielded apparent resistivity and phase response curves for each site. The acquired MT data were processed from the raw time-series via commercial remote-reference software from Phoenix Geophysics (based on Jones & Jödicke 1984). Remote reference processing was computed with concurrent sites and also with the remote site, however processed data from the remote site did not always improve the final response significantly for the MT sites. The best response (i.e. the smoothest transfer function with the best responses in the low signal frequency band between 1 and 5 kHz) from the far and local remote site was selected as the final response of the individual site. Both magnetic field and electric field references were processed, however the magnetic field data were used as it gave superior responses. Following processing, the physical consistency of the resistivity and phase estimates at each station were checked using the D+ algorithm (Parker 1980; Parker & Whaler 1981; Beamish & Travassos 1992; Parker & Booker 1996), using the WinGlink® software package. Outliers and inconsistent data points were removed from the data set, prior to data analysis and inversion. Site locations with both AMT and BBMT data had their responses merged to yield one response for the site. An example of a high-frequency MT response can be seen in Fig. 3. With the full five-component measurements now available, it was possible to analyse the dimensionality and directionality of the observations with different methods used within the MT community. Figure 3. View largeDownload slide Response curves of an example AMT site. Apparent resistivity (ρa, top) curves, phase curves (ϕ, middle), and induction arrows (bottom) in Parkinson convention are shown. The site location is marked with a pink circle on Fig. 2. Figure 3. View largeDownload slide Response curves of an example AMT site. Apparent resistivity (ρa, top) curves, phase curves (ϕ, middle), and induction arrows (bottom) in Parkinson convention are shown. The site location is marked with a pink circle on Fig. 2. 5 DATA ANALYSIS Prior to analysing and modelling of the acquired MT data, the presence of the Atlantic Ocean around the survey area was investigated to quantify the sharp electrical contrast between the ocean and land. The skin-depth of a typical BBMT sounding can easily reach great depths (several tens of km) depending on the actual conductivity structure, thus the ocean has a significant impact on the observed responses when the separation distance between the coast and the survey area is smaller than the depth of interest. This makes it difficult to determine reliable geoelectrical structures at deeper regions of the subsurface. In order to quantify the impact of the Atlantic Ocean, forward modelling was used to create synthetic responses from a pre-determined 3-D resistivity model. A simple 3-D model of São Miguel Island with a seawater resistivity value of 0.3 Ωm extending to a depth of 6000 m and a conservative uniform land resistivity of 100 Ωm was designed. The results from this forward modelling, indicate that the ocean does not impact on data with periods less than 1 s, therefore for the purpose of this study, only data up to a period of 1 s are used in the analysis and modelling. Analysis of the geoelectric dimensionality of MT data has acquired relevance due to the advancement of numerical codes that have made it possible, with confidence, to model and invert the data using 1-D, 2-D or 3-D approaches. Simultaneously, one can gain insights into the direction of the geoelectrical strike at various depths, which in turn can be used to qualitatively correlate with different geological processes and subsurface geological structures. A dimensionality study focussing on the rotational invariants of the MT tensor was presented by Weaver et al. (2000), however, in real and noisy data, the invariants are rarely zero, yielding geoelectric dimensionality of noisy responses to be interpreted as 3-D. Weaver et al. (2000) overcame this difficulty by implementing a threshold value, τ, beneath which the invariants are set to zero. The work of Martí et al. (2005) and Martí et al. (2010), that led to the development of WALDIM (Martí et al. 2009), a FORTRAN application to perform dimensionality analysis of MT data, found that a value of τ between 0.1 and 0.2 determines the dimensionality very well, and that the geoelectric strike direction can be estimated, when 2-D structure is inferred through a Monte Carlo approach with the addition of Gaussian noise. Along with classifying the dimensionality into bands of periods for each site, the WALDIM application has been extended to identify hints of anisotropy through the assessment of the invariants and strike directions from a pattern of sites and frequencies. According to Martí et al. (2010), these indicators of anisotropy occur in the inferred 2-D cases when there is an inconsistency in strike directions and when relationships between the tensors do not imply isotropic media. In this study, a τ of 0.15 was used and the dimensionality of each site was estimated within each period decade band. The results, along with the estimates of geoelectrical strike direction whenever quasi-2-D cases are found can be seen in Fig. 4. Figure 4. View largeDownload slide Results of dimensionality analysis and geoelectric strike estimations (purple coloured arrows) for four bands using the WALDIM code (Martí et al. 2009), superimposed on the main structural lineaments in the Furnas region as derived from Carmo et al. (2015). Note the dominant change in dimensionality of the geoelectrical subsurface from 1-D (purple squares) at the highest frequency range (10 000 to 1000 Hz) to a more 3-D (green squares) environment at the deepest regions (10 to 1 Hz). Hints of anisotropy are also estimated and can be seen to be coincident with steep topography and adjacent to main structural features. Legend (following the definitions of Martí et al. (2009): 1-D, one dimensional; 2-D, two dimensional; 3-D/2-Dtw, two -dimensional data affected by distortion (only twist); 3-D/2-D, general case of galvanic distortion over a 2-D structure; 3-D, affected or not by galvanic distortion; 3-D/2-Ddrt, 3-D or 2-D with diagonal regional tensor; 3-D/2-D or 3-D/1-D, galvanic distortion over a 1-D or 2-D structure with a non-recoverable strike direction; Anis1, homogeneous anisotropic medium; Anis 2, anisotropic body within a 2-D medium. Figure 4. View largeDownload slide Results of dimensionality analysis and geoelectric strike estimations (purple coloured arrows) for four bands using the WALDIM code (Martí et al. 2009), superimposed on the main structural lineaments in the Furnas region as derived from Carmo et al. (2015). Note the dominant change in dimensionality of the geoelectrical subsurface from 1-D (purple squares) at the highest frequency range (10 000 to 1000 Hz) to a more 3-D (green squares) environment at the deepest regions (10 to 1 Hz). Hints of anisotropy are also estimated and can be seen to be coincident with steep topography and adjacent to main structural features. Legend (following the definitions of Martí et al. (2009): 1-D, one dimensional; 2-D, two dimensional; 3-D/2-Dtw, two -dimensional data affected by distortion (only twist); 3-D/2-D, general case of galvanic distortion over a 2-D structure; 3-D, affected or not by galvanic distortion; 3-D/2-Ddrt, 3-D or 2-D with diagonal regional tensor; 3-D/2-D or 3-D/1-D, galvanic distortion over a 1-D or 2-D structure with a non-recoverable strike direction; Anis1, homogeneous anisotropic medium; Anis 2, anisotropic body within a 2-D medium. The results of the dimensionality analysis indicate that the Furnas Caldera region is quite complex, with dimensionality ranging from 1-D to 3-D for the four frequency bands. At the highest frequency band (10000 Hz – 1000 Hz), 1-D structure is dominant across the whole survey area, complementing the overlapping XY and YX modes in the MT responses. As the longer-period data probes deeper into the subsurface, the geoelectrical structure becomes more 3-D and this can be clearly seen in the 10 Hz – 1 Hz panel. Hints of anisotropy can be seen throughout the data set in all frequency bands, however the majority of these sites appear to be located in regions of steep topography (particularly north of Furnas Lake), and also in structurally controlled regions. The WALDIM code was also used to quantify the geoelectrical strike direction within the data set where the dimensionality showed 2-D structure and where hints of anisotropy were sensed. The general trend of this strike is NE–SW, with the average strike direction for each decade ranging from 42° to 55°. The overall complexity of the subsurface structure beneath Furnas can be further illustrated by plotting maps of the phase tensor ellipses from all sites as a function of period. The MT phase response becomes more sensitive to structure at greater depth and less sensitive to structure at shallow depth as the period increases. In contrast, the magnitude or apparent resistivity response will be distorted by the galvanic effect of near-surface conductivity heterogeneities at all periods. However, galvanic effects do not significantly affect the phase tensor (Caldwell et al. 2004) and phase tensor ellipse maps (Fig. 5) provide a method of visualizing spatial gradients in the conductivity structure directly. The principal axes of the phase tensor ellipse, ϕmax and ϕmin, indicate the horizontal directions of the maximum and minimum induction current, which reflects lateral variations in the resistivity structure. The variation of the orientation of the major axis of the ellipse can help identify changes in the dimensionality of the subsurface structure. In the case of 1-D structure, ϕmax = ϕmin and therefore the ellipse will be a circle. In the 2-D case, ϕmax ≠ϕmin and the phase tensor is an ellipse but also symmetric, however when a 3-D subsurface is probed, the ellipse is not symmetric. Four maps at frequencies that are approximately halfway within each decade band along with induction arrows are plotted in Fig. 5 in which the ellipses are scaled so that the long axis (representing ϕmax) is of uniform length, and are colour coded according to the value of ϕmin. At the highest frequency ranges (centred at 5263 Hz and 532 Hz), a significant amount of the ellipses in the survey area are quite circular, inferring the near surface conductivity structure to be essentially 1-D. This compliments the results from the WALDIM dimensionality study. The only exceptions to these are a number of sites that show a small degree of ellipticity which may be a response to the steep topography. At these regions, the magnitude of ϕmin is relatively low (colder colours) inferring a geoelectrical subsurface that is resistive. The magnitude of ϕmin shows that much higher phase values—an indicator of increasing conductivity (warmer colours) with depth—occur over the whole survey area at longer periods which can be seen in the frequency ranges (centred at 49 and 5.6 Hz). A remarkable conductive area is clearly shown just beneath the northeastern edge of Furnas Lake (Figs 5c and d). Induction arrows can be used to illustrate the relationship between the measured vertical and horizontal magnetic field components. Adopting the convention proposed by Parkinson (1962), the real part of such arrows point towards conductive bodies or concentrations of induced currents. Analysis of the induction arrows confirm that there is a conductive body beneath the region of Furnas Lake (Figs 5b and c), however, at the longest periods (Fig. 5d) the induction arrows provide a strong indication of the presence of a larger-scale conductivity structure located beyond the eastern and the northeastern part of Furnas Lake. Figure 5. View largeDownload slide Maps of phase tensor ellipses and real part of the induction arrows at four frequencies. The colour of the phase tensor ellipses (axes normalized by ϕmax) indicates ϕmin. The induction arrows are plotted using the Parkinson convention. Figure 5. View largeDownload slide Maps of phase tensor ellipses and real part of the induction arrows at four frequencies. The colour of the phase tensor ellipses (axes normalized by ϕmax) indicates ϕmin. The induction arrows are plotted using the Parkinson convention. Following this analysis, it is clear that the conductivity structure associated with the hydrothermal system beneath Furnas is 3-D and the observed MT data need to be inverted accordingly. Furthermore, the steep topography associated with the caldera can only be accounted for by true 3-D inverse modelling. 6 3-D MODELLING As discussed in Section 5, the potential influence of the Atlantic Ocean on the data set is absent for periods ≤1 s. The 3-D inversions thus were performed including only high-resolution topography and Furnas Lake bathymetry data employing the parallel version of the Modular system for Electromagnetic inversion code ModEM (Meqbel 2009; Egbert & Kelbert 2012; Kelbert et al. 2014). The 3-D mesh consisted of 119 × 119 × 149 (plus 10 air layers) in the north–south, east–west and vertical directions, respectively. The centre of the mesh comprised a uniform mesh of 89 × 89 × 149 cells, with a horizontal cell size of 50 m × 50 m. The lateral extents of the padding cells, 15 planes in each of the north, south, east and west directions, increased by a factor of 1.5. In the vertical direction, the cell thickness was of 7.5 m for 698.264 m < z ≤ 817.5 m; subsequent layer thicknesses were increased by a factor of 1.2. The input data were decimated to 17 periods per site with regular distribution for periods from 0.0001 to 1 s. In total, 39 sites were used in the inversion process. The 3-D inversions were run with full impedance tensor elements (Zxx, Zxy, Zyx and Zyy) and vertical magnetic transfer functions. Error floors were set to an absolute value of 3 per cent of (∣Zxy × Zyx∣)1/2. For the vertical transfer functions, a constant error of 0.015 was used. The starting model resistivity was set to 100 Ωm. The high-resolution bathymetry of the lake was included as a priori information and kept fixed during the inversion. A resistivity value of 63 Ωm was used for the lake (Andrade et al. 2016). Model regularization employed a smoothing parameter of 0.4 for both horizontal directions, and 0.3 for the vertical direction. A commonly used metric for the fit of the modelled to the observed data is the normalized root mean square nRMS, defined as   \begin{equation} {\rm nRMS} = \sqrt{\left( {{N_d} - 1} \right)_{}^{ - 1} \sum \nolimits _{j = 1}^{N_d} \left( \frac{d_{{\rm obs},j} - d_{{\rm calc},j}}{\sigma _j} \right)^2 } {\, }, \end{equation} (4)where dobs,j and dcalc,j are observed and calculated data, respectively, and σj represents the data errors for all N data points. A value of this metric close to 1 is commonly interpreted as a data fit within the range of observational error, that is, in the case of this study the error floors, which somewhat reduces the statistical significance. For the final inversion model, 112 iterations were required to fit the full impedance tensor and the tipper data, with an nRMS of 1.8. As shown in Fig. 6, this value is strongly influenced by very few large-residual sites, as is characteristic for a non-robust metric of the least-squares type. The fit to the majority of the observations becomes much clearer from the data fit curves and detailed maps of nRMS, which were calculated for each data components and for each frequency band at each station, are presented in Supporting Information Sections S1 and S2, respectively. The 3-D inversion results are summarized in Figs 7 and 8, and discussed in the following section. Figure 6. View large Download slide Map of nRMS misfit values for each station. A total nRMS of 1.8 was achieved. Maps of variation in nRMS for each frequency band and for all the elements of the impedance tensor are shown in Supporting Information Section S2. Figure 6. View large Download slide Map of nRMS misfit values for each station. A total nRMS of 1.8 was achieved. Maps of variation in nRMS for each frequency band and for all the elements of the impedance tensor are shown in Supporting Information Section S2. Figure 7. View largeDownload slide Horizontal slices of the final resistivity model at various depth. Panels (a) and (b) depict the depth extent of the shallow conductor beneath the fumarolic field on the northern edge of Furnas Lake. Panels (c) and (d) indicate the limits of the main conductor, C2, beneath the inner caldera at Furnas. Figure 7. View largeDownload slide Horizontal slices of the final resistivity model at various depth. Panels (a) and (b) depict the depth extent of the shallow conductor beneath the fumarolic field on the northern edge of Furnas Lake. Panels (c) and (d) indicate the limits of the main conductor, C2, beneath the inner caldera at Furnas. Figure 8. View largeDownload slide Extracted 2-D slices from our final 3-D resistivity model. (a) P1 is an NS profile passing through the main fumarolic region on the northern edge of Furnas Lake towards the most recent site of volcanism, the 1630 Dome, (b) an EW profile through the centre of the survey imaging beneath Furnas Lake and Pico do Gaspar and (c) a profile passing through the centres of the main domes in the inner caldera, Gaspar Dome and the 1630 Dome in the south of the profile. Figure 8. View largeDownload slide Extracted 2-D slices from our final 3-D resistivity model. (a) P1 is an NS profile passing through the main fumarolic region on the northern edge of Furnas Lake towards the most recent site of volcanism, the 1630 Dome, (b) an EW profile through the centre of the survey imaging beneath Furnas Lake and Pico do Gaspar and (c) a profile passing through the centres of the main domes in the inner caldera, Gaspar Dome and the 1630 Dome in the south of the profile. 7 RESULTS Depth slices through the 3-D resistivity inversion volume are shown in Fig. 7 and selected cross-sections through various regions of interest are shown in Fig. 8. The overall structure comprises a resistive surface layer (R1) of varying thickness with resistivities in the region of 100–1000 Ωm. It overlays zones of high conductivity (C1 and C2) beneath the main caldera, ranging from depths of 100 m to a general depth of at approximately 500 m. Almost all of the MT soundings recorded in the campaign are characterized by relatively high near-surface resistivity, in the region of 100–200 Ωm, with the exception of one site (namely FUR004, Supporting Information Section S1), adjacent to the main fumarolic fields, with a resistivity range in the order of tens of Ωm at the near surface. The shallow subsurface of the entire area within this survey (R1) is caldera infill, dominated by material derived from volcanic activity and subsequent erosion. The deposits include coarse lapilli, fine ashes and pumice blocks in the upper deposits, however due to the lack of any borehole data the full extent of this stratigraphy is unclear. The pore space is filled with water of meteoric origin, probably similar to the water within Furnas Lake (≈60 Ωm) (Andrade et al. 2016). A shallow, very localized, conductive zone (C1) is found in the northern region of Furnas Lake, at ≈100 m beneath the surface, adjacent to a mapped fault, that is steeply dipping (≈70°), known as Furnas Fault 2 (FF2, see Carmo 2013; Fig. 5). This conductive zone extends vertically ≈100 m, with a resistivity of 0.5–10 Ωm, which explains the behaviour of the high-frequency induction vectors, indicating a concentration of conductive material at this location and depth (Figs 5b and c). This zone may not be completely controlled by FF2 as the orientation of the lens like structure is not parallel to FF2 strike. Geoelectric strike direction estimates, derived from our WALDIM analysis, infer ranges of 40°–55° (Fig. 4), however, the strike of FF2 is ≈100° (Carmo et al. 2015). Note that the orientation of C1 is almost parallel to the steep caldera wall which may suggest a relation to the limits of caldera infill, or tectonic features related to the caldera wall. We interpret this conductive zone as a shallow, cap-like structure that partially seals the ascent of vapour feeding the surface hydrothermal manifestations directly above, possibly influenced/controlled by the fault, FF2. Fig. 9 compares the resistivity sections obtained by AMT and ERT methods approximately along an NE–SW, 1.5 km-long profile (P4) along with the soil CO2 degassing results from Viveiros et al. (2010). The corresponding ERT profile transects the FF2 structure on the northern edge of Furnas Lake. Data were acquired using a multielectrode ABEM system (Terrameter LS) using a roll-along Wenner-Schlumberger configuration (64 electrodes with a spacing of 10 m) with an injected current of 100–400 mA. The resulting apparent resistivity data including topography information were then inverted using the RES2DINV® software (Loke & Barker 1996; Loke 2016). Both geophysical methods clearly and coherently image the shallow conductive zone C1, that is found directly beneath the highest concentration of CO2 degassing in the region of 400–900 m on P4. This illustrates the consistency of using both electrical and EM geophysical methods in imaging subsurface features and isolating pathways beneath high CO2 degassing output in this part of the caldera. Due to the lack of coverage of data around Furnas village, which also exhibits localized zones of high CO2 degassing, we cannot conclude if there is a distinct network of shallow conductive zones distributed beneath the village feeding other zones of degassing. Figure 9. View largeDownload slide Comparison of AMT and ERT results along a profile P-4 (Fig. 8) across the fumarole area. Top: Satellite image of surface and topography, overlain by the CO2 flux distribution from Viveiros et al. (2010). The dashed white line marks the location of the main fault in the area, FF2. Centre: inversion results for an ERT profile starting at the NW rim of the lake following a roughly NNE direction. Bottom: slice through the 3-D AMT model approximately along the ERT profile above. Both methods agree well in their overlap. Figure 9. View largeDownload slide Comparison of AMT and ERT results along a profile P-4 (Fig. 8) across the fumarole area. Top: Satellite image of surface and topography, overlain by the CO2 flux distribution from Viveiros et al. (2010). The dashed white line marks the location of the main fault in the area, FF2. Centre: inversion results for an ERT profile starting at the NW rim of the lake following a roughly NNE direction. Bottom: slice through the 3-D AMT model approximately along the ERT profile above. Both methods agree well in their overlap. Beneath the localized conductor C1 we find a thin, more resistive region (R2) of approximately 100 m thickness that increases in resistivity with depth. It is approximately at sea level and separates C1 from the dominant conductive zone beneath the survey area, C2. This large conductive zone has an upper boundary at 150 m b.s.l., approximately 500–600 m beneath the surface. It has a very conductive upper layer in the region of 1–10 Ωm, and lies directly beneath C1, however the upper boundary of C2 varies spatially in our inversion, but this may result from poor data coverage especially in the region of Furnas Lake. The highest conductivity found within C2 lies directly beneath C1, which could indicate they are linked, however sensitivity analysis of our 3-D model in this region could not confirm that C1 and C2 are separate zones. Unfortunately, there are no borehole data in the Furnas region so is not possible to confirm or refute this hypothesis. The contrast between C2 and R2 can be clearly seen in 2-D slice, P1, that is orientated north-south along the eastern edge of Furnas Lake through our survey area (Fig. 8a). It shows C2 thinning towards the south and becoming more resistive but interestingly shows no conductive body beneath the 1630 Dome, the remnant of the most recent event of volcanism. C2 is clearly visible beneath Gaspar Dome and can be seen in both Profiles P2 (E–W across the survey area, Fig. 8b) and P3 (passing through the two volcanic domes, Fig. 8c). The top of this anomalous feature is found at a depth of ≈250 m b.s.l., that is, ≈600 m beneath the surface, with a thickness of at least 500 m. We cannot confirm the eastern boundary of C2, however C2 appears to connect the main fumarolic fields at Furnas with Gaspar Dome, but thins towards the south. The resolution of the resulting 3-D model is significantly reduced ≈1 km below the surface, underneath the C2 anomaly, though it cannot be excluded that this is due to the shielding effects of the high-conductivity zones above. The absence of C2 in the south of the survey implies that saline intrusion from the Atlantic Ocean is not being sensed by our high-frequency MT survey at these depths. Earlier modelling studies, prior to analysis and inversion, confirmed the ocean effect to be absent at these high-frequencies. 8 DISCUSSION The electrical conductivity of the subsurface is a crucial parameter in the characterization of volcanic and hydrogeothermal settings and has been exploited widely by MT surveys in various settings around the world in recent decades (Muñoz 2014, and references therein). Regions of high conductivity within the shallow Earth, can be attributed to partial melt related to magma reservoirs, the presence of graphite, metallic material, or conducting fluids. Magma reservoirs, in particular in view of the high melt fraction required, are highly improbable in the depth range found in this study as there is no volcanological evidence pointing to this explanation. Carbon deposits in the form of graphite conductors are equally improbable because the depth ranges of the conductive zones do not favour this explanation, in particular when the required large volume of connected conducting material is considered. Franco (2015) reports disseminated pyrite found in boreholes in a geothermal system at the northern flank of the near-by Fogo (Água de Pau) Volcano, ranging over the whole depth of the geothermal system, though no estimates of their quantity are given. The contribution of pyrite to the conduction process will be minor, but cannot be excluded. For hydrothermal systems, several possible explanations for the high conductivities exist: hydrothermal systems may appear conductive because of the fluids involved (which would need to be highly saline), or the alteration products caused by recent and ancient fluid movements, for example, a conductive clay cap overlaying the fluid system itself. Both may occur side by side in a given hydrothermal system. The resistivity structure determined by our 3-D inversion of high-frequency MT data, maps the structure(s) of the hydrothermal system beneath the Furnas Caldera and we will discuss this hypothesis in detail in the following. Our 3-D inversion of the upper hydrothermal system at Furnas yielded two highly conductive zones C1 and C2, respectively. Though a robust feature of our models, sensitivity studies of locally perturbed models show that we cannot be certain that they are fully separated. The data fit is not substantially changed by introducing a conductive connection between the two regions, however, their representative resistivity characteristics at different depths imply a complex system of pathways, which feeds the fumaroles at the surface above. For the purpose of this discussion we will treat them as separate zones. In volcanic areas, a conductor with resistivities in the region of 0.5–10 Ωm has been attributed to the presence of hydrothermal water or alteration minerals (Ussher et al. 2000; Muñoz 2014), however both involve contrasting permeability structures to be present. A layer of smectite clay typically exists around the reservoir in a geothermal system. Smectite and smectite-illite clays are the predominant products derived from rock alteration processes that form in the 50–200 °C zone over and adjacent to volcanic hosted systems (Browne 1978). These clay-rich zones are important components for geothermal hydrology, quite often partially sealing the reservoir, where high permeability is required for hydrothermal fluids to migrate within the system, usually through an extensive network of fractures. Typically, values of 1–10 Ωm are observed at clay-rich zones, whereas for the adjacent hydrothermal reservoir, conductivity ranges between 5 and 100 Ωm are believed to be typical (Wright et al. 1985; Pellerin et al. 1996). Thus, the hottest regions of a geothermal system are often characterized by higher resistivity than observed in zones of hydrothermal alteration (Ussher et al. 2000; Muñoz 2014). Geothermal manifestations (fumaroles) are found at the surface directly above the high conductivity zones C1 and C2 at the northern edge of Furnas Lake. Recent work by Woitischek et al. (2017), characterizing the origin of hydrothermal waters on São Miguel by their chemical and isotopic composition, confirm that the hydrothermal waters found at Furnas are generally of meteoric origin. At the surface they are slightly saline with resistivities near 10 Ωm, while the resistivity of the water within Furnas Lake is about 60 Ωm (Andrade et al. 2016). It is unknown how deep the meteoric waters migrate and which role they play in the full system. In the case of C1, the presence of clay material (smectite) is required along with the presence of water to account for the high conductivity values, sensed at these depths between 100 and 200 m beneath the surface. Petrophysical reasoning based on the work of Revil et al. (2002) quantifying the electrical properties of altered volcaniclastic materials can be used to estimate the percentage of smectite needed to account for a conductivity value, once an estimate of the temperature is known. Fig. 10 shows the range of smectite required to account for the conductivity of C1 for a broad range of temperatures and fluid conductivities. The petrophysical model assumed is described in Appendix. From the results we conclude that a mass fraction of approximately 20 per cent smectite is sufficient to explain the observed bulk conductivity of C1, where the remaining 80 per cent comprise fluid and the unaltered original volcaniclastic material. Figure 10. View largeDownload slide Calculation of smectite content based on the studies of Waxman & Smits (1968), Revil et al. (1998, 2002) and Flóvenz et al. (2005) for the temperature range 80–200 °C. The conductivity of the shallow conductor C1 (upper grey zone) structure cannot be explained by ionic conduction in the fluid assuming Archie's law (lower grey zone) even when comparatively high temperatures and salinities are used. Taking into account the surface conductivity of smectite allows us to explain the behaviour of the C1 zone assuming a moderate mass fraction of this clay mineral. An explanation of the underlying estimate is given in Appendix. Figure 10. View largeDownload slide Calculation of smectite content based on the studies of Waxman & Smits (1968), Revil et al. (1998, 2002) and Flóvenz et al. (2005) for the temperature range 80–200 °C. The conductivity of the shallow conductor C1 (upper grey zone) structure cannot be explained by ionic conduction in the fluid assuming Archie's law (lower grey zone) even when comparatively high temperatures and salinities are used. Taking into account the surface conductivity of smectite allows us to explain the behaviour of the C1 zone assuming a moderate mass fraction of this clay mineral. An explanation of the underlying estimate is given in Appendix. Unfortunately, there are no local borehole data to constrain our observations, however borehole data from geothermal exploration at the NW flank of Fogo Volcano (Ribeira Grande, approximately 20 km away from our study area) may provide some clues, even though the morphological position and hydrodynamical characteristic of that hydrothermal system differs from Furnas, the geochemical composition was found to be quite similar (Caliro et al. 2015). The work of Franco (2015) at Fogo (Água de Pau), identified two zones of hydrothermal alteration beneath a shallow zone of relatively unaltered rocks, based on borehole results. Evidence from this study documents a smectite layer at comparable depth to the conductive zone C1, that is observed in our study and the author assigns a temperature range 200 °C to 220 °C. This observation complements our petrophysically based estimates that confirm the requirement of smectite to yield the high conductivity of C1. We propose that the conductivity of C1 is controlled by the presence of smectite clay material, forming a partial seal of a localized vapour driven aquifer, at a temperature of up to 200 °C. Near Furnas Lake it is perturbed by the fault zone FF2. The horizontal resistivity slices in Fig. 7 show that at shallow depths, the strike of the anomalous zone deviates from FF2, with a significant clockwise rotation of the axis at larger depth when we reach C2. No low resistivity zones are found along the eastern continuation of this fault zone. This structure is consistent with the idea that the continuous formation of the low-permeability seal zone inhibits direct ascent, with preferential fluid or steam pathways changing with depth. Beneath C1 we find the thin, more resistive zone R2. Sensitivity analysis on the 3-D modelling cannot rule out that both C1 and the deeper, more extensive conductive zone C2 are connected. The resistivity change may be explained by a change in clay alteration mineralogy that can also act as a proxy for changes in temperature. The study of Arnason et al. (2008) documents that in the temperature range from 220 °C to 240–250 °C, low temperature zeolites gradually disappear and smectite is transformed into chlorite, yielding a mixed-clay zone, where smectite and chlorite coexist. Thus, R2 probably represents a mixed clay zone. This zone may exhibit a wide range of resistivities (Arnason et al. 2008; Hersir & Árnason 2015), from 10 Ωm to more than 1000 Ωm. The resistivity is strongly influenced by residual smectite content, less by porosity and pore/fracture characteristics for a low-salinity system as assumed in this study. Weisenberger et al. (2016) recently have shown that surface conduction, which critically depends on the cation exchange capacity (CEC, see Appendix), is considerably lower in chlorite-rich zones. The emergence of C2, the dominant conductive zone beneath Furnas may indicate the top of the chlorite-dominated zone, where all of the smectite has disappeared from the overlying mixed clay zone. Arnason et al. (2008) shows that the chlorite zone commences at 250 °C, and that at higher temperatures 260–270 °C, epidote becomes the dominant mineral in the so-called chlorite-epidote zone, where conductivities may greater than 100 Ωm (Ussher et al. 2000). However, these temperature estimates apply to freshwater systems, while brine fluids would shift the zoning to higher values. We conclude here that the more recent temperature estimates of ≈275 °C from Caliro et al. (2015) for the hydrothermal fluids feeding Furnas Lake fumaroles, are associated with the deeper reservoir C2. Ocean water (0.3 Ωm) can be eliminated as a source of the high conductivity within the system from the geochemical signature found by Woitischek et al. (2017). Furthermore, the observed conductive zone C2 does not extend towards the south, and we do not find a connection to the ocean there. One may speculate that the lack of southward connection originates from the recent explosive volcanism within the caldera, which has been migrating from the NE to SW since at least 1.9 ka ago (Guest et al. 2015) and may have destroyed earlier fluid pathways. It has to be noted, however, that ocean water is not the only potential source of salinity as the fluid associated with the C2 zone may derive from mixing with highly saline magma-derived fluids (Nicholson 1993; Yardley & Bodnar 2014). However, from our measurements it is not possible to further constrain fluid composition. If we assume that electrical conduction occurs mainly through the fluid phase filling the reservoir, we need the hydrothermal system to be in fluid phase. Fig. 11 shows vertical resistivity profiles extracted from the 3-D model for the sites nearest to the fumarole area, together with the boiling curves for different values of salinities, calculated using the method of Haas (1971) as modified by Canet et al. (2011). The position of C2 (centre) in the resistivity profiles at ≈500 m beneath the surface would imply a temperature of 250 °C, independent of the assumed salinity. This agrees well with the thermal structure found at the Ribeira Grande system (Rangel 2014; Franco 2015), where similar temperatures of about 240 °C were found. The recent geothermometric estimate of 275 °C (Caliro et al. 2015) is significantly higher than the value derived here, but still falls into the depth range of the C2 conductor. Figure 11. View largeDownload slide Boiling curves for different salinities between 1 and 8 wt per cent, calculated using the method described in Haas (1971) and Canet et al. (2011) for hydrostatic conditions and near-adiabatic ascent. Ocean water has an average salinity value of 3.5 wt per cent. The orange curves are vertical profiles at the sites nearest to the fumarole area, extracted from the final 3-D modelling results. The vertical blue line marks the recent temperature estimate of ≈270 °C from Caliro et al. (2015). Assuming that the conductor C2 is related to the fluid-steam interface, an estimate of reservoir temperature at about 250 °C could be derived from our models. This agrees notably well with the temperatures found in the geothermal wells in the Cachaços-Lombadas sector of the Ribeira Grande geothermal field at the NW flank of the Fogo (Água de Pau) Volcano. Note that temperatures close to the surface, where mixing with meteoric waters is important, may be considerably lower than the boiling curves suggest. Figure 11. View largeDownload slide Boiling curves for different salinities between 1 and 8 wt per cent, calculated using the method described in Haas (1971) and Canet et al. (2011) for hydrostatic conditions and near-adiabatic ascent. Ocean water has an average salinity value of 3.5 wt per cent. The orange curves are vertical profiles at the sites nearest to the fumarole area, extracted from the final 3-D modelling results. The vertical blue line marks the recent temperature estimate of ≈270 °C from Caliro et al. (2015). Assuming that the conductor C2 is related to the fluid-steam interface, an estimate of reservoir temperature at about 250 °C could be derived from our models. This agrees notably well with the temperatures found in the geothermal wells in the Cachaços-Lombadas sector of the Ribeira Grande geothermal field at the NW flank of the Fogo (Água de Pau) Volcano. Note that temperatures close to the surface, where mixing with meteoric waters is important, may be considerably lower than the boiling curves suggest. The complex conductivity image of the Furnas Caldera can possibly be explained by the following hypothetical process: an extended hydrothermal system within a large caldera develops, leading to the development of the typical associated alteration zones (clay minerals). Within this caldera, a succession of smaller eruptions occur (e.g. Gaspar and 1630 Domes), which perturb the pre-existing hydrothermal system locally, and add to the sedimentary filling of the original caldera. Subsequently, a new generation of hydrothermal systems develop, destroying or replacing the older ones. As pointed out by Guest et al. (2015), the area of the Furnas Caldera has a documented history of frequent eruptions since the formation of the Povoção Caldera at least 30 ka ago. 9 CONCLUSIONS This MT study has yielded, for the first time, a 3-D electrical resistivity model of the shallow hydrothermal system beneath Furnas Volcano (Azores, Portugal) and has been used to create a revised conceptual model of this dynamic area. The inversion of high frequency data identified highly conductive zones (0.5–10 Ωm) beneath the main central caldera, that feed the surface manifestations of volcano-hydrothermal activity that are characterized by boiling temperature fumaroles, steaming ground, thermal springs, and areas of intense diffuse degassing of CO2 and 222Rn. As no borehole data or temperature profile for this area is available, our petrophysical reasoning (see Appendix) led to the conclusion that the shallowest conductive zone (≈100 m beneath the surface) comprise at least 20 per cent of smectite. This high clay fraction could partially seal the deeper reservoir comprising meteoric water feeding the fumaroles on the surface above. New insights into the structural regime controlling this shallow conductive aquifer have been achieved. The geoelectrical strike of the shallow conductor is not aligned with the strike of the main fault in the area, which was previously interpreted to directly control the orientation of the feeding aquifers. The geoelectrical strike of the conductor is parallel to the steep caldera wall adjacent to it, which may infer that morphology of the caldera, especially beneath the infilling from various phases of caldera collapses and subsequent erosion is a factor in its structural character, however the main fault acts as a pathway from this aquifer towards the surface. The deeper, more dominant conductor imaged by MT outlines the broader extent of the Furnas hydrothermal system within the caldera. The extent of this conductive zone has not been fully mapped by our study due to access issues adjacent to Furnas village and the presence of a caldera lake. This conductor however does not extend towards the south, implying it is not associated with salt water intrusion from the Atlantic Ocean, confirming the recent geochemical analysis of the waters at Furnas to be mainly meteoric in origin. This conductive zone extends beneath Gaspar Dome but is not evident beneath the most recent centre of volcanism, the 1630 Dome, implying the current hydrothermal system is a result of various volcanic events. Borehole and thermal gradient data from geothermal exploration studies at Fogo (Água de Pau), an adjacent volcano, provide clues in characterizing the hydrothermal system in our study. We infer the deeper conductive zone beneath Furnas is dominated by chlorite, inferring a temperature of at least 250 °C. Imaging the proposed magma reservoir that fuels the volcano-hydrothermal system at Furnas, is not within the scope of this study. The inclusion of longer period MT data, combining topography and bathymetry data in the inversion, will be required to successfully image this deeper structure. ACKNOWLEDGEMENTS This multiyear study was primarily financed by DIAS. We appreciate the funding granted to AR (ISTerre) through the INSU project ‘Improving imaging of hydrothermal systems by joint inversion of ERT-AMT data. Case study of Furnas Caldera (Azores, Portugal)’. The fieldwork would not have been possible without the help of Ernesto Sousa, Victor Sousa, Bruno Medeiros, Lucía Rodríguez, Joana Eleutério, Annika Rödder, Damien Jougnot and Niels Grobbe. The Direção Regional do Ambiente, Serviços de Recursos Hídricos e Ordenamento do Território (Ponta Delgada) is thanked for providing the high-resolution bathymetry for the Furnas Lake. The tectonic maps were produced with the data from Azores Regional Government through a PhD grant to RC from Fundo Regional da Ciência e Tecnologia, and through Serviço Regional de Proteção Civil e Bombeiros dos Açores in the scope of the scientific and technical protocols to guarantee the Azores Seismovolcanic Surveillance and the Emergency Planning Studies of Centro de Informação e Vigilância Sismovulcânica dos Açores (CIVISA). The Irish Centre for High Performance Computing (ICHEC) is thanked for availing the Fionn cluster to carry out the numerical computations. Finally, G. Egbert, A. Kelbert, and N. Meqbel are gratefully thanked for making their ModEM code available to us. We thank Léa Levy for her comments on temperature stability range of secondary minerals, in particular smectite. Last, but not least, we thank the two reviewers, J. Hübert and J. Peacock, for their constructive and helpful comments. REFERENCES Aizawa K. et al.  , 2005. Hydrothermal system beneath Mt. Fuji volcano inferred from magnetotellurics and electric self-potential, Earth planet. Sci. Lett.  235( 1-2), 343– 355. https://doi.org/10.1016/j.epsl.2005.03.023 Google Scholar CrossRef Search ADS   Aizawa K., Ogawa Y., Ishido T., 2009. Groundwater flow and hydrothermal systems within volcanic edifices: Delineation by electric self-potential and magnetotellurics, J. geophys. Res.  114( B1), B01208. https://doi.org/10.1029/2008JB005910 Andrade C., Viveiros F., Cruz J.V., Coutinho R., Silva C., 2016. Estimation of the CO2 flux from Furnas volcanic lake (Sao Miguel, Azores), J. Volcanol. Geotherm. Res.  315 51– 64. https://doi.org/10.1016/j.jvolgeores.2016.02.005 Google Scholar CrossRef Search ADS   Archie G.E., 1942. The electrical resistivity log as an aid in determining some reservoir characteristics, Trans. Am. Inst. Min. Metall. Pet. Eng.  146 54– 62. Arnason K., Karlsdottir R., Eysteinsson H., Flovenz O.G., Gudlaugsson S.T., 2008. The resistivity structure of high-temperature geothermal systems in iceland, in Short Course III on Exploration for Geothermal Resources, Lake Naivasha, Kenya, October 24 - November 17  UNU-GTP, (http://www.unugtp.is, last accessed 15 September 2017). Beamish D., Travassos J., 1992. The use of the D+ solution in magnetotelluric interpretation, J. Appl. Geophys.  29( 1), 1– 19. https://doi.org/10.1016/0926-9851(92)90009-A Google Scholar CrossRef Search ADS   Blanco I., García A., Torta J.M., 1997. Magnetic study of the Furnas volcano (Azores), Ann. Geofis.  XL 341– 358. Browne P.R.L., 1978. Hydrothermal alteration in active geothermal fields, Annu. Rev. Earth Planet. Sci.  6 229– 250. https://doi.org/10.1146/annurev.ea.06.050178.001305 Google Scholar CrossRef Search ADS   Buforn E., Udías A., 2010. Azores-Tunisia, a tectonically complex plate boundary, in Advances in Geophysics  vol. 52, chap. 3, pp. 139– 182, ed. Dmowska R., Elsevier. Google Scholar CrossRef Search ADS   Cagniard L., 1953. Basic theory of the magneto-telluric method of geophysical prospecting, Geophysics  18 605– 635. https://doi.org/10.1190/1.1437915 Google Scholar CrossRef Search ADS   Caldwell T.G., Bibby H.M., Brown C., 2004. The magnetotelluric phase tensor, Geophys. J. Int.  158 457– 469. https://doi.org/10.1111/j.1365-246X.2004.02281.x Google Scholar CrossRef Search ADS   Caliro S., Viveiros F., Chiodini G., Ferreira T., 2015. Gas geochemistry of hydrothermal fluids of the S. Miguel and Terceira Islands, Azores, Geochim. Cosmochim. Acta  168 43– 57. https://doi.org/10.1016/j.gca.2015.07.009 Google Scholar CrossRef Search ADS   Camacho A.G., Montesinos F.G., Vieira R., 1997. A three-dimensional gravity inversion applied to Sao Miguel Island (Azores), J. geophys. Res.  102 7717– 7730. https://doi.org/10.1029/96JB03667 Google Scholar CrossRef Search ADS   Canet C., Franco S.I., Prol-Ledesma R.M., González-Partida E., Villanueva-Estrada R.E., 2011. A model of boiling for fluid inclusion studies: application to the Bolaños Ag-Au-Pb-Zn epithermal deposit, Western Mexico, J. Geochem. Explor.  110 118– 125. https://doi.org/10.1016/j.gexplo.2011.04.005 Google Scholar CrossRef Search ADS   Carmo R., 2013. Estudos de neotectónica na ilha de S. Miguel, uma contribuição para o estudo do risco sísmico no arquipélago dos Açores., PhD thesis  University of the Azores. Carmo R., Madeira J., Ferreira T., Queiroz G., Hipólito A., 2015. Volcano-tectonic structures of São Miguel Island, Azores, in Volcanic Geology of S. Miguel Island (Azores Archipelago)  chap. 6, pp. 65– 86, eds Gaspar J.L., Guest J.E., Duncan A.M., Barriga F., Chester D.K., Geological Society of London. Google Scholar CrossRef Search ADS   Carvalho M., Forjaz V., Almeida C., 2006. Chemical composition of deep hydrothermal fluids in the Ribeira Grande geothermal field (São Miguel, Azores), J. Volcanol. Geotherm. Res.  156( 1-2), 116– 134. https://doi.org/10.1016/j.jvolgeores.2006.03.015 Google Scholar CrossRef Search ADS   Chave A., Jones A., 2012. The Magnetotelluric Method: Theory and Practice  Cambridge Univ. Press. Google Scholar CrossRef Search ADS   Cole P., Queiroz G., Wallenstein N., Gaspar J., Duncan A., Guest J., 1995. An historic subplinian/phreatomagmatic eruption: the 1630 AD eruption of Furnas volcano, São Miguel, Azores, J. Volcanol. Geotherm. Res.  69( 1–2), 117– 135. https://doi.org/10.1016/0377-0273(95)00033-X Google Scholar CrossRef Search ADS   Cruz J.V., França Z., 2006. Hydrogeochemistry of thermal and mineral water springs of the Azores archipelago (Portugal), J. Volcanol. Geotherm. Res.  151( 4), 382– 398. https://doi.org/10.1016/j.jvolgeores.2005.09.001 Google Scholar CrossRef Search ADS   Egbert G.D., Kelbert A., 2012. Computational recipes for electromagnetic inverse problems, Geophys. J. Int.  189( 1), 251– 267. https://doi.org/10.1111/j.1365-246X.2011.05347.x Google Scholar CrossRef Search ADS   Flóvenz O.G., Spangenberg E., Kulenkampff J., Árnason K., Karlsdóttir R., Huenges E., 2005. The role of electrical interface conduction in geothermal exploration, in Proceedings World Geothermal Congress , Antalya, Turkey, 24– 29 April 2005. Forjaz V., 1996. Geothermal reservoirs of Furnas Laboratory Volcano, in 2nd Workshop on European Laboratory Volcanoes  Santorini, Greece. Franco A., 2015. Subsurface geology and hydrothermal alteration of Cachaços-Lombadas sector, Ribeira Grande geothermal field, (São Miguel island, Azores), Tech. Rep. UNU-GTP-2015-10, United Nations University, Geothermal Training Programme, (http://www.unugtp.is, last accessed 15 September 2017). Guest J., Pacheco J.M., Cole P.D., Duncan A.M., Wallenstein N., Queiroz G., Gaspar J.L., Ferreira T., 2015. The volcanic history of Furnas Volcano, São Miguel, Azores, in Volcanic Geology of São Miguel Island (Azores Archipelago)  chap. 6, pp. 125– 134, eds Gaspar J.L., Guest J.E., Duncan A.M., Barriga F., Chester D.K., Geological Society of London. Google Scholar CrossRef Search ADS   Guest J.E., Gaspar J.L., Cole P.D., Queiroz G., Duncan A.M., Wallenstein N., Ferreira T., Pacheco J.-M., 1999. Volcanic geology of Furnas Volcano, São Miguel, Azores, J. Volcanol. Geotherm. Res.  92( 1–2), 1– 29. https://doi.org/10.1016/S0377-0273(99)00064-5 Google Scholar CrossRef Search ADS   Haas J.L., 1971. The effect of salinity on the maximum thermal gradient of a hydrothermal system at hydrostatic pressure, Econ. Geol.  66 940– 946. https://doi.org/10.2113/gsecongeo.66.6.940 Google Scholar CrossRef Search ADS   Hersir G., Árnason K., 2015. Resistivity of rocks, in Short Course X on Exploration for Geothermal Resources  organized by UNU-GTP, GDC and KenGen, at Lake Bogoria and Lake Naivasha, Kenya, 2015 November 9–December 1, (http://www.unugtp.is, last accessed 15 September 2017). Hoover D., Amaral R., Broker M., 1983. Preliminary report on audio-magnetotelluric survey on São Miguel island, Azores, Portugal., Open-File Report 82-441, USGS, Menlo Park CA, USA. Hoover D., Rodrigues Da Silva A., Pierce H.A., Amaral R., 1984. Application of audio-magnetotelluric surveys on São Miguel island, Azores, Portugal, in Transactions of the Geothermal Research Council , no. 8, pp. 499– 503, Menlo Park CA, USA. Jeffery A.J., 2016. Petrogenesis and contrasting eruption styles of peralkaline silicic magmas from Terceira and São Miguel, Azores, PhD thesis , Keele University, UK. Jones A., Jödicke H., 1984. Magnetotelluric transfer function estimation improvement by a coherence-based rejection technique, in SEG Technical Program Expanded Abstracts 1984 , pp. 51– 55. Kelbert A., Meqbel N., Egbert G.D., Tandon K., 2014. ModEM: A modular system for inversion of electromagnetic geophysical data, Comput. Geosci.  66 440– 53. https://doi.org/10.1016/j.cageo.2014.01.010 Google Scholar CrossRef Search ADS   Loke M., 2016. ‘Tutorial: 2-D and 3-Delectrical imaging surveys’. Available at: http://www.geotomosoft.com, last accessed 15 September 2017. Loke M., Barker R.D., 1996. Rapid least-squares inversion of apparent resistivity pseudosections by a quasi-Newton method, Geophys. Prospect.  44( 1), 131– 152. https://doi.org/10.1111/j.1365-2478.1996.tb00142.x Google Scholar CrossRef Search ADS   Machado F., 1959. Submarine pits of the Azores plateau, Bull. Volcanol.  21 109– 116. https://doi.org/10.1007/BF02596510 Google Scholar CrossRef Search ADS   Martí A., Queralt P., Jones A.G., Ledo J., 2005. Improving Bahr's invariant parameters using the WAL approach, Geophys. J. Int.  163 38– 41. https://doi.org/10.1111/j.1365-246X.2005.02748.x Google Scholar CrossRef Search ADS   Martí A., Queralt P., Ledo J., 2009. WALDIM: A code for the dimensionality analysis of magnetotelluric data using the rotational invariants of the magnetotelluric tensor, Comput. Geosci.  35( 12), 2295– 2303. https://doi.org/10.1016/j.cageo.2009.03.004 Google Scholar CrossRef Search ADS   Martí A., Queralt P., Ledo J., Farquharson C., 2010. Dimensionality imprint of electrical anisotropy in magnetotelluric responses, Phys. Earth planet. Inter.  182( 3-4), 139– 151. https://doi.org/10.1016/j.pepi.2010.07.007 Google Scholar CrossRef Search ADS   Meqbel N.M., 2009. The electrical conductivity structure of the Dead Sea Basin derived from 2D and 3D inversion of magnetotelluric data, PhD thesis , Freie Universität Berlin, Germany. Miranda J.M., Luis J.F., Abreu I., Victor L.A.M., Galdeano A., Rossignol J.C., 1991. Tectonic framework of the Azores Triple Junction, Geophys. Res. Lett.  18( 8), 1421– 1424. https://doi.org/10.1029/91GL01607 Google Scholar CrossRef Search ADS   Miranda J.M., Luis J.F., Lourenço N., Fernandes R.M.S., 2015. The structure of the Azores triple junction: implications from São Miguel Island, in Volcanic Geology of São Miguel Island (Azores Archipelago), chap. 2 , pp. 5– 13, eds Gaspar J.L., Guest J.E., Duncan A.M., Barriga F., Chester D.K., Geological Society of London. Google Scholar CrossRef Search ADS   Montesinos F.G., Camacho A.G., Vieira R., 1999. Analysis of gravimetric anomalies in Furnas caldera (São Miguel, Azores), J. Volcanol. Geotherm. Res.  92( 1-2), 67– 81. https://doi.org/10.1016/S0377-0273(99)00068-2 Google Scholar CrossRef Search ADS   Moore J.N., Allis R.G., Necmčok M., Powell T.S., Bruton C.J., Wannamaker P.E., Raharjo I.B., Normanc D.I., 2008. The evolution of volcano-hosted geothermal systems based on deep wells from Karaha-Telaga Bodas, Indonesia, Am. J. Sci.  308 1– 48. https://doi.org/10.2475/01.2008.01 Google Scholar CrossRef Search ADS   Moore R.B., 1990. Volcanic geology and eruption frequency, São Miguel, Azores, Bull. Volcanol.  52 602– 614. https://doi.org/10.1007/BF00301211 Google Scholar CrossRef Search ADS   Muñoz G., 2014. Exploring for geothermal resources with electromagnetic methods, Surv. Geophys.  35( 1), 101– 122. https://doi.org/10.1007/s10712-013-9236-0 Google Scholar CrossRef Search ADS   Nicholson K., 1993. Geothermal Fluids  Springer-Verlag, Berlin Heidelberg, ISBN: 978-3-642-77846-9. Google Scholar CrossRef Search ADS   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. https://doi.org/10.1186/BF03352120 Google Scholar CrossRef Search ADS   Parker R.L., 1980. The inverse problem of electromagnetic induction: existence and construction of solutions based on incomplete data, J. geophys. Res.  85( 8), 4421– 4428. https://doi.org/10.1029/JB085iB08p04421 Google Scholar CrossRef Search ADS   Parker R.L., Booker J.R., 1996. Optimal one-dimensional inversion and bounding of magnetotelluric apparent resistivity and phase measurements, Phys. Earth planet. Inter.  98( 3-4), 269– 282. https://doi.org/10.1016/S0031-9201(96)03191-3 Google Scholar CrossRef Search ADS   Parker R.L., Whaler K., 1981. Numerical methods for establishing solutions to the inverse problem of electromagnetic induction, J. geophys. Res.  86 9574– 9584. https://doi.org/10.1029/JB086iB10p09574 Google Scholar CrossRef Search ADS   Parkinson W., 1962. The influence of continents and oceans on geomagnetic variations, Geophys. J. R. astr. Soc.  6 441– 449. https://doi.org/10.1111/j.1365-246X.1962.tb02992.x Google Scholar CrossRef Search ADS   Pedone M., Viveiros F., Aiuppa A., Giudice G., Grassa F., Gagliano A.L., Francofonte V., Ferreira T., 2015. Total (fumarolic + diffuse soil) CO2 output from Furnas volcano, Earth Planets Space  67 doi:10.1186/s40623-015-0345-5. https://doi.org/10.1186/s40623-015-0345-5 Google Scholar CrossRef Search ADS   Pellerin L., Johnston J.M., Hohmann G.W., 1996. A numerical evaluation of electromagnetic methods in geothermal exploration, Geophysics  61( 01), 121– 130. https://doi.org/10.1190/1.1443931 Google Scholar CrossRef Search ADS   Petiau G., Dupis A., 1980. Noise temperature coefficient and long-time stability of electrodes for telluric observations, Geophys. Prospect.  28( 05), 792– 804. https://doi.org/10.1111/j.1365-2478.1980.tb01261.x Google Scholar CrossRef Search ADS   Rangel G., 2014. Temperature model and tracer test analysis for the Ribeira Grande geothermal sield, (São Miguel island, Azores), Tech. Rep. UNU-GTP-2014-29, United Nations University, Geothermal Training Programme, (http://www.unugtp.is, last accessed 15 September 2017). Revil A., Cathles III L.M., Losh S., 1998. Electrical conductivity in shaly sands with geophysical applications, J. geophys. Res.  103( B10), 23 925–23 936. https://doi.org/10.1029/98JB02125 Revil A., Hermitte D., Spangenberg E., Cochemé J.J., 2002. Electrical properties of zeolitized volcaniclastic materials, J. geophys. Res.  107( B8), 2168, doi:10.1029/2001JB000599. https://doi.org/10.1029/2001JB000599 Google Scholar CrossRef Search ADS   Revil A., Johnson T.C., Finizola A., 2010. Three-dimensional resistivity tomography of Vulcan's forge, Vulcano Island, southern Italy, Geophys. Res. Lett.  37( 15), L15308, doi:10.1029/2010GL043983. Google Scholar CrossRef Search ADS   Revil A. et al.  , 2011. Hydrogeology of Stromboli volcano, Aeolian Islands (Italy) from the interpretation of resistivity tomograms, self-potential, soil temperature and soil CO2 concentration measurements, Geophys. J. Int.  186( 3), 1078– 1094. https://doi.org/10.1111/j.1365-246X.2011.05112.x Google Scholar CrossRef Search ADS   Sherburn S., Bannister S., Bibby H., 2003. Seismic velocity structure of the central Taupo Volcanic Zone, New Zealand, from local earthquake tomography, J. Volcanol. Geotherm. Res.  122( 1), 69– 88. https://doi.org/10.1016/S0377-0273(02)00470-5 Google Scholar CrossRef Search ADS   Sigmundsson F., Tryggvason E., Alves M.M., Alves J.L., Pálsson K., Ólafsson H., 1995. Slow inflation of the Furnas volcano, São Miguel, Azores, suggested from initial leveling and Global Positioning System measurements, Geophys. Res. Lett.  22( 13), 1681– 1684. https://doi.org/10.1029/95GL01604 Google Scholar CrossRef Search ADS   Silva C., Viveiros F., Ferreira T., Gaspar J.L., Allard P., 2015. Diffuse soil emanations of radon and hazard implications at Furnas Volcano, (São Miguel Island (Azores), in Volcanic Geology of São Miguel Island (Azores Archipelago)  chap. 15, pp. 5– 13, eds Gaspar J.L., Guest J.E., Duncan A.M., Barriga F., Chester D.K., Geological Society of London. Google Scholar CrossRef Search ADS   Simpson F., Bahr K., 2005. Practical Magnetotellurics  Cambridge University Press, Cambridge, UK. Google Scholar CrossRef Search ADS   Tikhonov A., 1950. The determination of the electrical properties of deep layers of the earth's crust, Dokl. Acad. Nauk. SSR  73 295– 297. Ussher G., Harvey C., Johnstone R., Anderson E., 2000. Understanding the resistivities observed in geothermal systems, in Proceedings World Geothermal Congress  Kyushu-Tohoku, Japan, 2000 May 28–June 10. Vanorio T., Virieux J., Capuano P., Russo G., 2005. Three-dimensional seismic tomography from P wave and S wave microearthquake travel times and rock physics characterization of the Campi Flegrei Caldera, J. geophys. Res.  110 B3, B03201, doi:10.1029/2004JB003102. https://doi.org/10.1029/2004JB003102 Google Scholar CrossRef Search ADS   Viveiros F., Cardellini C., Ferreira T., Caliro S., Chiodini G., Silva C., 2010. Soil CO2 emissions at Furnas volcano, São Miguel Island, Azores archipelago: Volcano monitoring perspectives, geomorphologic studies, and land use planning application, J. geophys. Res.  115( B12), B12208, doi:10.1029/2010JB007555. https://doi.org/10.1029/2010JB007555 Google Scholar CrossRef Search ADS   Vozoff K., 1972. The magnetotelluric method in the exploration of sedimentary basins, Geophysics  37 98– 141. https://doi.org/10.1190/1.1440255 Google Scholar CrossRef Search ADS   Wallenstein B., Duncan A., Chester D.K., Rui Marques D., 2007. Fogo volcano (São Miguel, Azores): a hazardous edifice, Géomorphologie : Relief, Processus, Environnement  13( 3), 259– 270. https://doi.org/10.4000/geomorphologie.2853 Google Scholar CrossRef Search ADS   Waxman M.H., Smits L.J.M., 1968. Electrical conductivities in oil-bearing shaly sands, Soc. Pet. Eng. J.  8 107– 122. https://doi.org/10.2118/1863-A Google Scholar CrossRef Search ADS   Weaver J.T., Agarwal A.K., Lilley F.E.M., 2000. Characterization of the magnetotelluric tensor in terms of its invariants, Geophys. J. Int.  141 321– 321. https://doi.org/10.1046/j.1365-246x.2000.00089.x Google Scholar CrossRef Search ADS   Weisenberger T.B., Ingimarsson H., Eyjólfsdóttir E.I., Lévy L., Hersir G.P., Ólafur G.Flóvenz, 2016. Validation of the influence of cation-exchange capacity on resistivity logs, in Proceedings European Geothermal Congress, Strasbourg, France, 19-24 Sept 2016 , (https://www.geothermal-energy.org/publications_and_services/conference_paper_database, last accessed 15 September 2017). Woitischek J., Dietzel M., Inguaggiato C., Böttcher M.E., Leis A., Cruz J.V., Gehre M., 2017. Characterisation and origin of hydrothermal waters at São Miguel (Azores) inferred by chemical and isotopic composition, J. Volcanol. Geotherm. Res.  346 104– 117. https://doi.org/10.1016/j.jvolgeores.2017.03.020 Google Scholar CrossRef Search ADS   Wright P.M., Ward S.H., Ross H.P., West R.C., 1985. State of the art geophysical exploration for geothermal resources, Geophysics  50( 12), 2666– 2696. https://doi.org/10.1190/1.1441889 Google Scholar CrossRef Search ADS   Yardley B.W.D., Bodnar R.J., 2014. Fluids in the continental crust, Geochem. Perspect.  3( 1), 1– 123. https://doi.org/10.7185/geochempersp.3.1 Google Scholar CrossRef Search ADS   Zandomeneghi D., Almendros J., Ibáñez J.M., Saccorotti G., 2008. Seismic tomography of Central São Miguel, Azores, Phys. Earth planet. Inter.  167( 1–2), 8– 18. https://doi.org/10.1016/j.pepi.2008.02.005 Google Scholar CrossRef Search ADS   SUPPORTING INFORMATION Supplementary data are available at GJI online. Figure S1. Apparent resistivity and phase responses plotted for all four components of the impedance tensor. While the solid circles represent the observed data, the solid lines represent the 3-D model responses. Figure S2. Apparent resistivity and phase responses plotted for all four components of the impedance tensor. While the solid circles represent the observed data, the solid lines represent the 3-D model responses. Figure S3. Apparent resistivity and phase responses plotted for all four components of the impedance tensor. While the solid circles represent the observed data, the solid lines represent the 3-D model responses. Figure S4. Apparent resistivity and phase responses plotted for all four components of the impedance tensor. While the solid circles represent the observed data, the solid lines represent the 3-D model responses. Figure S5. Maps of data fit (nRMS) for the high-frequency bands 10,000 1,000 Hz and 1,000 100 Hz. Figure S6. Maps of data fit (nRMS) for the low-frequency bands 100 10 Hz and 10 1 Hz. Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper. APPENDIX: CALCULATION OF SMECTITE MASS FRACTION Here, we summarize the parameter settings used for the calculation of smectite content following the studies by Waxman & Smits (1968), Revil et al. (1998, 2002) and Flóvenz et al. (2005). In Fig. 10, we took a broad interval of fluid conductivity, temperature, and porosity based on the borehole data and modelling results from São Miguel Island (Carvalho et al. 2006), the borehole data obtained by Flóvenz et al. (2005) from a geothermal system in Iceland that we consider as a comparable system to Furnas Volcano, and the laboratory results of Revil et al (2002) on altered samples (see Table A1). The purpose of this study is to demonstrate that the fluid conductivity alone cannot explain the observed value of bulk resistivity of the C1 body; a certain amount of clays characterized by high surface conductivity has to be assumed to explain our inversion results. We suggest the presence of smectite as it has the highest surface conductivity: the mass fraction of about 20 per cent of smectite is enough to explain the observed bulk conductivity of the conductive zone C1 and would still be coherent with high permeability suggested in the fault zone (FF2) from geochemical observations. An additional evidence is given by borehole data in Ribeira Grande (Franco 2015) where the smectite-layer was observed at a depth that is comparable to the depth of the C1 conductive zone. The smectite can be used as indicators of temperature; they are found at a certain temperature interval only and are gradually replaced by less conductive chlorites at a temperature around 200–230 °C (e.g. Franco 2015). The transition from the low to high resistivity below the C1 body (the R2 resistive body or alternatively, the gradual increase in resistivity in the C2 body) could signify this temperature transition. Table A1. Physical properties used for calculation of smectite mass fraction (Fig. 10). Density of rock (borehole data from Flóvenz et al. 2005)  ρ  (kg m−3)  Eq. (A3)  Fluid conductivity  σf  (S m−1)  0.02–3  Hittorf number for NaCl brines (Revil et al. 2002)  t(+)    0.38  Cation exchange capability (Revil et al. 2002)  CEC  (C kg−1)  87 000  Porosity ((borehole data from Flóvenz et al. 2005)  ϕ    0.02–0.25  Cementation exponent (Flóvenz et al. 2005)  m    2.75  Tortuosity (Flóvenz et al. 2005)  a    0.7  Temperature range    (°C)  80–200  Temperature coefficient (eq. A4)  υf    0.02  Temperature coefficient (eq. A4)  υs    0.04  Apparent mobility of counterions at 25 °C (Revil et al. 2002)  β(s)  (m2 s−1 V−1)  5.2 × 10−9  Density of rock (borehole data from Flóvenz et al. 2005)  ρ  (kg m−3)  Eq. (A3)  Fluid conductivity  σf  (S m−1)  0.02–3  Hittorf number for NaCl brines (Revil et al. 2002)  t(+)    0.38  Cation exchange capability (Revil et al. 2002)  CEC  (C kg−1)  87 000  Porosity ((borehole data from Flóvenz et al. 2005)  ϕ    0.02–0.25  Cementation exponent (Flóvenz et al. 2005)  m    2.75  Tortuosity (Flóvenz et al. 2005)  a    0.7  Temperature range    (°C)  80–200  Temperature coefficient (eq. A4)  υf    0.02  Temperature coefficient (eq. A4)  υs    0.04  Apparent mobility of counterions at 25 °C (Revil et al. 2002)  β(s)  (m2 s−1 V−1)  5.2 × 10−9  View Large A number of parameters characterizing the pore fluid or the rock are used in the equations listed below: the fluid conductivity σf, the surface conductivity σs; their ratio is a so-called Dukhin number ζ = σs/σf. The t(+) is the Hittorf number, the fraction of electrical current transported by cations in the electrolyte (see Table A1 for values). The F = aϕ−m is a modified intrinsic formation factor (Archie 1942), depends on the porosity of rock ϕ, tortuosity a, and the cementation exponent m. At low salinities (ζ ≥ t(+)), the bulk conductivity σ is calculated with the Waxman–Smits model (Waxman & Smits 1968; Revil et al. 2002) as a combination of the surface and fluid conductivity contributions   \begin{equation} \sigma = \sigma _s + \frac{\sigma _f}{F}\left( 1 - t_{(+)}\right). \end{equation} (A1) The surface conductivity σs is a function of the cation exchange capacity CEC,   \begin{equation} \sigma _s = \frac{2}{3} \rho \beta _s {\rm CEC} \end{equation} (A2)where the ρ is the density, while βs is the surface mobility of the counterions (Table A1). Smectite bulk density can be calculated following the regression determined in the work of Flóvenz et al. (2005) as   \begin{equation} {\rho _{\rm{smectite}}} = 2.95 - 0.032\phi. \end{equation} (A3) The dependence of σf and σs on temperature can be approximated linearly by   \begin{equation} \sigma _{f,s}(T) = \sigma _{f,s}(T_0)\left[ 1 + \upsilon _{f,s}\left( T - T_0\right) \right], \end{equation} (A4)where T is the temperature, T0 is its reference value (20 °C), and υf, s are the regression coefficients for fluid and surface conductivities, respectively. Following Revil et al. (1998), at higher salinities the bulk electrical conductivity can be expressed as the sum of conductivity contributions by cations (σ(+)) and anions (σ(−)) as   \begin{equation} \sigma = \sigma _{(+)} + \sigma _{(-)}. \end{equation} (A5)These contributions are given by   \begin{eqnarray} \sigma _{(+)} &=& \frac{\sigma _f}{F} \left[ F \zeta + \frac{1}{2} ( t_{(+)} - \zeta )\right. \nonumber\\ &&\left.\times \left( 1 - \frac{\zeta }{t_{(+)} } + \sqrt{\left( 1 - \frac{\zeta }{t_{(+)}} \right) ^2 + \frac{4F}{t_{(+)}} \zeta } \right) \right], \end{eqnarray} (A6)and   \begin{equation} \sigma _{(-)} = \frac{\sigma _{f}}{F} \left( 1 - t_{(+)}\right). \end{equation} (A7) © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

Journal

Geophysical Journal InternationalOxford University Press

Published: Apr 1, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off