One of the difficulties in trying to infer reservoir seismic properties from the laboratory results is that separate phenomena are likely convolved. Variations in fluid temperature change its contribution to the P-wave velocity. At the same time, the properties of the rock matrix vary with temperature as a function of changes in the equilibrium assemblage of hydrothermal minerals. Experiments of Jaya *et al.* [26] take place on much too short a timescale for equilibrium alteration to occur and the experiments of Boitnott [11] which investigate the effect of variations in porosity and variations in the concentration of hydrothermal clays all take place at a single temperature.

Experiments conducted by Kiddle [22] which investigate the effect of temperature and pressure on rock from Montserrat were conducted on dry rock samples from the surface and shallow subsurface and are not representative of saturated altered reservoir rocks. Experiments on saturated samples were also conducted but these were performed at room temperature and pressure. These experiments showed little effect of temperature on P-wave velocity but a large velocity variation was measured between andesite and andesite breccia with P-wave velocities for breccia being significantly lower. The high velocity zones seen in the velocity perturbation maps in the appendix is clearly related to the dense andesite cores associated with the volcanic centers [13,22]. However, there is a significant variation in P-wave velocity in the flank deposits. We suggest that the flank deposits with the anomalously low velocities, particularly the anomaly in the southwest of the island which is co-located with a verified geothermal system, are associated with hydrothermal alteration.

The lateral variation in the P-wave velocity of the flank deposits coupled with the fact the seismic velocity anomalies peak at a depth of 1750 m or so (within the zone of potential phyllic alteration within the reservoir) rather than being anomalous from the surface to depth supports the hypothesis that the velocity variations are dominated by changes in the hydrothermal mineral assemblage.

#### *3.2. Correlation between Seismic Anomaly and Temperature*

Figure 7 suggests, as was postulated in Ryan *et al.* [9], that there is a robust correlation between reservoir temperature and seismic velocity. To test this hypothesis, we recast the data to explore the correlation between the two. For this analysis, we combine the data from wells MON-1 and MON-2 to obtain the maximum number of data points under the assumption that the wells, which are approximately 500 m apart, sample similar regions of the reservoir. This assumption is supported by the similarity between the temperature profiles and the normalized seismic velocity anomalies for both wells. The value of the percentage P-wave velocity perturbation is multiplied by í1 to create positive values for the seismic velocity anomaly.

**Figure 7.** Correlation plot for normalised seismic velocity anomaly (*Vpn*) and reservoir temperature. Data for wells MON-1 and MON-2 are combined. Three different correlation regimes (R1a, R1b and R2) are suggested by the data. One point has been rejected as an outlier. This point which relates to the bottom of well MON-1 is thought to be slightly distorted.

To produce this correlation plot, the resolution of the temperature data was reduced to match the resolution of the resampled seismic velocity model. The temperatures were averaged over 250 m intervals. The data suggest three dominant regimes (*Vpn* refers to normalised seismic velocity anomaly and *T* refers to reservoir temperature in °C):

(1) A slow increase in seismic velocity anomaly with temperature between 29 and 192 °C characterized by the following Equation (2) which fits the data with an *R*2 value of 0.95.

$$V\_{pn} = 0.003T - 0.08\tag{2}$$

(2) A rapid increase in seismic velocity anomaly with temperature between 192 and 219 °C characterized by Equation (3) which fits the data with an *R*2 value of 0.90.

$$V\_{pn} = 0.021T - 3.6\tag{3}$$

(3) A rapid decrease in seismic velocity anomaly with temperature between 219 and 263 °C characterized by Equation (4) which fits the data with an *R*2 value of 0.96.

$$V\_{pn} = -0.015T + 4.4\tag{4}$$

The decision to parameterize the correlation data in the way described above is purely empirical and is strictly only valid for the local area and for the specific set of tomographic inversion parameters used. The linear fits were chosen for their simplicity and because they provide a good fit to the data. The correlation could have been parameterized using different functions.

The observed correlation is between temperature data obtained from well logs and the magnitude of *Vpn* obtained from inversion of active seismic data [13]. As with all geophysical inversions, the inverse model here is non-unique [32] as discussed in Section 2.2. The exact form of the inverse model is dependent on the damping parameters used to stabilize the inversion. The values were chosen to maximize the structural information contained within the seismic data whilst minimizing the effects of errors inherent in the data which make it difficult to obtain a solution to the ill-conditioned inverse problem [32].

The correlations between *Vpn* and temperature we obtain are not solely a function of subsurface geology but also depend on the inversion parameters chosen. A related consequence of the inversion process is the tendency for geological anomalies to be "smeared out" which limits the ability of the inversion to resolve small anomalies. Checkerboard test results shown in the appendix show that the resolving power of the data remains good down to a depth of 3 km. Structures at the sub-kilometer scale are resolvable at this depth in the region of the geothermal reservoir. The checkerboard test data (appendix) also illustrates how the initial model is "distorted" by the inversion process. The intensity of the anomaly decreases and the anomaly spreads out in space. Whilst the intensity of the anomaly varies, its general shape, although blurred tends to be stable. For this reason, we use the normalised seismic velocity anomaly in our analysis to reduce the effect of the particular choice of damping parameters on the magnitude of velocity anomaly.

Our favoured interpretation of the variation in seismic velocity with temperature is that it is related to the assemblage of hydrothermal minerals at different temperatures with clays being formed between 70 and 220 °C which progressively plasticise the rock matrix and reduce bulk and shear moduli of the rock and thus reduce P-wave velocity. Above 220 °C propyllitic alteration occurs which increases the proportion of framework silicate minerals such as quartz, epidote and adularia. This change increases the rigidity of the rock and increases the elastic moduli, again increasing P-wave velocity.

Alternate explanations for the observed low velocity anomalies such as variations in lithology, porosity, fluid saturation and cementation are also possible. None of these alternatives, however, directly explain the variation of the anomaly with depth, in particular why it peaks at a depth of around 1750 m in the region of the two wells. Alternate explanations cannot be totally rejected without further investigation, particularly of the cuttings from the existing wells. Core from future wells would be particularly valuable in future studies.

#### *3.3. Creating a 3-D Temperature Model*

Having determined that a correlation between seismic velocity anomaly and temperature exists, the next step is to use this correlation to estimate the subsurface temperature field from the seismic tomographic data. Knowing how subsurface temperatures vary in 3D will shed light on subsurface fluid flow pathways that would be of interest in geothermal energy exploitation.

Assuming our hypothesis explaining the cause of seismic velocity variation is sufficiently correct, it would be straightforward to convert the seismic anomaly data to temperature estimates if they were uniquely known. As discussed in Section 2.2, inversion of the seismic arrival data leads to "distortion" of the measured anomaly leading to variations in the shape and intensity of the modeled anomaly. If we consider the seismic anomaly data as being made up of a series of 1-D vertical columns through the model space the "distortion" will vary as we move from the edge of the anomaly to the center. In particular, the intensity of the anomaly will appear to decrease as it is "smeared out" at the edges. However, using the assumption that the maxima of the anomaly should be relatively undisturbed even if its absolute value varies and that the maximum is related to a unique process *i.e.*, a change in hydrothermal mineral assemblage which occurs at a fixed temperature in this geothermal reservoir we propose using the normalised seismic velocity anomaly (*Vpn*) to "undistort" the seismic data in the region close to the wells.

The seismic velocity anomalies as measured at the locations of wells MON-1 and MON-2 which are 500 m apart have very similar, somewhat Gaussian, shapes and similar amplitudes (Figure 3). If we look at 1-D anomalies as we move away from the well locations to the edges of the anomaly we see that the shape of the anomaly remains similar but the amplitude decreases (Figure 8). We suggest compensating for this amplitude reduction in the region close to the wells, at which reservoir conditions are similar, by using *Vpn*. We deem this a logical approach because conditions are similar in the region close to the measurement point, and the seismic velocity anomaly has a similar cause and shape. We postulate that the anomaly maximum occurs at a fixed temperature and is caused by a change in the regime of hydrothermal alteration. We therefore suggest that the maximum occurs at the same temperature in locations close to the well. Although this is not a fully rigorous approach, it has the benefit of simplicity and transparency. Our assumption, that nearby regions are similar, should be most robust closest to the zone where temperatures were measured and becomes less certain with distance from that location.

**Figure 8.** Variation in seismic velocity anomaly as one moves from south to north near the location of well MON-1.The graphs relate to vertical columns through the seismic tomography perturbation model going from 5000 to 9500 m·N at 6000 m·E. Each block in the resampled model domain is a cube with side 250 m. The horizontal axis indicates the magnitude of the anomaly in %. Note also how the magnitude of the anomaly increases from the southern edge to a maximum and then decreases again.

Figure 8 shows the variation in the seismic anomaly from the tomographic model going from south to north near well MON-1. This line correlates most closely to the N-S cross-section at 6250 m·E through the estimated temperature model in Figure 9. Figure 3 shows the bell-shaped anomaly seen in the tomographic model in the vicinity of the wells. The peak of this curve is associated in our interpretative framework with a change in hydrothermal mineral assemblage above 220 °C which causes an increase in the rigidity of the rock. An interesting feature of Figure 8 is that we can see the peak of the curve indicating the change in hydrothermal mineral assemblage shallowing from south to north from about 2000 m depth to just under 1500 m. This shallowing is reflected in the estimated temperature model of Figure 9 at 6250 m·E. Figure 9 shows a shallowing of the high temperature region of the model at around 7500 m·N.

To complete this analysis efficiently we created a MATLABTM (Mathworks, Natick, MA, USA) script to automatically analyse the seismic anomaly model. The steps carried out by the script are outlined below:

(4) For each column, the anomaly is analysed to determine its similarity to the master anomaly's roughly Gaussian shape. This is accomplished using a weighted least squares algorithm to fit an offset Gaussian function to the anomaly (Appendix).

Figures 9, 11 and 12 show the results of the temperature estimation algorithm. Although the results of the estimation are, according to our theoretical framework, most justifiable at the areas closest to the wells the algorithm was allowed to operate over the entire seismic anomaly model domain. Results in regions far from the wells must be looked at skeptically.

Looking at the seismic velocity perturbation model (Figure A1) we can see that temperature estimates have been made in the areas corresponding to the three regions in which there are slow seismic anomalies in the northwest, northeast and in the region of interest which is known to host a geothermal system, the southwest. The regions to the northwest and southwest cannot be interpreted with any degree of confidence since they are far from the wells that provide the correlating data. However, the analysis shows that the seismic anomalies in these three regions are similar to each other in the details of their structure.

The method described in this study presents the possibility of a method for estimating subsurface temperatures in a geothermal system with unprecedented resolution. Such a method would be of great use in geothermal exploration and well targeting. The interpretation framework described here requires further work to verify its validity but if it proves to be valid it would deliver another interpretation framework similar to the interpretation scheme used in interpreting resistivity anomalies in high temperature reservoirs [33]. The areas of the model domain where no temperature estimates were made do not necessarily indicate low temperature. It simply means that the form of the seismic anomaly in these regions was dissimilar to that of the master anomaly and no interpretation was made.

**Figure 9.** (**a**–**f**) N-S Cross-sections through the southwestern region of the temperature model estimated from the seismic velocity anomaly data. This region corresponds to seismic anomaly associated with the known geothermal system. Black line indicates the track of well MON-1 and the red line indicates the track of well MON-2. SHV (Soufrière Hills Volcano); SGH (St. George's Hill). The colour scale indicates temperature in degrees Celsius.

**Figure 10.** Fit of seismic anomaly data to an offset Gaussian function. The equation constants are ǻ = í0.57, *k* = 8.76, ȝ =10.26, ı = 3.95, the rms error is 0.80 (see appendix for details). A perfect fit is not necessary. The purpose of the fitting is simply to gauge general similarity of shape and to estimate the amplitude of the anomaly. This anomaly was at location 4500 m·E and 8750 m·N in the model.

Focusing on the southwest anomaly we see that the estimation algorithm makes testable predictions about the geothermal system. The model indicates that the regions of shallowest high temperature occur beneath St. George's Hill. This is consistent with an upflow in this area. The zone of high shallow temperatures can be seen particularly well in Figure 11 at the 1500 m depth. The shallowing of the high temperatures beneath St. George's Hill is seen clearly in Figure 9 in the 6250 m east cross-section.

This finding which indicates the existence of a geothermal system with an upflow beneath St. George's Hill away from the Soufrière Hills volcano is supported by earthquake data and resistivity data described in Ryan *et al.* [9] and is consistent with early results from the well logging.

The temperature anomalies indicated to the NW and NE and seen in Figure 11h are not corroborated by any evidence of active geothermal systems in these areas. However, no extensive studies have been carried out in these areas to date. A possible explanation for these zones is that they are related to zones of relict hydrothermal alteration which has modified the seismic velocities but is no longer related to a real temperature anomaly. This illustrates the constant necessity for multiple sources of information when interpreting geophysical data.

**Figure 11.** (**a**–**g**) Planar views through the SW region of the temperature model estimated from the seismic velocity anomaly data. This region corresponds to the seismic anomaly associated with the known geothermal system; (**h**) Planar view through the entire temperature model at a depth of 2000 m. Apparent temperature anomalies to the NW and NE may only be indicative of relict hydrothermal alteration in these locations or some other process altogether. The NW trending black line indicates a zone of structural weakness identified by Wadge and Isaacs [34]. The NE trending black line indicates a fault plane identified from earthquake hypocenters beneath St. George's Hill [9]. The magenta dots relate to the locations of wells MON-1 and MON-2. The cyan dot relates to the location of the alkali-chloride hot spring at Hot Water pond. GBH (Garibaldi Hill); SGH (St. George's Hill); GM (Gage's Mountain); CP (Chances Peak); SHV (Soufrière Hills Volcano); RM (Roche's Mountain); SSH (South Soufrière Hills); RE (Roche's Estate); CH (center Hills); SH (Silver Hills). The colour scale indicates temperature in degrees Celsius.

**Figure 12.** (**a**–**f**) E-W Cross-sections through the southwestern portion of the temperature model estimated from the seismic velocity anomaly data. This region corresponds to seismic anomaly associated with the known geothermal system. Black line indicates the track of well MON-1 and the red line indicates the track of well MON-2. SGH (St. George's Hill). The colour scale indicates temperature in degrees Celsius.

#### **4. Conclusions and Further Work**

The main conclusion of this study is that there is persuasive evidence to suggest that seismic velocities in a high temperature geothermal reservoir are strongly influenced by temperature dependent hydrothermal mineral assemblages. We suggest that this influence may allow seismic velocity anomalies to be used to estimate reservoir temperatures. Further work, however, is required to fully reject alternate possibilities.

Sub-surface temperature estimation suggests the existence of an upflow zone beneath St. George's Hill in the geothermal system in the southwest of Montserrat. Estimated "temperature anomalies" in the northwest and northeast of the island are possibly due to relict hydrothermal alteration related to extinct geothermal systems.

The current interpretation scheme for the seismic velocity anomaly data is based on a pattern of hydrothermal alteration commonly seen in high temperature geothermal systems [15]. However, further work on the petrology of samples from Montserrat's geothermal wells is required to determine whether the details of hydrothermal alteration in this case are consistent with the seismic interpretation.

The methods used to correlate temperature log data to the seismic anomaly data, and particularly to "undistort" the seismic anomaly data, are empirical and crude. A more rigorous approach should be sought in future work.

In a similar vein, it may be possible to obtain a more theoretically robust method for correlating seismic anomalies to hydrothermal alteration by using Gassman's equation [27] in tandem with experimental results on the physical properties of hydrothermally altered rocks.

#### **Acknowledgments**

We thank the government of Montserrat for permission to publish the results of this study. We are also very grateful to those who supported the SEA-CALIPSO study including NSF, NERC, Discovery Channel, British Geological Survey, Foreign and Commonwealth office, PASSCAL Instrument centre and the Montserrat Volcano Observatory. Although originally for another purpose this data set has proved to be a unique and valuable addition to research in geothermal exploration; the kind of happy accident which often occurs in science. We are also grateful to Pat Browne and Bridget Lynne for many informative discussions on hydrothermal alteration in geothermal systems. The Institute of Earth Science and Engineering is also acknowledged for providing us with time to work on this study. We gratefully acknowledge the three anonymous reviewers whose thoughtful work helped to greatly improve the clarity of this manuscript.

#### **Author Contributions**

Graham Ryan formulated the project's original concept and wrote the manuscript. He is also responsible for developing the algorithms used to analyse the well temperature log and tomographic model data, and wrote the scripts used to generate the 3-D temperature model. Eylon Shalev was responsible for creating the *Vp* tomographic model along with data for the *Vp* perturbation model, ray hit and checkerboard test plots. He also reviewed the manuscript.

#### **Appendix**

#### *A.1. Seismic Tomography*

**Figure A1.** P-wave tomography results from Shalev *et al.* [13] displayed as perturbation from the average velocity at each depth. Blue represents faster velocities and red represents slower velocities. (**a**) map view slice through the target volume at 2000 m depth; (**b**) E-W cross-sections through target volume at 8750 m north; and (**c**) 6750 m north; (**d**) N-S cross-section through target volume at 5250 m east; and (**e**) 6250 m east. The NW trending black line indicates a zone of structural weakness identified by Wadge and Isaacs [34]. The NE trending black line indicates a fault plane identified from earthquake hypocenters beneath St. George's Hill [9]. The magenta dots relate to the locations of wells MON-1 and MON-2. The cyan dot relates to the location of the alkalichloride hot spring at Hot Water pond. GBH (Garibaldi Hill); SGH (St. George's Hill); GM (Gage's Mountain); CP (Chances Peak); SHV (Soufrière Hills Volcano); RM (Roche's Mountain); SSH (South Soufrière Hills); RE (Roche's Estate); CH (center Hills); SH (Silver Hills). The colour scale indicates % velocity perturbation.

**Figure A1.** *Cont.*

**Figure A2.** Number of observations used in the seismic tomography inversion shown for each on-shore seismometer. Red dots indicate Reftek seismometers. Green dots indicate Texan miniature seismic recorders and blue dots indicate seismometers which are part of the permanent monitoring network of the Montserrat Volcano Observatory.

**Figure A3.** Location of seismometers (red triangles) used in this study located on-island and off-shore. The green line indicates the track of the vessel towing the air-gun seismic source. The green and purple dots mark the positions of wells MON-1 and MON-2, respectively. The blue outline surrounding the island indicates the extent of a shallow submarine platform.

**Figure A4.** Ray hit counts for the tomographic inversion domain including. The panels show the number of rays that pass through 250 m × 250 m × 250 m boxes at (**a**) 0 m; (**b**) 500 m; (**c**) 1000 m; (**d**) 2000 m and (**e**) 3000 m depth. Only the center portion of the tomography box is plotted, with the outline of Montserrat in black. Colors indicate the log10 of the number of rays that track through each box.

**Figure A5.** Checkerboard test results for anomalies at (**b**) 0 km; (**c**) 0.5 km; (**d**) 1 km; (**e**) 2 km; (**f**) 3 km and (**g**) 4 km depths; (**a**) is the source pattern for the checkerboard test (shown at 3 km depth). The *x* and *y*-axes indicate distance in km east and north respectively. The colour scale indicates percentage velocity perturbation. Synthetic noise with a zero-mean Gaussian distribution and 0.02 s standard deviation was added to the travel-time data. Note the variation in the amplitudes of anomalies.

**Figure A4.** *Cont.*

#### *A.2. Gaussian Fitting Algorithm*

In order to determine the similarity of 1-D columns through the seismic velocity anomaly model to the seismic anomalies observed at the locations of wells MON-1 and MON-2, we exploited the similarities of the seismic anomalies as a function of depth to offset Gaussian functions following the equation:

$$W\_{\rho} = k \exp \left( \left( d - \mu \right)^{2} \Big/ \int\_{\sigma^{2}} \right) + \Delta \tag{A1}$$

where *Vp* is the magnitude of the seismic anomaly, *k* is a scaling factor, *d* is depth below sea level in meters, ȝ is the depth at which the peak *Vp* magnitude occurs, ı is the characteristic width of the Gaussian function and ǻ is a scalar offset.

In order to obtain a simple linear function we assume that ǻ is a known value and simply minimize the remaining Gaussian function once it has been subtracted:

$$V\_{\rho} = k \exp{-\left(\left(d - \mu\right)^{2} \bigg\prime\_{\sigma^{2}}\right)}\tag{A2}$$

To obtain a linear equation, we take the natural logarithm of both sides of the equation to obtain a polynomial in *d*:

$$
\ln V\_{\rho} = \left( \begin{matrix} -\frac{1}{\rho^2} \\ \sigma^2 \end{matrix} \right) d^2 + \left( \begin{matrix} 2\mu \\ \sigma^2 \end{matrix} \right) d + \left( \ln k - \left\lfloor \mu^2 \Big/ \sigma^2 \right\rfloor \right) \tag{A3}
$$

Treating the variables *d* and *d2* as linear variables, we can do a simple weighted least squares parameter estimation using Equation (A4).

In order to bias the fitting algorithm to fit the peaks of the anomalies at the expense of the tails, we employed a weighted least squares algorithm. In which the errors varied monotonically from 1/10th to 1/100th of the average value of ln*Vp* from the initial low value to the peak value of the anomaly and then back to 1/10th. This error information is captured in the diagonal matrix C where the leading diagonal contains the squared reciprocals of the corresponding error values:

$$
\left\langle \mathbf{x}\_{\mathbf{w}} \right\rangle = [\mathbf{A}^T \mathbf{C} \mathbf{A}] \mathbf{A}^T \mathbf{C} \mathbf{b} \tag{A4}
$$

where is the vector of linear best-fit coefficients of Equation (A3), A is the matrix of linear variables in Equation (A3) and b is the vector of ln*Vp* values. The parameters ȝ, ı and *k* can be recovered from simple rearrangement of the elements of . To generate the best-fit parameters from data from a 1-D seismic anomaly, we simply use Equation (A4) using the values of ln *Vp*, *d* and *d2*. A description of weighted least squares estimation can be found in any book on elementary linear algebra, e.g., [35]. *wx wx*

The RMS error of the fit of the data to the Gaussian model is obtained using the Equation (A5):

$$RMSError = \left[\frac{\sum\_{1}^{n} \left(V\_{\text{pmax}} - V\_{\text{pft}}\right)^{2}}{n}\right]^{\frac{1}{2}}\tag{A5}$$

where *Vpmeas* is the value from the seismic velocity anomaly model and *Vpfit* is the estimate derived from the best fit model. *n* is the number of data points in the 1-D column.

#### **Conflicts of Interest**

The authors declare no conflict of interest.

#### **References**

## **Determining the Maximum Depth of Hydrothermal Circulation Using Geothermal Mapping and Seismicity to Delineate the Depth to Brittle-Plastic Transition in Northern Honshu, Japan**

#### **Yota Suzuki, Seiichiro Ioka and Hirofumi Muraoka**

**Abstract:** This paper defines the maximum possible vertical extent of hydrothermal circulation in granitic crust, and thus the maximum depth within which geothermal reservoirs can be encountered. To evaluate prospective geothermal fields we constructed a geothermal database in northern Honshu, Japan that includes 571 points of thermal data of existing wells and hot springs. Depth-temperature curves were normalized by the Activity Index for three-dimensional extrapolation and a depth contour map of the 380 °C isotherm was plotted as an assumed brittle-plastic transition for granitic crust. Shallower-depth anomalies of the brittle-plastic transition on this map are closely coincident with the Quaternary volcanoes and their prospective geothermal fields. It should be noted that the bottom of the spatial distribution of seismicity in the volcanic fields shows strong correlation to the 380 °C isotherm. This result indicates reliability of the subsurface three-dimensional thermal map and suggests that the 380 °C isotherm strongly constrains the bottom surface of seismicity, fracturing and hydrothermal convection in granitic crust.

Reprinted from *Energies*. Cite as: Suzuki, Y.; Ioka, S.; Muraoka, H. Determining the Maximum Depth of Hydrothermal Circulation Using Geothermal Mapping and Seismicity to Delineate the Depth to Brittle-Plastic Transition in Northern Honshu, Japan. *Energies* **2014**, *7*, 3503-3511.

#### **1. Introduction**

Since the unprecedented disastrous Great East Japan Earthquake at 11 March 2011, geothermal energy is being re-evaluated as a more crucial energy resource, not only as green energy, but also as indigenous energy in northern Japan. When compared to other renewable energy sources, the main strength of geothermal energy among the variety of renewable energy sources is that it can be produced stably in all weather conditions, whereas its main weakness is that it uses invisible underground resources prolonging development lead time as well as increasing the initial investment risks and costs. In order to explore subsurface geothermal resources, geological, geochemical and geophysical exploration techniques have been developed, but one of the efforts that should be done is the estimation of subsurface thermal structure from thermal data such as existing wells and hot springs.

Fortunately 27,219 hot spring sources are distributed all over Japan [1], and numerous hot springs are concentrated in Aomori Prefecture, a northern end of Honshu Island, Japan (Figure 1). Most of these Japanese hot springs are artificially developed by shallow wells. We constructed a geothermal database in Aomori Prefecture that includes available information on geothermal exploratory wells, hot spring wells, and natural hot springs.

**Figure 1.** Map showing the location of the Aomori Prefecture, a northern end of Honshu Island, Japan. Red triangles show the main Quaternary volcanoes and calderas.

The neighboring Iwate Prefecture and Akita Prefecture have three installed geothermal power plant units, respectively, including the oldest Matsukawa plant in Japan that has been working for an almost half century since 1966, but unfortunately no geothermal power plants have been developed in Aomori Prefecture so far. The geothermal database we have set up, will assist geothermal development in Aomori Prefecture before detailed geothermal exploration surveys are carried out in individual geothermal fields.

The geothermal database can be applied to draw a variety of subsurface thermal structures. As the spatial density of the point data is still relatively limited, it will be adequate to apply the data to a broad but fundamental structure. One of the objectives would be to map the depth of the brittle-plastic transition because it is synonymous with the bottom of the permeability and comprises the broadest as well as most fundamental structure delineating a bottom of the entire hydrothermal convection sphere.

This paper describes a geothermal database in Aomori Prefecture and draws a map of the depth of the brittle-plastic transition based on some depth-temperature normalization techniques that allow three-dimensional extrapolation in up-flow zones. Then, this map of the depth of the brittle-plastic transition is verified by the bottom of seismicity which is far independent destruction phenomena from the thermal structure concerned. We used the Generic Mapping Tools 4.5.9 for drawing most of figures in this paper [2].

#### **2. Geothermal Database**

Our geothermal database consists of 35 points of geothermal exploratory well data, 510 points of hot spring well data and 26 points of natural hot spring data. Our geothermal database was based on two main data sources, both of which are open to the public. One is "*Record of Hot Spring Geology in Aomori Prefecture*" published by the Department of Health and Welfare, Aomori Prefectural Office in 1997 [3]. The other is "*Atlas of Hydrothermal Systems in Japan*" published from the Geological Survey of Japan, the National Institute of Advanced Industrial Science and Technology (AIST) in 2007 [4].

The "*Record of Hot Spring Geology in Aomori Prefecture*" describes 434 points of hot spring data in Aomori Prefecture by their locations, well depths, discharge temperatures, discharge rates, major chemical constituents and geological columns along the wells. Among them 26 hot springs are natural hot springs with no artificial drilling. We used the 26 points of natural hot spring data and 395 points of hot spring data developed by wells of which both the well depth and discharge temperature are known. We digitized all the paper-based data from the report.

The "*Atlas of Hydrothermal Systems in Japan*" describes 7203 points of hot spring data and 3066 well data points from all over Japan for drawing a variety of geothermal maps. We used 115 points of hot spring data and 35 points of geothermal exploratory well data from Aomori Prefecture. They were also originally paper-based data, but a digital-version on CD-ROM was published later [5].

#### **3. Estimate of Subsurface Geothermal Structure**

Determining accurate subsurface three-dimensional thermal structure is one of the ultimate goals on geothermal exploration. However, depth-temperature curves obtained by temperature logging in geothermal wells show usually complicated patterns due to the permeability variation of given geological formations with depth and three-dimensional geothermal structure that cannot easily be acquired. Then, normalization techniques for depth-temperature curves are necessary. One of the normalization techniques is the Activity Index originally proposed for evaluating temperature ranks of geothermal fields by Hayashi [6].

Geothermal exploration mainly concerns up-flow zones of the hydrothermal convection system rather than down-flow zones. The highest temperature curves in up-flow zones are normally limited by the boiling point temperature curve [7]. The lowest temperature curves are normally limited by the linear thermal conduction curve with an average continental geothermal gradient 30 °C/km. The Activity Index (*AI*) is defined by the following equation:

$$AI = \frac{a}{b} \times 100 = \left(1 - \frac{Tb - Tm}{Tb - Tg}\right) \times 100\tag{l}$$

where *Tm* is the maximum temperature at the observed depth, *Tb* is the boiling point temperature curve of pure water at the observed depth, and *Tg* is an average geothermal gradient (30 °C/km) at the observed depth. In addition, *a* and *b* show the *Tm í Tg* and *Tb í Tg* respectively. The concept of the Activity Index is graphically represented (Figure 2).

**Figure 2.** A concept of the Activity Index. Six solid curves show *AI* = 0, *AI* = 20, *AI* = 40, *AI* = 60, *AI* = 80 and *AI* = 100, respectively. Red circles show 571 points of thermal data in Aomori Prefecture.

The *AI* 100 curve exactly coincides with the boiling point temperature curve of pure water and is here calculated in the approximation equation by Haas [7]. When the Activity Index can be determined at a given depth, we can extrapolate the normalized temperature curve to the arbitrary depth using the Equation (1) [6]. Therefore, we can easily plot the subsurface three-dimensional thermal structure by this normalization technique. When the temperature logging data are available, the borehole temperature is used as *Tm* at the given depth. However, only discharge temperature is available on most of hot spring wells. *Tm* was approximated by the discharge temperature on these hot springs. It is convenient that the discharge temperature in Celsius degree of natural hot springs at the ground surface can be adopted as the same numbers of *AI* as discussed by Hayashi [6]. As natural hot springs represent the surface manifestation of a hydrothermal up-flow, their temperatures tend to form a cluster around 100 °C. However, natural hot springs are only 26 pieces among the entire 536 hot springs in this database and no biases are found in the normalization by the Activity Index. Chemical geothermometry is useful to estimate reservoir temperatures, but we did not use the method because it does not provide the depth information. As a result, we can use 571 point data for *AI* mapping in Aomori Prefecture (Figure 2).

Pure water reaches a critical point temperature near the depth of 3500 m. However, the critical point of water dramatically shifts to the higher temperature with increasing salinity in the natural brine system so that we simply extrapolated the *AI* curves to the depth of 15 km. Then, we can plot the subsurface three-dimensional thermal structure of upper crust in Aomori Prefecture.

We here try to plot the brittle-plastic transition which is critically important to delineate the bottom surface of hydrothermal convection, particularly to extract the shallower apices of the isothermal surface associated with high-level magmatic intrusions. Deformed rocks within fault zones indicate that the transition from purely brittle to purely plastic deformation processes may occur over a relatively broad range of temperature and pressure [8]. This range is often denoted as a semi-brittle region between the brittle-ductile transition (BDT) and the brittle-plastic transition (BPT) [9]. For quartz and feldspar rich rocks in continental crust this range occurs at an approximate temperature in the range from 250 to 400 °C.

The temperature of the *in-situ* brittle-plastic transition in an actual borehole was too high to be measured by the conventional temperature logging tools but it was successfully measured by a variety of techniques in the well WD-1a which was drilled into the depth of 3729 m in the Kakkonda geothermal field, Iwate Prefecture, northern Japan in 1995 [10]. The bottom-hole temperature exceeded 500 °C (Figure 3) [10,11]. The temperature profile includes a temperature inflection point at 380 °C and at the depth of 3100 m where a boiling point curve above the depth shows a hydrothermal convection zone and a steep gradient thermal conduction curve below the depth shows a non-fracture and non-hydrothermal convection zones [10,11]. The temperature inflection point at 380 °C is thus considered to be the brittle-plastic transition (BPT) on granitic crust [12] (Figure 3). On the other hand, the maximum strength zone was estimated at the temperature 340 °C and at the depth of 2000 m by the high density of fractures, measurements of the differential strain curve analysis (DSCA) [13] on oriented cores and the theory of the strength envelope of lithosphere [12]. This maximum strength zone can be denoted as the brittle-ductile transition (BDT; Figure 3). The brittle-ductile transition (BDT) is important for a nest of earthquakes and fractures because of the preferred stress concentration on the maximum strength zone [12]. From a geothermal point of view, the brittle-plastic transition (BPT) is crucial because it is a bottom surface of seismicity, fracturing and hydrothermal convection.

**Figure 3.** Synoptic models of the brittle-plastic transition (BPT) and the brittle-ductile transition (BDT). (**a**) Depth-temperature profile along the well WD-1a at the Kakkonda geothermal field, Iwate, Japan [10,11]. (**b**) Depth-strength envelope along the well WD-1a at the Kakkonda geothermal field, Iwate, Japan [12]. Ȝ shows pore pressure.

Thus, the most reliable temperature of the brittle-plastic transition is 380 °C on granitic crust [10–12], and we tried to draw the depth of the 380 °C isothermal surface from the Activity Index. Figure 4 shows the calculated results of the assumed brittle-plastic transition. Most of the depth of the brittle-plastic transition showed deeper than 12 km below sea level. However, the prominent shallow depth anomalies coincide with the Quaternary clustered volcanoes such as Hiuchi-dake and Osore-zan in the Shimokita Peninsula to the north and Hakkoda-san to the south. The map is consistent with known prospective geothermal fields.

**Figure 4.** Map of the depth of the 380 °C isotherm below sea level as a brittle-plastic transition in Aomori Prefecture, Japan. Solid dots show 571 points of thermal data in Aomori Prefecture.

#### **4. Discussion: Verification by Seismicity**

Shallow-depth anomalies of the 380 °C isotherm are exactly consistent with the areas of the Quaternary volcanoes and their prospective geothermal fields (Figures 1 and 4). The shallowest apex is observed at a depth about 5 km in the Hakkoda-san clustered volcano area that seems not necessarily shallow as a magmatic heat source region. It is due to the regional-scale smooth contouring from the scarce random point data and 0.2 min (about 370 m) gridding. The map would provide the subsurface three-dimensional thermal structure and seems useful for subsurface geothermal assessments. However, further verification would be expected from a scientific point of view. To verify the reliability of the subsurface three-dimensional thermal structure on this map seismicity is useful because the seismic phenomenon is primarily independent from thermal phenomena.

We used the earthquake catalog provided by the Japan Meteorological Agency (JMA) [14]. The JMA earthquake data were collected in Aomori Prefecture during the period from January 1990 to July 2012. As these data include deep subduction zone earthquakes, we selected earthquake data with their epicenters to be shallower than the depth of 20 km. To verify the reliability of the map (Figure 4), the spatial distribution of the earthquakes is compared with the thermal structure (Figure 5). The thermal point data such as wells and hot springs are restricted on the land areas. Therefore, we selected earthquake data on the land area too.

**Figure 5.** Comparison between the 380 °C isotherms and seismicity on cross sections. (**a**) Hakkoda volcano area; the upper shows the east-west cross section and the lower shows the north-south cross section (Figure 4); (**b**) Shimokita Peninsula area; the upper shows the east-west cross section and the lower shows the north-south cross section (Figure 4). Red triangles show the main Quaternary volcanoes and calderas. Solid lines show the 380 °C isotherm and open circles show epicenters of earthquakes enlarging the diameter with magnitude during January 1990 to July 2012.

Figure 5 shows the estimated 380 °C isotherm produced from the Activity Index. Hakkoda volcano area and Shimokita Peninsula area show the prominent anomalies that are characterized by the dramatically shallower depth of the brittle-plastic transition than other areas. The bottom of the

spatial distribution of seismicity shows the strong correlation to the 380 °C isotherm. This result indicates that the subsurface three-dimensional thermal map in Figure 4 is basically reliable. This result also suggests that the 380 °C isotherm strongly correlates to the bottom surface of seismicity, fracturing and hydrothermal convection. Geothermal database and the Activity Index are useful to draw the subsurface three-dimensional thermal structure.

#### **5. Conclusions**

Based on the geothermal database for available thermal data such as wells and hot springs in Aomori Prefecture, northern Honshu, Japan, depth-temperature curves were normalized by the Activity Index for three-dimensional extrapolation and a depth contour map of the 380 °C isotherm was plotted as an assumed brittle-plastic transition for granitic crust. Shallower-depth anomalies of the brittle-plastic transition on this map are closely coincident with the Quaternary volcanoes and their prospective geothermal fields. It should be noted that the bottom of the spatial distribution of seismicity in the volcanic fields shows strong correlation to the 380 °C isotherm. This result indicates reliability of the subsurface three-dimensional thermal map and suggests that the 380 °C isotherm strongly constrains the bottom surface of seismicity, fracturing and hydrothermal convection in granitic crust.

#### **Acknowledgments**

This work was supported by the Japan Society for the Promotion of Science (JSPS) Kakenhi (Grant-in-Aid for Scientific Research B) Grant Number 26281053. We thank the Japan Meteorological Agency (JMA) and the Ministry of Education, Culture, Sports, Science and Technology in Japan (MEXT) for providing the earthquake epicenter data. This paper benefited from the comments and suggestions of three anonymous reviewers.

#### **Author Contributions**

Yota Suzuki is a principal investigator of this work so that he has made the geothermal database in Aomori Prefecture, the first manuscript and all the figures. Seiichiro Ioka advised on hydrology and the data quality of the geothermal database. Hirofumi Muraoka co-operatively provided many papers on the geothermal database and the Kakkonda project.

#### **Conflicts of Interest**

The authors declare no conflict of interest.

#### **References**

1. Ministry of the Environment, Japan. The State of Use of Hot Springs in Japan in the Fiscal Year 2012. Available online: http://www.env.go.jp/nature/onsen/data/riyo\_h24\_1.pdf (accessed on 1 May 2014). (In Japanese)

### **3D Geothermal Modelling of the Mount Amiata Hydrothermal System in Italy**

#### **Paolo Fulignati, Paola Marianelli, Alessandro Sbrana and Valentina Ciani**

**Abstract:** In this paper we build a subsurface model that helps in visualizing and understanding the structural framework, geology and their interactions with the Mt. Amiata geothermal system. Modelling in 3D provides the possibility to interpolate the geometry of structures and is an effective way of understanding geological features. The 3D modelling approach appears to be crucial for further progress in the reconstruction of the assessment of the geothermal model of Mt. Amiata. Furthermore, this model is used as the basis of a 3D numerical thermo-fluid-dynamic model of the existing reservoir(s). The integration between borehole data and numerical modelling results allows reconstructing the temperature distribution in the subsoil of the Mt. Amiata area.

Reprinted from *Energies*. Cite as: Fulignati, P.; Marianelli, P.; Sbrana, A.; Ciani, V. 3D Geothermal Modelling of the Mount Amiata Hydrothermal System in Italy. *Energies* **2014**, *7*, 7434-7453.

#### **1. Introduction**

The Mt. Amiata volcano-geothermal area, located in Southern Tuscany (Figure 1), is characterized by very high heat flux (up to 600 mW/m2 ) and hosts two high-temperature geothermal fields (Bagnore and Piancastagnaio) that, together with Larderello-Travale, are the only economically exploited fields in Italy for power production.

Mt. Amiata is the largest Tuscan volcano, consisting of dacitic, rhyodacitic and minor olivine-latitic lavas that erupted in a period ranging from 300 to 190 ka [1,2]. It is located on a structural high delimitated to the East by the normal (border) faults of the western side of Siena-Radicofani Neogene basin. Mt. Amiata volcanic products are mainly effusive and of silicic composition. Volcanics cover an area of about 90 km2 (Figure 1) lying on "Tuscan" units (metamorphic rocks of Palaeozoic age forming the basement, Triassic-Jurassic carbonate evaporitic successions and Cretaceous-Oligocene terrigenous formations), Oligocene-Mesozoic Ligurian units (shale, sandstones and marles), Palaeogene Subligurian units and Neogene sediments [3–8]. The volcanic activity was triggered by the activation of regional transcurrent faults following the Mio-Pliocene extensional phase that induce an important crustal thinning. This latter and the associated tectonic favoured the emplacement of a magma body at 6–7 km depth [9–11]. The relatively shallow intrusion of low-density magma was at the beginning of the uplift (up to 1000 m) of this area during the Quaternary [6,11]. The main fault system in the Mt. Amiata volcanic complex has an ENE-WSW direction. This direction represents the main fissure system along which magma raised from the shallow magma chamber. Nearly all the lava domes are aligned along this axial fissure system [1] known as the Mt. Amiata fault (see the partly coinciding location of faults and eruptive centres in Figure 1). Other transcurrent fault systems occur North and South of the volcano. The active geothermal fields as well as Hg-(Sb) mineralization of Mt. Amiata area are linked to these fractured systems and to the thermal anomaly induced by the shallow intrusion (Figure 1).

**Figure 1.** Geological sketch map of the area of Mt. Amiata and its surroundings [8]. The whole area coincides with the *X*-*Y* extents of the 3D geological model; the rectangle in red coincides with the *X*-*Y* extents of the numerical model.

Geothermal exploration of the Mt. Amiata area started at the beginning of the 50s, and the first overview of the geothermal system was carried out in 1970 [3]. Previous studies were developed for mining activities linked to cinnabar production [12]. All available data (more than 100 drillings from 500 to about 4000 m depth, as well as many geophysical, geological, structural, geochemical and hydrogeological studies), derived from the exploration and exploitation activities mainly carried out for geothermal purposes, allow developing a well-constrained geothermal model.

The goal of this paper is to build a fit-for-purpose 3D geological model of Mt. Amiata geothermal system. Modelling in 3D provides the possibility to interpolate the geometry of geological structures and is an effective way of understanding the subsoil geology. As a consequence, a 3D modelling approach appears to be crucial to make further progress in the understanding of the geothermal-geological system of Mt. Amiata. Furthermore, this geological model is used as the basis of a 3D numerical thermo-fluid-dynamic simulation of the volcanic geothermal system that led to an integrated geothermal model of the area.

#### **2. Mt. Amiata Geothermal System**

Mt. Amiata is a classic volcanic-intrusive geothermal system. This is characterized by clay-rich sedimentary units (Ligurian alloctonous nappes and Neogene units) that represent a very efficient impermeable cover. Two confined water-dominated geothermal reservoirs occur. The shallower one (500–1500 m of depth) is hosted in carbonate-evaporitic rocks belonging to the Tuscan Nappe, in particular in the "Calcare Cavernoso" formation and locally in fractured carbonate units. Calcare Cavernoso is a vacuolar carbonate breccia characterized by dolostone porphyroclasts (often dissolved) cemented by carbonate cement resulting from alteration of gypsum evaporite layers. This formation is often very permeable due to fracturing and porosity. The mean temperatures found in this reservoir range from 150 to 220 °C. The reservoir hosts a two-phase water dominated aquifer. In correspondence of structural highs (traps) this reservoir presents pressurized gaseous caps in which gas and steam accumulate [13]. The gaseous phase is mainly constituted by CO2, in accordance with the very high CO2 degassing that characterizes the whole Mt. Amiata area [14,15]. Shallow reservoir outcrops in two different areas (North East and South of Mt. Amiata). These probably represent the main recharge areas of the geothermal fields.

The deeper reservoir (below 2500–3000 m of depth) is hosted in thermometamorphic fractured phyllites and metasandstones of the Tuscan Metamorphic Complex, modified by thermometamorphic processes [16]. These belong to the termometamorphic aureole of the recent silicic intrusion. In the deeper reservoir temperatures up to 300–350 °C are actually measured. The fluid in this reservoir is a two-phase too [3,10]. This represents the main reservoir industrially exploited at present [13] with a total power production of 88 MWe.

Below the deeper reservoir, reflection seismic profiles reveal the occurrence of a geopressurized system evidenced by discontinuous "bright spot"-type reflections, named K-horizon [17–21], at depth between 4 and 7 km. This is assumed to be the root of the hydrothermal systems in the Mt. Amiata area. The temperatures estimated for the K-horizon are considered to be above 400 °C [17,22]. This horizon is present also outside the geothermal region but at much greater depth [18]. The interpretation of K horizon is still debated: some authors [18,21] suggested, basing on the earthquake hypocentres distribution, that the K-horizon corresponds to a ductile-brittle boundary; others [22] suggested that this ductile-brittle transition occurs in coincidence with high fluid pressure; others [23] explain the K-horizon as the roof of over-pressurized reservoirs on the basis of the interpretation of deep crustal seismic lines.

The Mt. Amiata intrusion, present at about 6 km of depth, represents the heat source of the hydrothermal system described above.

#### **3. 3D Geological and Geothermal Model**

#### *3.1. Available Data*

Geological models in three dimensions integrate different kinds of data (*i.e.*, geological maps, cross-sections, seismic profiles, gravity profiles, logs of deep wells *etc.*) much more easily than traditional techniques and allow the three-dimensional (3D) structure to be represented more coherently and exhaustively [24–26]. The available data implemented in this work are reported as follow:

The geological maps used in this work derive from the recent 1:10,000 maps of the CARG (CARtografia Geologica) project [27]. The cross sections are taken from CARG project and from literature [3,4,6,8,27–29].

More than 100 deep geothermal boreholes (from 153 to 4836 m of depth) occur (UNMIG, Ufficio Nazionale Minerario per gli Idrocarburi e le Georisorse geothermal database) in the Mt. Amiata area. Geothermal wells (Figures 1 and 2) are mostly located in correspondence with the active geothermal fields of Bagnore and Piancastagnaio. Well data concern stratigraphy, temperatures, pressure of the fluid phase, flow rates, technical drilling profiles, depth of the water table. Unfortunately no geophysical borehole logs are available, with exception of temperature logs, because in some cases they were not carried out and in others the data were not published by the owner companies. Many other shallower wells (from 100 to 750 m of depth) were drilled in the past for mining activities. All data concerning these wells were also used in this work.

Data from geophysical investigation, derived from past geothermal and hydrocarbon (this latter in the Radicofani graben) exploration, are available (Figure 3). The geophysical data used in this work are: the map of the resistive substratum [3] and seismic profiles [4,20,29,30]. Regional gravity data were also incorporated in the study [31–33].

#### *3.2. Data Integration Toward a 3D Geological Model*

All data were implemented in the 3D modelling Petrel 2011 platform (Schlumberger Limited, Houston, TX, USA). Data include both hard copies (digitalized and georeferenced) and GIS (geographical information system) format data. Geological data were elaborated in the geodatabase. Geological units were grouped on the basis of their "geothermal" characteristics, essentially in terms of cover type and reservoir type formations, taking into account their permeability (Figure 4).

Six units were derived as follow: (1) the volcanic complex: fractured and permeable lavas; (2) the Neogene and Quaternary mainly impermeable sedimentary deposits: clay, sands and minor conglomerate; (3) the Ligurian units mainly impermeable: shale, sandstones and marls; (4–5) the Tuscan Nappe, subdivided into two sub-units: (4) an upper sub-unit denoted as TN 2, represented by the mainly impermeable terrigenous successions and (5) a lower sub-unit that is denoted as TN 1 and represented by the fractured and permeable carbonate-evaporitic successions (shallow geothermal reservoir); (6) the Palaeozoic basement, formed by the Tuscan Metamorphic Complex constituted by low permeability phyllites and permeable thermometamorphosed fractured phyllites (this latter is the deep geothermal reservoir). The geological map of the area (Figures 1 and 2) and the geological cross sections (Figure 5) [3,4] were elaborated on the basis of the above reported criteria, digitalized and uploaded in the model.

**Figure 2.** 3D view of geological map spread on DEM (digital elevation model) of the Mt. Amiata area (vertical dimension is exaggerated by factor of 2). *Z*-axis in meters, *X*-axis and *Y*-axis: the coordinate system used is WGS 84/UTM zone 32N (EPSG: 32632). Arrow points to the North. Legend as reported in Figure 1.

**Figure 3.** Map of the Mt. Amiata area with the location of the deep boreholes and the traces of geological, seismic and gravimetric sections.

**Figure 4.** Sketch of the geological units used in this work (right) and the complete geological succession taken from literature. Stratigraphic succession in the Monte Amiata area (modified after [34,35]: see the references for further information). Ligurian and Subligurian Units; Ophiolitic Unit; APA: "Argille a Palombini" Fm.; CCL: "Calcari a calpionelle" Fm.; DSA: "Diaspri" Fm.; PRN, GBB, BRG: Ophiolites (serpentinites, gabbros, basalts); "S. Fiora" Unit; MLL: "Monte Morello" Fm.; AVR: "Santa Fiora" Fm.; PTF: "Pietraforte" Fm.; Subligurian Unit; ACC: "Canetolo" Fm.; Tuscan Nappe: MAC: "Macigno" Fm.; STO: "Scaglia Toscana" Fm.; MAI: "Maiolica" Fm.; DSD: "Diaspri" Fm.; POD: "Marne a Posidonia" Fm.; LIM: "Calcare Selcifero" Fm.; RSA: "Calcare Rosso Ammonitico" Fm.; MAS: "Calcare Massiccio" Fm.; RET: "Calcare a Rhaetavicula contorta" Fm.; CCA: "Burano and Calcare cavernoso" Fm.; Tuscan metamorphic units; "Monticiano-Roccastrada" Unit (MRU); MRU3: "Verrucano" Group; MRU2: Palaeozoic Phyllite-Quartzite Group.

Also the well stratigraphy was generalized on the basis of the simplified stratigraphy described above (see Figure 4 right). As regard the data from geophysical investigations, they were implemented in the geodatabase as follows:

Interpreted 2D reflection seismic data were added to geodatabase and used for constraining geometries and characteristics of the subsoil units. In Figure 3 traces of interpreted seismic lines, derived from [4–6,10,13,19–22,29,30], are reported. In particular, seismic lines define the structural setting of the main geological units, evidencing reverse faults (main thrusts) involved in Miocene compressive phase of North Apennine chain formation, normal faults formed during the successive Pliocene-Pleistocene tensional phase.

Geoelectrical data, covering wide areas of the volcano and carried out during the early exploration phase of the Mt Amiata geothermal fields [3], are used to define geometries and thickness of resistive substratum, coinciding with carbonate and evaporite formations (shallow geothermal reservoir).

**Figure 5.** Fence diagram built from geological cross sections (vertical dimension is exaggerated by factor of 2). *Z*-axis in meters; *X*-axis and *Y*-axis: the coordinate system used is WGS 84/UTM zone 32N (EPSG: 32632). Arrow points to the North.

Gravity data reported in the Bouger anomaly map were elaborated through 2.5 D numerical modelling and implemented in geodatabase to constrain geometries and thickness of subsoil formations, including deep intrusive low-density bodies and thickness of Neogene sedimentary basin. All geophysical interpretations and modelling were carefully calibrated using deep and shallow wells, previously implemented in our geodatabase.

The 3D geological model is 28 km × 20 km wide, with its longer wide oriented EW, and has a vertical extension of approximately 4 km. A 3D grid, which forms the skeleton of the geological model, was then generated and used for the property modelling stage, in which 3D surface (such as lithology, resistivity and temperature) were created.

The main horizons, corresponding to the top of the simplified geological units discussed above, were obtained by using the "convergent interpolation algorithm" built in Petrel platform (*i.e.*, see Figure 6).

**Figure 6.** 3D view of geological map spread on DEM and of the horizon corresponding to the top of the first geothermal reservoir in the Mt Amiata area (vertical dimension is exaggerated by factor of 2). *Z*-axis in meters; *X*-axis and *Y*-axis: the coordinate system used is WGS 84/UTM zone 32N (EPSG: 32632). Arrow points to the North. Legend as reported in Figures 1 and 4.

Shallower horizons are more accurately located and modelled because more wells encounter them. The reconstruction of deeper horizons is more problematic because there are fewer data; obviously, the horizons are more accurately modelled in the areas of Bagnore and Piancastagnaio active geothermal fields where density of deep wells and geophysical data is considerably higher compared to other areas (Figures 1, 2 and 6). The outputs of the 3D model contain information about depths of stratigraphic units and other geological information about faults, erosive surfaces and geothermal information such as reservoir volumes, used for geothermal potential assessment, and to constrain the conceptual model of the geothermal system. In correspondence of the Mt. Amiata volcano, geoelectric survey data were used to reconstruct surfaces between cover and reservoir units. Seismic reflection lines and gravimetric modelling are used to constrain the surface shapes where deep wells are not present. The top surface of Ligurian cover unit is well constrained in correspondence of geothermal fields (Figure 7) where stratigraphy of deep boreholes occurs. Particular attention was also dedicated to the reconstruction of surfaces delimitating the shallow reservoir and the metamorphic basement (Figure 7). The abrupt deepening of all surfaces toward East reflects the occurrence of the normal master fault, delimitating the Western side of Radicofani Pliocene sedimentary basin (Figures 1, 6 and 7). In this particular situation low angle normal fault offset reaches 200–400 m. Radicofani basin (eastward) have about 2400 m of sediments in its depocentre, as indicates by deep wells.

**Figure 7.** 3D subsurface horizons of the study area that are incorporated in the model, as reconstructed from available surface and subsurface data (vertical dimension is exaggerated by factor of 2). In figure three of these modelled horizons are shown: (**a**) green = top of cover formations (Ligurian units); (**b**) Blue = top of shallow geothermal reservoir; (**c**) Purple = top of metamorphic basement (see Figure 4 for details).

**Figure 7.** *Cont.*

#### **4. Numerical Modelling**

The output of the 3D geological modelling was used as input data for the numerical simulation of Mt. Amiata hydrothermal system. The geometry of the top of deeper geothermal reservoir, which is hosted within the Palaeozoic phylladic basement (unit 6 of Figure 4), was defined on the basis of the boreholes that encountered it. No data constrain the base of this reservoir; therefore in the numerical model we assumed a thickness of 500 m.

The Mt. Amiata hydrothermal system was simulated using the numerical simulator TOUGH2 (Lawrence Livermore Laboratories, Berkeley, CA, USA), a general-purpose code for modelling multi-dimensional, multiphase/multi-component flow and heat transport in porous and fractured media [36]. The modelling was performed using the "equation of state" module (EOS1) for pure water. Corey curves were adopted for the relative permeability. Capillarity, adsorption, and double porosity were not considered. The simulation domain has an area of 23 km × 17 km (391 km2 , Figure 1) and a depth of 5500 m below sea level (*i.e.*, the model height ranges from 7160 to 5750 m, considering local topographic relief). The model is centred on the active geothermal fields where available data allow constraining a reliable physical model. The numerical grid is subdivided into 22 layers and consists of 8228 cells (Figure 8); each horizontal layer consists of 374 (22 × 17) cells with variable sizes. A cell size of 1 km × 1 km is used. Along the domain boundary, where great detail is not necessary, cell size of 2 km × 2 km has been adopted. The layer thickness varies from a maximum of a few thousand meters to only 60 m as a function of the horizons derived from 3D geological modelling.

The physical parameters used in the model (permeability, porosity, density, thermal conductivity under saturated condition, specific heat) are listed in Table 1 and based on literature values [37–40].

**Figure 8.** 3D simulation grid. Cap rock corresponds to units 2, 3, 4 of Figure 4; Shallow reservoir corresponds to unit 5 of Figure 4; Phylladic basement, Deep reservoir and Basement rock correspond to unit 6 of Figure 4.

**Table 1.** Petrophysical parameters of main rock types.

Temperature and pressure are considered time invariant in the cells at the top and bottom of the simulation grid. A fixed temperature of 10 °C and atmospheric pressure were set for cells along the upper boundary. A variable temperature between 400 and 500 °C was assumed for the cells in the bottom layer, in accordance with the temperatures assumed for K-horizon. The lateral boundaries are considered no-flow insulating boundaries. The simulation of the natural state of the geothermal system (steady state) covers a period of 200,000 years but in some simulations the steady state was reached before. During the simulation some parameters (mainly permeability and temperatures at the bottom layer) were tuned to match the temperature distribution at depth. To verify the reliability of the simulated values of temperatures, simulated well temperature profiles were compared with temperature profiles of wells drilled in the area. A good match between simulated and measured temperature in deep wells was obtained (Figure 9), suggesting that a satisfactory simulation was achieved. The shape of the profiles from temperature modelling is linear or stepped. A vertical shape of the profiles indicates the occurrence of a convective heat transfer (reservoir zones), whereas an oblique shape is suggestive of a mainly conductive heat transfer (cover and impermeable units). The simulated temperature profiles of PC 30, BAGN 13 bis and BAGN 3 bis deep boreholes are indicative on the occurrence of both geothermal reservoirs of Mt. Amiata area (deep and shallow) and this is in agreement with the well logs. The simulated temperature profile of Nibbio 8 borehole is in good agreement with the measured temperatures at the level of the deeper reservoir whereas it does not agree at the level of the shallower one. The simulated profile suggests a convective heat transfer whereas the measured temperatures indicate a mainly conductive heat transfer. This disagreement is probably due to scarce permeability of the carbonate units of the Tuscan Nappe in the area of Nibbio 8 well; this would locally prevent the efficient hydrothermal fluid circulation and the consequent convective heat transfer. Capannacce 1 and PC 24 boreholes reach only the shallower geothermal reservoir.

**Figure 9.** Simulated and observed temperature vertical profiles along hole for selected deep geothermal wells in Bagnore and Piancastagnaio active geothermal fields of the Mt. Amiata area. On the left of each temperature profile the simplify stratigraphy (see Figure 4) of the well is reported.

The main results of numerical modelling are shown in Figure 10. At the level of first geothermal reservoir (0 m sea level), numerical modelling results evidence the occurrence of three high temperature zones centred in the Piancastagnaio (East), Bagnore (West) and Nibbio (South-East) geothermal areas. The higher temperatures observed in the Piancastagnaio area (East) are due to the normal fault system that characterizes the Eastern sector of Mt. Amiata. On the other hand, the highest temperature and highest productivity of the deep boreholes recorded in the area are found in the Piancastagnaio geothermal field. At level of the second and deeper geothermal reservoir (í2000, í3000 m sea level) the only two significant high temperature zones occur in correspondence of the Piancastagnaio and Bagnore geothermal fields (Figure 10). This agrees with the temperature measured in deep geothermal boreholes.

**Figure 10.** Thermal state of Mt. Amiata area at 0 m sea level (sl) (**a**), í1000 m sl (**b**), í2000 m sl (**c**), í3000 m sl (**d**) respectively, as deduced from numerical modelling. Dotted circles indicates the two active industrially exploited geothermal fields of Bagnore (left) and Piancastagnaio (right).

**Figure 10.** *Cont.*

#### **5. Discussion and Conclusions**

The use of a 3D subsurface modelling represents an effective aid to better comprehend the geologic framework of the Mt. Amiata geothermal system. After combining the information, we received a 3D geological/stratigraphic model and outlined a conceptual model of the geothermal system of Mt. Amiata area.

The Mt. Amiata geothermal system is formed by two distinct water dominated reservoirs [3,13,41]. The shallower one is hosted in the Mesozoic, carbonate-evaporitic formations of the Tuscan Nappe at depths ranging between 500–1500 m, the second, and deeper one in the Palaeozoic, metamorphic basement at depths below 2500–3000 m. These two reservoirs are characterized by different temperatures (150–220 °C for the shallower one, 300–350 °C for the deeper one) and, although they are in hydrostatic equilibrium, they are separated by a low permeable layer as evidenced by the temperature profiles of deep wells (Figure 9).

3D integrated modelling allows the reconstruction of the shallow geothermal reservoir geometries (Figure 7). Figure 7a evidences the irregular shape of the top surface of Ligurian cover units and highlights the sharp deepening of this surface that occurs in correspondence of the western master fault of the Radicofani basin. In Figure 7b the top surface of carbonate fractured shallow reservoir is shown. This has a complex geometry due to the tectonic framework of the area, particularly characterized by an important uplift caused by the emplacement (between 130 and 300 ka) of the low-density intrusion that constituted the volcano feeding system. In the northern side of Mt Amiata, the reservoir outcrops in Bagni S. Filippo area (Figure 1), this represents one of the recharge zones of the geothermal fields. The carbonate units deepen abruptly (more than 30° dip westward) reaching about 1000 m of deepening westward (Figure 11). Shallow reservoir outcrops also in Selvena-Castell'Azzara hills and deep northward toward the volcano edifice (Figures 7b and 11). This represents another probable recharge area of the geothermal systems. The top surface of the reservoir shows structural highs in correspondence with the Piancastagnaio, Bagnore and Poggio Nibbio geothermal fields and Bagni S. Filippo low-temperature hydrothermal system, where thermal and gas pools are widespread. These highs correspond to highs either in the Palaeozoic basement and intrusive body and K-horizon. The 3D geological model takes in evidence that the thickness of the shallow geothermal reservoir is strongly variable (150–1000 m) with the thinner portions that often correspond to the structural highs (Figure 7b). In Figure 7c the top surface of the metamorphic basement is represented. The dome shape of the metamorphic basement influences the thickness of the carbonate-evaporitic shallow reservoir that is laminated in correspondence of the metamorphic basement highs. This is due to the displacement of the ductile carbonate-evaporitic units of the Tuscan Nappe as a consequence of the emplacement and uplift of the low-density intrusion below Mt. Amiata volcano.

**Figure 11.** Isobaths (m sl) of the top of the shallow geothermal reservoir in the Mt. Amiata area.

It is impossible to reconstruct with the same accuracy the top of the deeper geothermal reservoir. This is due not only to the paucity of borehole data that encountered this horizon but also because the top of the deeper reservoir does not correspond with a specific geologic unit but it depends on the occurrence of thermometamorphic and metasomatic processes that affected the phylladic basement in the past, modifying its rheology. This is analogous to what observed in the similar Larderello-Travale geothermal system [16].

The occurrence of the active geothermal systems in the south-central portion of the Mt. Amiata may suggest that the intrusion (the heat source of the geothermal systems) is probably asymmetric with respect to the Mt. Amiata volcano. Further evidence supporting this is represented by the occurrence of considerable mercury mineralization [12,42] that is linked to the hydrothermal activity of the area and mainly distributed in the central-southern part of Mt. Amiata (Figure 1).

The subsequent combination of borehole data and numerical modelling results (in the areas where borehole data of temperature are not present), implemented in the Petrel geological model platform, allows reconstructing the temperature behaviour at the top of the shallower geothermal reservoir that is an industrial target for geothermoelectric power production (Figure 12). This is an important issue accomplished thank to this approach. In previous studies on Mt. Amiata geothermal system [13], numerical modelling experiments reconstructed the isotherms along horizontal slices in the subsoil of the geothermal system, as we also carried out in Figure 10. This of course does not take into account the 3D behaviour of geological units that are cut by horizontal slices (Figure 10). The isotherms reported in Figure 12 show the highest values in correspondence of Piancastagnaio area, where they reach temperatures in excess of 230 °C, and Bagnore area where temperatures up to 180 °C are shown. These very high temperatures could be related to highs of the intrusive body and, particularly in the Piancastagnaio area, to the fractured gauge zone of the normal master fault delimiting the western side of the Radicofani basin, which favours the rise of hot fluids from depth toward the shallow geothermal system.

Figure 12 shows that, in the northern and southern portion of the area, the temperatures sharply decrease. This occurs in correspondence with the outcropping of the evaporite-carbonatic units of the Tuscan Nappe (Figures 1 and 11) that form the shallow geothermal reservoir. This is in agreement with the fact that these outcrops are assumed as the main hydrogeological recharge areas of the Mt. Amiata geothermal system.

The obtained 3D geothermal model of Mt. Amiata area can be used for a first theoretical evaluation of the geothermal potential of the shallow reservoir (where we have more data for constraining its volume) [43]. As an output of the 3D model, a volume of 293 km3 has been estimated for the shallow geothermal reservoir, considering the area of the numerical modelling simulation (23 km × 17 km), and an average temperature of the reservoir of 148 °C. The estimate of geothermal energy stored underground above the mean annual temperature (15 °C), obtained by the volume method [43], is 2840 GWy (thermal). Assuming a mean value of the recovery factor of 10%, the exploitable energy of the first reservoir can be quantified in 284 GWy.

The quality of the elaborated model and the possibility of using it for prediction purposes will be improved when data from new deep wells will be available. An integrated approach with 3D tools makes geologic information more accessible to a wider audience of non-geologic specialists and decision-makers. The model is also suitable for the planning of new boreholes and further exploration and exploitation steps to proceeds with the usage of geothermal energy in Mt. Amiata area.

**Figure 12.** Map of the isotherms at the top of the shallower geothermal reservoir as elaborated by the combining of 3D geological and numerical modelling (up). Map of the isobaths of the top of the shallower geothermal reservoir from 3D geological model (down).

#### **Acknowledgments**

We thank SORGENIA GEOTHERMAL Srl (Milan, Italy) for access to the data and for permission to publish them and Matteo Ceroti and Sefano Scazzola for encouraging this work. We are in debt with two anonymous referees whose comments and suggestions significantly contribute to improve the quality of this work.

#### **Author Contributions**

Alessandro Sbrana had the original idea for the study and, with all co-authors contributed to the conceptualization and writing of the paper.

#### **Conflicts of Interest**

The authors declare no conflict of interest.

#### **References**

### **The Marsili Volcanic Seamount (Southern Tyrrhenian Sea): A Potential Offshore Geothermal Resource**

#### **Francesco Italiano, Angelo De Santis, Paolo Favali, Mario Luigi Rainone, Sergio Rusi and Patrizio Signanini**

**Abstract:** Italy has a strong geothermal potential for power generation, although, at present, the only two geothermal fields being exploited are Larderello-Travale/Radicondoli and Mt. Amiata in the Tyrrhenian pre-Apennine volcanic district of Southern Tuscany. A new target for geothermal exploration and exploitation in Italy is represented by the Southern Tyrrhenian submarine volcanic district, a geologically young basin (Upper Pliocene-Pleistocene) characterised by tectonic extension where many seamounts have developed. Heat-flow data from that area show significant anomalies comparable to those of onshore geothermal fields. Fractured basaltic rocks facilitate seawater infiltration and circulation of hot water chemically altered by rock/water interactions, as shown by the widespread presence of hydrothermal deposits. The persistence of active hydrothermal activity is consistently shown by many different sources of evidence, including: heat-flow data, gravity and magnetic anomalies, widespread presence of hydrothermal-derived gases (CO2, CO, CH4), 3 He/4 He isotopic ratios, as well as broadband OBS/H seismological information, which demonstrates persistence of volcano-tectonic events and High Frequency Tremor (HFT). The Marsili and Tyrrhenian seamounts are thus an important—and likely long-lasting-renewable energy resource. This raises the possibility of future development of the world's first offshore geothermal power plant.

Reprinted from *Energies*. Cite as: Italiano, F.; De Santis, A.; Favali, P.; Rainone, M.L.; Rusi, S.; Signanini, P. The Marsili Volcanic Seamount (Southern Tyrrhenian Sea): A Potential Offshore Geothermal Resource. *Energies* **2014**, *7*, 4068-4086.

#### **1. Introduction**

A large amount of data on the geological, geophysical and geochemical features of submarine volcanic activity and seamounts have been acquired in the last three decades. The close proximity of magma chambers to the seafloor, in conjunction with tectonic activity due to plate motion, deformation and cooling of erupted lavas results in convective circulation of dense, cold seawater through the cracked and fissured upper portions of the lithosphere; this circulation promotes the formation of venting sites that release hot hydrothermal fluids and dissolved elements [1–4].

Submarine hydrothermal activity has been studied so far as an energy source for free-living and symbiotic chemosynthetic bacteria, which form the base of the food chain in these unique habitats [5]. Because of their huge, long-lasting recharge and high-temperature/high-pressure characteristics, the submarine hydrothermal fluids are now investigated also as a potentially exploitable geothermal energy source. Candidate areas for offshore geothermal exploitation have been identified in the Gulf of California, the Juan de Fuca Ridge, the Japan Sea, the Okhotsk Sea, the Andaman Sea and the Tyrrhenian Sea since late 1970s [6]. However, the technology at that time was not still achievable at competitive costs beyond environmental, legal and institutional problems which had to be overcome. Nowadays, the continuous and growing developments in oil and gas exploration and exploitation techniques allow for an easier and economically competitive approach to the investigation and the energy potential assessment of submarine hydrothermal systems. However, before any reliable quantification for energetic exploitation of any offshore geothermal reservoir, a multidisciplinary submarine exploration has to be done, including geological, geophysical geochemical, oceanographic and biological investigations.

At the beginning of the last century, Italy became the first country worldwide to exploit high-temperature geothermal resources for power generation, with the result that, nowadays, Larderello-Travale/Radicondoli and Mt. Amiata (Bagnore and Piancastagnaio) are exploited geothermal fields [7]. They are both set in a continental extensional tectonic environment (lithospheric thickness between 20 and 30 km), characterised by deep and shallow volcanic systems and high heat-flow values (regional value of 120 mW/m2 , with maxima up to 1 W/m2 [8]). In both the geothermal areas, magmatic bodies provide the necessary heat source to deep and shallow reservoirs hosted in local metamorphic, carbonate and anhydritic formations [9,10]. In order to improve the Italian geothermal energy production, in 1970s and 1980s an intensive exploration program (geological, geophysical and geochemical surveys, as well as drilling activities) took place in Latium (about 100-deep wells were drilled in the Latera caldera, Vico Lake, Cesano, Bracciano Lake and Alban Hills), Campania (Phlaegrean Fields and Ischia Island) and Sicily (Vulcano and Pantelleria Islands) [7]. Several unfavourable physical, chemical and logistic features strongly limited the use of these potential geothermal fields. In particular, the geothermal exploitation was hindered by the low permeability of the reservoirs, the high salinity and acidity of the hot fluids, the sluggish and low amount of groundwater recharge, the high potential explosivity of these magmas (rhyolites, trachytes and phonolites) implying a high volcanic hazard as well as the strong urbanisation and/or tourist use of those areas.

The Southern Tyrrhenian Sea might represent the future target for geothermal energy exploration and exploitation due to the generally high heat-flow [11,12] and the widespread presence of seamounts for which the presence of hydrothermal deposits has been well documented. Among them, the Marsili seamount exhibits some features that make it a potential site hosting exploitable submarine geothermal systems.

This paper accounts for the recently collected information of geological, geochemical and geophysical features of the Marsili seamount (Figure 1), the largest European and Mediterranean volcanic edifice, making the seamount a likely large geothermal energy resource that could significantly improve the Italian geothermal power generation, nowadays providing 5.5 TWh per year that cover only 1.6% of national mean electricity production [13].

**Figure 1.** Location of the study area. MV: Marsili Volcano; MB Marsili Basin; ST: volcanoes of Stromboli; VU: volcanoes of Vulcano; FL: Filicudi Island; SIS: Sisifo submarine Volcanoes; LM: Lametini Seamounts; PA: Palinuro Seamount; VB: Vavilov Basin. The map on the right is modified after [14].

#### **2. Results**

This section summarizes the results of a large number of multidisciplinary data collected by several scientific cruises carried out by Italian Research Institutions during the last 20 years and describes the new results collected during two campaigns carried out in 2007 and 2010.

#### *2.1. Geothermal Potential of Southern Tyrrhenian Basin and Marsili Seamount*

The Southern Tyrrhenian Sea is a back-arc basin developed from Miocene to Present in the frame of the coeval formation of Apennine-Maghrebides chain, structured above the subducting north-western Ionian oceanic slab [15,16]. Its evolution has been characterised by volcanic activity induced by a wide tectonic extension from North-West (Sardinia) to South-East (Aeolian Arc) [17–19]; deep and shallow seismicity occurring in this area has been also well documented (e.g., [20–22]). The best morphologic evidences of these geodynamic processes are the two oceanic crust floored Vavilov and Marsili sub-basins and the homonymous seamounts, placed at north-west and south-east portions of the Southern Tyrrhenian basin, respectively [14,23]. Other geodynamic-related morphologic features are the numerous volcanic structures forming an arc all along the south-eastern margin of the basin: the Aeolian Arc, representing the emerged part; Palinuro, Glabro, Alcione and Lametini, representing the northern submerged arm; Eolo, Enarete and Sisifo, representing the north-western submerged arm (Figure 2; [23]).

**Figure 2.** Bathymetry of South-Eastern Tyrrhenian Basin. M: Marsili; MB: Marsili Basin; P: Palinuro; A: Alcione; L: Lametini; Sc: Stromboli Canyon; Eo: Eolo; En: Enarete; S: Sisifo (modified after [23]).

The heat-flow anomalies of Southern Tyrrhenian Sea [11,12] are always high in correspondence with the volcanic structures; the regional mean value is higher than 120 mW/m2 and comparable to the Italian onshore volcanic district. More precisely, the highest values have been recorded close to the Vavilov (140 mW/m2 ) and Marsili (250 mW/m2 ) seamounts. Heat-flow rates as high as 300 mW/m2 and 500 mW/m2 were detected on the uppermost and central portion of Vavilov and Marsili seamounts respectively [24]. Those values fit with the geophysical data, collected in the last decades such as:

All the above-mentioned data strongly suggest the following inferences:

• A low explosivity due to the presence of lesser involved magmas with a lower amount of dissolved water is expected.

In this frame, Marsili seamount, presenting the highest aforementioned anomalies, is here proposed as the best prone site of sustained submarine geothermal potential. Marsili seamount is a recent volcanic structure (1–0.1 Ma), mainly composed of basalts and, to a lesser extent, andesites and trachy-andesites with calc-alkaline affinity [32–34]. This structure has been interpreted as an inflated, small-scale spreading centre, since bio- and magneto-stratigraphic data show features comparable to those typical of mid-ocean ridges [35–38].

Dekov and Savelli [39], summarising all the observations since the late 1960s on the rocks sampled from Marsili and surrounding volcanic centres, affirmed that these areas are affected by hydrothermal fluid circulation; in their model, cold seawater enters into fractured rocks and then is superheated by magmatic bodies at crustal depths. Clues of active venting of hydrothermal fluids over the Marsili seamount have been provided by the injection of magmatic-type volatiles into the sea water as demonstrated by the anomalous 3 He content and the whole chemical composition of the gases dissolved in water column [40]. However, direct thermal measurements on Marsili seamount are still lacking. This latter discovery strongly points to Marsili still being hydrothermally active, in agreement with the previous mineralogical data, and supported by a significant contribution by juvenile fluids.

#### *2.2. Morphology*

The physiographic features of Marsili seamount are shown in Figure 3. The Marsili seamount rises 3500 m from the abyssal plain to 489 m minimum depth. The volcanic edifice is about 60 km long and about 20 km wide; two other smaller structures run parallel to the volcano and are located on its left and right part, with heights on the order of several hundred metres. Marsili is elongated mainly along a NNE-SSW axis; however, this extensional axis is not perfectly linear but it shows a sigmoid trend, with the southernmost and northernmost axial directions both north-eastern trending, while the central axis is closer to the N-S direction. This deformational feature is commonly associated with strike-slip faulting under volcanoes [41], but it can also be related to long kilometric normal faults orthogonal to the volcanic system spreading direction, with transverse faults at the tips that accommodate the regional strain [42]. In both cases, this observation implies a strong tectonic control of the whole volcanic edifice. The volcano summit is characterised by a narrow crest, 20 km long and 1-km wide, over the 1000 m isobaths, cut by linear structures, mainly disposed parallel to the extensional axis. They appear as small ridges, up to 100 m high and 750 m long, not fault bounded, formed by alignments of small monogenic volcanic mounds. These segments may result from the feeding of magma to the seafloor along dikes that produce, or follow, discontinuous, en-echelon crack systems. Steep bathymetric gradients separate the crest portions from the deeper volcano flanks, forming very long and lesser sloped scarps extending to the basin. In several places, these lower steep scarps terminate with gently dipping terraces, elevated for several hundreds of meters from the abyssal Marsili basin. Spherical cones are present on the flanks of the volcano with diameters ranging between a few kilometres to several hundred metres.

These flanks are also cut by several, well evidenced kilometric valleys; the largest of them is located near to the central and top portion of Marsili. It has an amphitheatre shape in the uppermost part edged by roughly vertical walls; such geomorphologic feature suggests a flank collapse of this portion of the Marsili flank, with landslide debris accumulated distally on the Marsili basin, a kind of process already found elsewhere [43,44].

**Figure 3.** Detailed bathymetry of Marsili seamount. Red dashed lines: linear structures; violet dotted lines: main circular cones and terraces; yellow dashed line: major landslide (modified after [45]).

#### *2.3. Magnetic Data*

Figure 4 shows the total field magnetic anomaly map, reduced to the magnetic pole, recorded above Marsili seamount. In agreement with literature data [35,37,38], positive magnetic anomaly maxima are located along the central sectors of the volcanic structure. At northern and southern tips the highest positive anomalies of Marsili are present, with maximum values around 1500 nT; in the central and highest portion of volcano and in the middle of these two former anomalies, the values are around 0 nT, reaching a minimum of í100 nT. Negative anomalies are located at the base of the western flank, with a mean value around í200 nT; the base of the eastern flank again shows negative values, achieving a local minimum of í500 nT in the central and southern portions, whereas the north-eastern part of this sector is characterised by small positive values.

**Figure 4.** Magnetic anomaly field (reduced to magnetic pole) of Marsili; values of the magnetic anomaly field in nT (modified after [45]).

#### *2.4. Gravity Data*

According to the map of Faye gravity anomalies of Marsili seamount (Figure 5), positive maxima, reaching more than 80 mGals, are present on the crest portion of volcano, above the 1000 m isobaths, and progressively decrease to about 20 mGals on the volcano base. On the whole, the sea-surface free-air anomaly field largely reflects the topography of the volcano; moreover, the relative low amplitude of the positive gravity anomalies suggests the presence of a low density zone below the seamount.

#### *2.5. Seismological Recordings*

The broadband OBS/H (Ocean Bottom Seismometer with Hydrophone), deployed on the Marsili flat top (39°16,383', 14°23,588') at a depth of 790 m, recorded more than 1000 seismo-volcanic and hydrothermal signals [3]. By comparing the signals recorded with typical volcanic seismic activity, the recorded signals were grouped into: Volcano-Tectonic type A (1 event), Volcano-Tectonic type B (817 event), High Frequency Tremor (159 events) and quasi-monochromatic Short Duration Events (32 event).

**Figure 5.** Faye anomaly field of Marsili seamount; values of the gravity field in mGal (modified after [45]).

An intriguing feature of the seismic noise is its spectral content presenting progressively growing energy levels in a broad-band frequency range from 4 to 60 Hz; diffuse spectral peaks are also present between 4 and 30 Hz. The seismic events are characterized by an indiscernible S-phase. On the basis of preliminary frequency content observations they can be distinguished in two main groups: about 720 events with frequency content between 4 and 10 Hz and about 80 events with very high frequency content, between 40 and 80 Hz [3].

#### *2.6. Hydrothermal Fluids*

In order to recover information on the type of vented fluids and to constrain their origin, we carried out water column surveys performing casts and tows across the seamount using a rosette equipped with a CDT (SeaBird 911+ Conductivity Temperature Depth device) and Niskin bottles. The surveys have been carried out during two scientific cruises in November 2007 (by R/V Urania) and July 2010 (by R/V Astrea).

Sea-water samples were collected to extract the dissolved gas phase by further laboratory procedures. The extracted gases were analyzed for the chemical composition calculated taking into account the solubility coefficients (Bunsen coefficients "ȕ") of each gas specie, the volume of the gas extracted (in cm3 ), the volume of the extracted water sample (in cm3 ) and the equilibration temperature as described in [46,47].

Table 1 lists the results (expressed in cm3 STP/LH2O, namely millilitres of gas per litre of water at standard temperature and pressure conditions of 1bar and 25 °C) showing that CO2, with concentrations ranging from 0.53 to 1.9 cm3 STP/LH2O, is by far the most abundant component of the dissolved gas phase besides oxygen and nitrogen, in the range 2.7–4 and 8.3–10 cm3 STP/LH2O, respectively. Among the other components CH4 is always detected in concentrations ranging from 3.9 × 10í5 to 2.1 × 10í4 cm3 STP/LH2O, two orders of magnitude above the equilibrium with the atmosphere (ASSW value equal to 1 × 10í6 cm3 STP/LH2O).

**Table 1.** Analytical results of the dissolved gas phase in samples collected during two different cruises (2007 and 2011). Hydrothermal-derived gases are released besides components from the atmosphere. Analytical results in cm3 STP/LH2O (see text for details). Depth in meters. Bdl = below detection limits. Data of ASSW (Air Saturated Sea Water) and from the bottom sample of a vertical cast carried out over the Tyrrhenian Abyssal plain are reported for comparison. Coordinates are reported as they were collected during the cruises.

#### **3. Discussion**

The collected data provide the first useful information on the characterisation of this huge offshore structure that might be a potential site for offshore geothermal exploitation. The large variety of morphological forms, *i.e*., linear structures, lava sheets, spherical cones, fault scarps and landslides, suggests that Marsili volcanic seamount has encompassed a long evolution, composed of different stages. The whole volcanic edifice shows evidence of mature volcanism, since almost all the volcanic forms look diffusely marked by further tectonics and erosion. The most mature forms, represented by circular cones and terraces, lava sheets, *etc.*, are predominantly located on the lower flanks of Marsili volcano; the youngest ones, typically represented by linear structures, are along the crest portion. Fault scarps are easily observable, pointing out the influence of tectonics on the whole edifice. As a consequence, the volcano likely has a well-developed network of fractures and cavities that play an important role in defining the permeability field of the volcano; these provide a pathway for the inflow of seawater and the consequent hydrothermal circulation and are primarily responsible for determining the magnitude of submarine geothermal resources.

The basaltic to basaltic-andesite rock compositions reported in [34] are in agreement with the strength of magnetic anomalies measured on the northern and southern portions of the Marsili volcanic edifice. Generally, freshly erupted basalts (also rapidly cooled) are strongly magnetic [48,49] because they have not been exposed to demagnetisation processes for a long time. However, the most intriguing feature of the magnetic anomaly field are the very low values in the central sector of the crest, testifying the presence of rocks with very low and/or a-magnetic properties, for which hydrothermal activity is considered as the main cause. The hydrothermal circulation of fluids may interact with the source rocks and reduce the magnetisation by breaking down the magnetic minerals [50,51]; this evidence corroborates with the occurrence of hydrothermal processes, possibly still active or recorded by geothermal deposits on the crest of Marsili volcano [39,52]. The pattern of the observed anomalies fits with a Curie isotherm mean depth located at around 4–5 km below Marsili crest, indicating a temperature of more than 600 °C at the volcano base and highlighting the possible presence of magmatic bodies. The high internal temperatures are also supported by the extremely high heat-flow measurements carried out on Marsili [11,24]. All those complementary data coherently suggest that Marsili seamount contains an intense and shallow heat source.

Although hydrothermal minerals are generally less dense than primary minerals, their effect on gravity is negligible. The free-air anomaly high observed on the top of Marsili seamount occurs over the topographic high over which the hydrothermal system is supposed to be located. So, the gravity anomalies are controlled by local geologic structure. Generally, the observed gravity anomaly field can be fitted only assuming a mean density of the Marsili volcanic structure of about 2 g/cm3 . Taking into account the petrographic features of Marsili rocks [34] as well as the magnetic data, such values can be attributed to rock porosity/permeability, possibly filled with aqueous and volatile phases. On this basis, it can be tentatively inferred that the Marsili volcano should have a significant porosity, possibly more than 10% by volume.

Caratori Tontini *et al.* [53] proposed a model for the volcano summit consisting of a large altered region, with a vanishing magnetisation of 0 A/m and a density of 2.0 g/cm3 , and an active, a-magnetic hot magma chamber 3 km below the summit, with density 2.3 g/cm3 based on volatile concentrations [54].

Seismic observations by means of the OBS/H deployed on the crest of Marsili seamount [3] fit with the presence of active hydrothermal discharge highlighted by the geochemical features of the dissolved gases (Table 1). The continuous high frequency seismic and acoustic noise recorded on the Marsili crest is in agreement with historic seismic noise recordings in geothermal areas. In particular, very shallow hydrothermal manifestations (active venting sites) usually produce signals with maximum energy in the frequency band over 50 Hz [55–57], as observed on Marsili seamount.

Moreover, the seismicity characterised by a frequency content between 4 and 10 Hz could be connected to hydrothermal fluids circulation, in analogy with long-period (LP) events (e.g., [58–60]), while those with 40–80 Hz frequency content could be generated by surface degassing phenomena ([3] and references therein).The results of the chemical and isotopic analyses of gas extracted from the sampled sea-waters allow us to confirm the existence of a dissolved gas phase different from the atmospheric. The results clearly show that besides the expected atmospheric gases, a significant content of CO2 marks all the gas samples. Figure 6 shows the contemporary presence of atmospheric components (represented by O2 and N2) as well as CO2 typically originated by an endogenic source.

**Figure 6.** O2-N2-CO2 ternary diagram. The samples plot on a straight mixing line between atmospheric components O2 and N2 and CO2 typically of magmatic/hydrothermal origin The Air Saturated Sea Water mark (grey star, ASSW) indicates where dissolved gases in ASSW (air-saturated sea water) sea water plot. Blue triangles = samples from above the area marked by physical anomalies (see text); grey crosses = samples from a tow-yow across the seamount; black filled circles = samples from a vertical cast over the Tyrrhenian batial plain. The numbers indicate the sampling depths.

The plot along the Air-Saturated Sea Water (ASSW)-CO2 mixing line highlights the injection of CO2-dominated volatiles. Samples marked by the blue triangles come from the area marked by the largest anomalies of physical parameters (magnetic, gravimetric, seismic). Those samples plot along a line with a constant CO2/N2 ratio showing and increasing trend of hydrothermal fluids. Data from waters sampled by a tow-yow across the whole seamount, show a lower, although significant CO2 content, due to their provenance from other areas clearly marked by the presence of hydrothermal activity although of less intensity. Figure 7 plots atmospheric gas (typically

represented by oxygen) together with CO2 and CH4 representative for hydrothermal-type components. The plot shows the occurrence of dissolution processes responsible for CO2 loss and enhancement of the less soluble species such as CH4. Both figures plot data collected during two cruises in 2007 (grey crosses) and 2011 (blue triangles). The latter samples have been collected over the area supposed to be the main geothermal source of the whole Marsili seamount where the low energy volcanic blasts might have occurred in recent times [61].

**Figure 7.** CH4-O2-CO2 triangular diagram. The plot shows the relative concentrations of hydrothermal-derived components (CO2-CH4) besides the atmospheric component (here represented by oxygen) normally dissolved in sea water. It is easy to recognize the large enrichment in deep-originated CO2 and CH4 to the respect to ASSW. The arrows show the trends of modification of the dissolved gas assemblage due to the injection of hydrothermal fluids and the gas/water interaction phenomena leading to the dissolution of the very soluble CO2 and enrichment of less soluble gas species like CH4. Symbols as in Figure 6.

The available scientific data set has been considered a good starting point for an Italian company (Eurobuilding S.p.A., Servigliano, Italy) to propose a cutting-edge project to drill the first offshore geothermal well on the Marsili seamount, [62,63]. The estimates they performed bring to an effective electrical power generation of 200 MW. The offshore well should be connected to steam turbines that together with condenser systems, power generators and tension elevators could be hosted on the power plant platform. A high capacity cable will drive the electricity from the platform to the power grid located on the Italian coast.

#### **4. Conclusions**

The comprehensive and multidisciplinary (geological, geophysical and geochemical) dataset provided results coherently pointing to the Marsili seamount as a possible target for offshore geothermal exploration and exploitation developments.

The main scientific key points supporting future evaluations of the Marsili seamount as a possible geothermal energy resource can be summarized as:

The possibility that a large geothermal energy reservoir exists about 100 km off the Italian coasts and that such a reservoir is able to provide a significant, long-lasting amount of exploitable energy, makes the Marsili seamount a potentially attractive target for private and public investors.

To summarize the case for developing the world's first offshore geothermal power plant in this region we consider a model of the edifice with constant average density of basalts equal to 2.67 g/cm3 , namely the difference between the gravimetric anomaly calculated from the model and that coming from the measurement leads to a density value between 1.7 and 2.3 g/cm3 . Considering also the petrologic and magnetic features of the rocks, the estimated density values can be a consequence of a high porosity of the rocks filled by volatiles. Assuming an average porosity of 10% and the total estimated volume of the geothermal reservoir >100 km3 , the volume of the geothermal fluids circulating in the central and upper part of the Marsili seamount, it results that about 10 km3 of exploitable hot fluids can be recovered [45].

A first exploration drilling is needed to verify the proposed results and to recover (for the first time) direct information on the physical, chemical and isotopic features of the deep hydrothermal vents. In fact a hydrothermal fluid with a temperature of 250 °C and pressure of 30 bar with a flow rate of 1 m3 /s allows the energy production of 8 MW, while fluids marked by a temperature of 375 °C and pressure of 250 bar allow an energy production ten times bigger for the same flow rate.

It is absolutely clear that new data and further multidisciplinary investigations are needed, including the assessment of the natural risks (volcanic and seismic) related the huge and almost unknown submarine volcano. The development of micro-biological studies as well as research and development activities aimed at planning and deploying seafloor observatories for long-term monitoring activities [65] are also considered as mandatory. The results will improve the scientific knowledge useful to better constrain the available and exploitable geothermal resources inside the volcano. The future geothermal explorations have at least to include:

It is noteworthy to consider that, although the offshore geothermal energy exploration and exploitation have not yet been considered a feasible option, nowadays the technologies required for on-site geothermal exploitation and production of electric energy are already available, and can be adopted in view of the depth of seawater, about 600–800 m, and the minimum distance from the Italian coasts, less than 100 km. Presently, some projects are evaluated in European countries, for instance in Iceland, to develop offshore geothermal power plants to use high-temperature fluids lying along the mid-ocean ridges. Iceland is also developing the world's longest subsea power cable (about 1000 km), to export geothermal energy to Scotland. An eventual Marsili submarine power cable would "only" be around 150 km long.

Once the potentially exploitable geothermal energy of the Marsili seamount is well-constrained, that submarine volcano may become home to the world's first offshore geothermal power plant, potentially capable of doubling the current geothermal power output of Italy [45].

The increasing request for electrical power, the energy prices and the increasing know-how of the utilization of this resource makes the Marsili seamount an attractive scientific and industrial target. Several other seamounts of the Tyrrhenian sea, belonging to the same geological framework, share Marsili characteristics, among the most important: high heat-flow anomalies, relative shallow depth of the crest and rock porosity.

#### **Acknowledgments**

The Authors wish to thank the Eurobuilding S.p.A. for the fruitful cooperation during the data collection. The authors are indebted with Diego Paltrinieri for his comments and discussions during the development of the scientific activities.

#### **Author Contributions**

Patrizio Signanini and Paolo Favali were scientific coordinators of the geophysical projects; Mario L. Rainone and Sergio Rusi collected geophysical data during the 2006 scientific cruise and took care of the editing of the paper; Angelo De Santis e Francesco Italiano coordinated the geophysical and geochemical (respectively) data elaboration. All the authors cooperated to the writing of the paper.

#### **Conflicts of Interest**

The authors declare no conflict of interest.

#### **References**

### **Geothermal Potential Evaluation for Northern Chile and Suggestions for New Energy Plans**

#### **Monia Procesi**

**Abstract:** Chile is a country rich in natural resources, and it is the world's largest producer and exporter of copper. Mining is the main industry and is an essential part of the Chilean economy, but the country has limited indigenous fossil fuels—over 90% of the country's fossil fuels must be imported. The electricity market in Chile comprises two main independent systems: the Northern Interconnected Power Grid (SING) and the Central Interconnected Power Grid (SIC). Currently, the primary Chilean energy source is imported fossil fuels, whereas hydropower represents the main indigenous source. Other renewables such as wind, solar, biomass and geothermics are as yet poorly developed. Specifically, geothermal energy has not been exploited in Chile, but among all renewables it has the greatest potential. The transition from thermal power plants to renewable energy power plants is an important target for the Chilean Government in order to reduce dependence on imported fossil fuels. In this framework, the proposed study presents an evaluation of the geothermal potential for northern Chile in terms of power generation. The El Tatio, Surire, Puchuldiza, Orriputunco-Olca and Apacheta geothermal fields are considered for the analysis. The estimated electrical power is approximately 1300 MWe, and the energy supply is 10,200 GWh/year. This means that more than 30% of the SING energy could be provided from geothermal energy, reducing the dependence on imported fossil fuels, saving 8 Mton/year of CO2 and supplying the mining industry, which is Chile's primary energy user.

Reprinted from *Energies*. Cite as: Procesi, M. Geothermal Potential Evaluation for Northern Chile and Suggestions for New Energy Plans. *Energies* **2014**, *7*, 5444-5459.

#### **1. Introduction**

Chile is a republic located in South America bordered by the Andes Mountains and the Pacific Ocean. The country is partitioned into 15 administrative regions. Chile is rich in natural resources such as copper, timber, nitrates, precious metals and molybdenum. It is the world's largest producer and exporter of copper and mining is the primary Chilean industry and an essential component of the Chilean economy. Although Chile has an abundance of copper and other mining resources, it has limited indigenous fossil fuels and over 90% of its fossil fuels needs must be met by imports [1].

In 2012, Chile's energy demand was approximately 60,000 GWh and it is projected to increase by 4.9% annually reaching a consumption of 136,000 GWh/year by 2030 [1,2]. The installed capacity for electric generation is 17,500 MWe, producing approximately 65,600 GWh/year [3]. Of the total installed capacity 62% corresponds to fossil-fuelled power plants, 34% to hydropower plants and only 4% to renewable resources [3]. The majority of the total electricity is consumed by the mining and industry sectors. The electricity market in Chile is comprised of two main independent systems:

SING—the Northern Interconnected Power Grid;

SIC—the Central Interconnected Power Grid.

SING has an installed capacity of approximately 3800 MWe, 100% of which is generated by imported fossil fuel. SIC has an installed capacity of 13,500 MWe, 53% of which is generated by imported fossil fuels, 42% by hydroelectric power and the remaining 4% by renewable resources [3].

Therefore, Chile is heavily dependent on imported fossil fuels and hydropower represent the main indigenous energy source of the country. Other renewable resources such as wind, solar, biomass and geothermal energy are poorly developed, but several corporations, such as the Centre for Renewable Energy (CER), have been working to ensure the optimal participation of renewable energies in Chile's energy matrix to contribute to the sustainable development of the country [4].

In 2011, The Global Energy Network Institute (GENI) suggested that a strategic energy plan for Chile was necessary to ensure both energy autonomy of the country and the transition from thermal power plants to renewable energy power plants [4]. This passage from conventional energy to renewable sources is indispensable both to reduce the amount of greenhouse gases emitted into the atmosphere and to reduce the dependence on imported fossil fuels. Under such circumstances, renewable energy resources in Chile should become more relevant. Geothermal energy is not yet exploited here, but among all renewable resources it has the greatest development potential. Although no geothermal power plants have been installed to date, a vigorous geothermal exploration program is under way [5]. A preliminary evaluation of northern Chile's geothermal power potential is approximately 2000 MWe, whereas the country's central-southern region is estimated to be between 1000 and 1500 MWe [6,7]. This makes Chile one of the most attractive countries for the utilisation of geothermal energy.

The purpose of this work is to review the geothermal potential evaluation for northern Chile in terms of power generation, considering the conventional geothermal resources of Surire, Puchuldiza, Apacheta, El Tatio and Irriputunco-Olca. This evaluation could provide useful information for the development of a new strategic energy plan applicable to the SING. Northern Chile has been chosen for the analysis because more data are available than for the central and southern regions, and because of the country's large energy reliance on imported fossil fuels.

#### **2. Regional Setting of the Study Area**

Northern Chile has a relatively homogeneous geological setting, consisting of Lower Miocene–Pleistocene ignimbrite deposits and andesitic–rhyolitic volcanic products overlying Middle Cretaceous–Upper Miocene volcano-sedimentary formations [8–13]. The latter hosts the main hydrothermal reservoirs that predominantly consist of andesitic lava and pyroclastic flows, conglomerates, breccias, sandstones, siltstones, limestones, marls and evaporites [10,14]. The main hydrothermal systems are located within the NS-, NW-trending grabens [8,10,15] on the western side of the Pliocene–Holocene Central Andean Volcanic Zone.

For each studied system (Figure 1), a brief geographical/geological setting with relative descriptions of the geothermal manifestations at the surface is reported. The geothermal fields of Colpitas and Larima are also described, but it was not possible to evaluate their geothermal potential due to a lack of public data.

**Figure 1.** Geothermal areas of northern Chile. The red box shows the studied area.

#### *2.1. Colpitas*

The Colpitas geothermal field is located in the northernmost part of Chile, in the Arica and Parinacota region, with an elevation from 4000 to 5200 m above sea level (a.s.l.). Geologically, the area is characterised by volcanic rocks and volcanoclastic deposits and volcano-sedimentary sequences ranging from Miocene to Holocene [16]. The stratigraphy is mainly comprised of the Lupica Formation (Upper Oligocene-Miocene), corresponding to the *basement* of the basin. It is constituted by rhyolitic ash flow tuffs, andesitic lavas and subvolcanic plugs as well as epiclastic sandstones. The Lupica Formation underlies, in the western part of the study area, the Huaylas Formation (Upper Miocene). The Huaylas Formation is mainly formed by lacustrine epiclastic sandstones and gravels, with interbedded ash layers and is itself overlain by the Lauca Ignimbrite (Upper Pliocene).The potential reservoir is likely located within Cenozoic volcanoclastic rocks. The thermal springs, located in the northern and southern part of the Colpitas field, have temperatures that range from 28 to 55 °C with a total flow <10 L/s [16]. The North Thermal springs are interpreted as being the most representative of deeply derived thermal waters [16]. Stable isotope composition indicate that thermal spring waters have been subject to evaporation, or more likely, varying amounts of mixing with groundwater brine that underlines the salt deposits in the basin [16]. Na/K ratio of the North Thermal Springs gives the hottest equilibration temperatures of approximately 235 °C. Many of the springs have moderate bubbling of what is likely to be CO2 and H2S gas [16]. Recently, a slimhole well was drilled to prove the existence of a geothermal reservoir and to evaluate the potential of the field [3,16].

#### *2.2. Surire*

The Surire hydrothermal system is located at an altitude of 4000–4300 m a.s.l., in the southeast part of the Salar de Surire. Many volcanoes rise above 5500 m around the salar, some lavas have flowed into it and are now partially buried by sediments. The volcanoes are dacitic and andesitic; the dacities generally occur as domes and the andesities form stratovolcanoes. In 1972, as many as 133 thermal discharges occurred in an approximately 15 km2 area [15]. Presently, most of the bubbling pools and thermal springs are located along the southern border of the salar and have temperatures between 20 and 80 °C [15]. The Database of Geothermal Resources in Latin America & the Caribbean [17] indicated a potential temperature of the geothermal reservoir of 110–234 °C and an electric power output of 50–60 MWe.

#### *2.3. Puchuldiza*

The Puchuldiza-Tuja hydrothermal system is located at an altitude of approximately 4100–4200 m a.s.l. and is 27 km SW of the active Isluga volcano that is characterised by permanent fumarolic activity [18]. The hydrothermal area of Puchuldiza is limited both to the north and south by several Plio-Pleistocene stratovolcanoes reaching altitudes higher than 5000 m. A Pleistocene-Holocene fracture system characterized by northeast-southwest faults affected the Plio-Pleistocene volcanoes and produced, in part, the surficial hydrothermal activity. The fluid discharges within the Puchuldiza area are controlled by the Churicollo, Puchuldiza and Tuja faults. Several thermal springs with low gas emissions surround the main emission areas [19]. The geothermal reservoir is hosted in the Utayane Ignimbrite and the Puchuldiza Formation. The permeability of the volcanic formations is due to cooling joints and tectonic fractures [6]. Geothermometry evaluations suggest that the fluid reservoir has a relatively high equilibrium temperature, up to 270 °C [18]. The Database of Geothermal Resources in Latin America & the Caribbean [17] indicated a potential electric power output of 190 MWe.

#### *2.4. Lirima*

The Lirima geothermal field is located at an altitude of 3900 m a.s.l., 25 km SW of the Sillajguay volcanic chain. The geology of the area is characterized by Mesozoic basement rocks constituted by clastic-carbonate sequences, volcanic and sedimentary rocks of Oligocene-Miocene and middle Miocene to Pleistocene volcanic edifices [20]. In the Lirima area have been recognized bubbling pools, along the western side of the field, and three main sites with thermal springs [21]. The thermal springs are characterized by temperatures between 38 °C and 80 °C, high Cl and B concentration, į18O enrichment, and relatively low Mg concentration; consistent with deep circulation from a geothermal reservoir, and low mixing degree. Minimum temperatures from water and gas geothermometers range from 200 to 200 °C.

#### *2.5. Apacheta*

Apacheta is located 105 km NE of Calama City and 55 km NW of the El Tatio hydrothermal system. A 180 m deep well (PAE-1) drilled by the Chilean National Mining Company (CODELCO) in 1998 produced steam measured at 88 °C [22]. Fluid discharges emitting superheated steam (up to 118 °C) [22] with high flow rates are found along the eastern flank of the 5150 m high Apacheta volcano. Currently, in the Apacheta geothermal field, a project (Cerro Pabellón) is underway for power production from geothermal resources. Project feasibility studies began in 2005, then four wells were drilled from 2009 to 2010, whose depths reached between 1300 and 2000 m. The results of the production and injection tests showed the presence of a liquid-dominated reservoir with a maximum-measured temperature of 260 °C. ENEL Green Power has planned the drilling of 13 wells to operate a 50 MWe power plant. This project represents the first commercial-scale geothermal plant in the country [7].

#### *2.6. Irriputunco-Olca*

The Irriputuncu-Olca field is characterised by the presence of the Irriputuncu and Olca Volcanos, located in the Chilean Altipiano at 4000–5000 m a.s.l., and in the vicinity of the copper mine. Irriputuncu is an active dacitic stratovolcano, with fumaroles at the top crater and one acid-sulphate hot spring at the base of the volcano. Two slim boreholes (800 and 1430 m in depth) measured a bottom hole temperature close to 150 °C and 195 °C (at 3350 and 3000 m a.s.l., respectively) [23]. Time domain electromagnetic (TEM) and Magnetotelluric (MT) data suggests the presence of a potentially deeper reservoir at approximately 220 °C [23]. Olca is an andesitic volcano, of which the TEM-MT data exhibit two conductive layers intercalated with resistive zones. The possible thickness of the reservoir is 2000–3000 m and for the surface area there are three estimates, conservative (7.5 km2 ), likely (15 km2 ) and optimistic (45 km2 ) [23]. Preliminary results, assuming up to 10 MWe/km2 and a reservoir temperature of 230–300 °C, suggest a potential for electric generation between 75 and 450 MWe [23].

#### *2.7. El Tatio*

El Tatio is located 100 km E of Calama City at an altitude of 4300 m a.s.l. Several thermal springs, fumaroles, geysers and boiling and mud pools are present. Hydrogeological models [24,25] indicate that meteoric waters infiltrate in recharge areas 15 km E of the field. The main hydrothermal reservoir is confined within the permeable Puripicar Formation and the Salado Member. An important secondary aquifer occurs in the Tucle Dacite subunit that is capped by the impermeable Tatio Ignimbrite [24]. In Figure 2, the simplified geological profile and circulation conceptual model are shown. The potential geothermal reservoir is hosted in the Puripicar Formation and Salado Member, although a temperature of approximately 170 °C was recorded in the permeable levels hosted in the Grupo volcanic de Tucle [6].

In Figure 3, a profile of the El Tatio Graben is shown. It crosses through the wells numbered 1, 4, 9 and 7, from NW to SE. In the boreholes, three permeable zones were detected. The permeability essentially originated by tectonic fracturing or rapid cooling of the volcanic bodies [9]. The temperatures of the three permeable zones ranged from 170, 230 to 260 °C, moving from the Grupo volcanico de Tucle to Puripicar Formation and Salado Member, respectively [9].

**Figure 2.** Simplified geological map of the El Tatio geothermal field, geological profile and circulation conceptual model. The dashed box represents the area in Figure 3 (modified from [9]).

**Figure 3.** The geological profile from NW-SE of the El Tatio graben through the boreholes numbered 1, 4, 9 and 7 (modified from [9]). For further details see text in the previous page; paragraph relative to El-Tatio.

#### *2.8. Geothermal Exploration in Chile*

Initial geothermal exploration in the Central Andean Volcanic Zone took place in late 1960s in response to increasing Chilean energy demands. At El Tatio, a pre-feasibility investigation, funded in 1967 by the Corporation for the Promotion of Development and the United Nations Development Program (CORFO/UNDP), was followed by geological, geophysical and geochemical surveys from 1968 to 1980. Six 600 m deep exploration wells, drilled between 1969 and 1971, encountered temperatures up to 250 °C. Seven production wells drilled in 1973 and 1974 disclosed three discrete reservoirs with temperatures up to 260 °C. Three of these wells produced an average of 14.7 kg/s (adequate for 6 MW each); two other wells produced less, but could still be capable of 5 MW each. An electric power output of 100 MWe for the El Tatio geothermal field was estimated [26]. At Puchuldiza, geological, geochemical and geophysical studies were performed by CORFO/UNDP (from 1968 to 1974) and by Japan International Cooperation Agency (JICA) (from 1978 to 1980) to evaluate the geothermal potential. Six wells were drilled in 1976 and a depth of 1200 m was reached. Temperatures up to 175 °C were measured at depth of 900 m [27,28].

In the Surire zone, geological and geochemical [15,29] investigations were performed by CORFO between 1972 and 1979. Reservoir temperatures up to 230 °C were estimated by geothermometric calculations based on the water chemistry of the thermal discharges [29]. Geological and geochemical studies were conducted but no wells were drilled.

Geothermal exploration in the Central Andean Volcanic Zone was abandoned in 1982 because of both the remote location of the hydrothermal systems and economic factors. After almost three decades, private and governmental companies have planned to conduct a new phase of geothermal exploration in the systems investigated between 1969 and 1982 as well as in other areas of northern Chile. New slim holes will be drilled at Puchuldiza, Polloquere, Pampa Lirima, Colpitas and Juncalito by Energía Andina, which plans to have a geothermal plant working by 2015 [30].

In central-southern Chile, the geothermal activity is related to the Pocuro fault system (33°–34° S) and to the Liquiñe-Ofqui Fault Zone (39°–46° S). Several slimholes have been drilled in Tinguiririca, Calabozos, Laguna del Maule, Chillán and Tolhuaca with a potential output estimated at 3–10 MW per well (Figure 4). At Tolhuaca, two holes have been drilled up to depths of 1200 m and a 50 MW geothermal power plant is planned to start production in 2013. The potential is estimated at 600 MWe to 950 MWe in this area [31]. Currently, 79 exploration and 7 exploitation concessions have been given to Chilean and foreign companies [3].

**Figure 4.** Geothermal areas of southern Chile.

#### **3. Methodology**

In the selected potential geothermal fields, the minimum and maximum electric energy supply (*Esu*) was evaluated considering an operation time (*OT*) of 8000 h/year. To evaluate the *Esu* it was necessary to estimate both heat (*Q*) and electric power (*We*).

The recoverable *Q* of examined geothermal reservoir was computed using the following equation:

$$
\underline{Q} = \underline{m} \cdot C\_{\overline{n}} \cdot (T - T\mathfrak{d}) \tag{1}
$$

where *m* is the mass flow rate of the geothermal system (kg/s), *Cw* (J/kg·K) is the specific heat capacity of the fluids contained in the geothermal reservoir, *T* is the reservoir temperature (K) and *T*0 is the reinjection temperature. The minimum reservoir temperature *T*í and the maximum reservoir temperature *T*+ were used to compute the minimum (*Q*í) and maximum (*Q*+) heat, respectively.

The production rate *m* of the examined geothermal system is obtained by multiplying the volume of the geothermal reservoir *V* (km3 ) and the specific productivity (*mw*) [32] derived from the flow well tests. The *We* was estimated using the following equation:

$$W\_{\mathfrak{e}} = \underline{Q} \cdot \mathfrak{n} \tag{2}$$