Analysis of the magmatic – Hydrothermal volcanic field of Tacora Volcano, northern Chile using travel time tomography

Journal of South American Earth Sciences 94 (2019) 102247

Contents lists available at ScienceDirect

Journal of South American Earth Sciences journal homepage: www.elsevier.com/locate/jsames

Analysis of the magmatic – Hydrothermal volcanic field of Tacora Volcano, northern Chile using travel time tomography

T

C. Paveza,b, D. Comteb,c,∗, F. Gutiérrezd, D. Gaytáne a

Departamento de Geología, Facultad de Ciencias Físicas y Matemáticas, Universidad de Chile. Plaza Ercilla 803, 13518, Santiago, Chile Advanced Mining Technology Center, Facultad de Ciencias Físicas y Matemáticas, Universidad de Chile. Avenida Tupper 2007, 8370451, Santiago, Chile c Departamento de Geofísica, Facultad de Ciencias Físicas y Matemáticas, Universidad de Chile. Avenida Blanco Encalada 2002, 8370449, Santiago, Chile d GeoExpedition, Las Palomas 25, Pirque, Santiago, Chile e INFINERGEO SpA, Chile b

A R T I C LE I N FO

A B S T R A C T

Keywords: Tacora volcano Passive seismic tomography Conceptual model Clay levels Clustering method

Tacora Volcano (17°43′S – 69°46′W) lies at the southernmost end of a 10 km-long volcanic lineament that extends between Chile and Perú. Around Tacora volcano, thermal manifestations are two active fumarolic fields located at the western flank of the stratovolcano and at the volcano summit, indicating active magma degassing in a shallow hydrothermal system. Beneath Tacora volcano is located the NW Challaviento reverse fault that belongs to the Incapuquio - Challaviento fault system of Middle Eocene age. To complement previous exploration results and conceptual modeling developed by INFINERGEO SPA, seventeen short period seismic stations were installed around Tacora Volcano, between August and December 2014. Using the P and S wave arrival times of locally recorded seismicity, a 3D velocity model was determined through a travel time tomography. We interpreted high Vp/Vs values as water-saturated areas, corresponding to the recharge zone of Tacora hydrothermal system. In addition, low values of ΔVp/Vp (%) and Vp/Vs ratio represent the location of a gas-saturated magmatic reservoir and circulation networks of magmatic-hydrothermal fluids. Low Vp/Vs volumes (magma reservoirs/high temperature hydrothermal fluids), the presence of fumarolic fields and surface hydrothermal alteration have a spatial correlation. This suggests a structural control of the Challaviento fault in the hydrothermal flow. Finally, we present a cluster analysis using the ΔVp/Vp (%) parameter. Through this analysis, we found a method for the identification of a key structure in depth composed by the magma reservoir (low Vp/Vs ratios, low ΔVp/Vp (%)), clay level areas (intermediate values of ΔVp/Vp (%)), and degassification zones (low values of ΔVp/Vp (%)) directly related with the surface thermal manifestations.

1. Introduction Determination of the internal structure of magmatic – hydrothermal systems using travel time seismic tomography (Koulakov, 2013; Gritto and Jarpe, 2014; Saccorotti et al., 2014; Pavez et al., 2016) is a low cost method widely used to study active volcanic areas (Lees, 2007). Using this method, it is possible to identify and map geological structures (Jousset et al., 2011), infer the hydrodynamic state of the area (Tryggvason et al., 2002) and locate spatially the magmatic reservoirs (Pavez et al., 2016). This method enables the generation of geometrical constraints for conceptual models that explains how the geothermal system works (Kuhn, 2004). Furthermore, different petro-physical parameters of the host rock such as temperature, composition and density, can be linked to the 3D velocity structure of P- and S- seismic

waves, Vp/Vs ratios and the percentage variation parameter ΔVp,s/Vp,s (%) (Lees, 2007). These variables could also indicate the presence of magmatic reservoirs, and distinguish between fluid and melt zones in active volcanoes, which is also used as a geothermal exploration tool in volcanic – hydrothermal systems (Koulakov, 2013). Particularly, ΔVp/ Vp (%) has been widely used to reflect volcanic features, such as volcanic deposits, magmatic intrusions, geological structures and lithological domains (Tanaka et al., 2002; Muksin et al., 2013). Its analysis has also been related with temperature changes and fluid density characterization (Poletto et al., 2018). Changes in the seismic velocities through the medium are commonly associated with lithological differences (i.e. magmatic intrusions, volcanic deposits, sedimentary sequences), or even the presence of fluids at different temperatures (Lees, 2007). For instance, in volcanic systems, high Vp/Vs (1.75–1.80)

∗ Corresponding author. Departamento de Geofísica, Facultad de Ciencias Físicas y Matemáticas, Universidad de Chile. Avenida Blanco Encalada 2002, 8370449, Santiago, Chile. E-mail address: [email protected] (D. Comte).

https://doi.org/10.1016/j.jsames.2019.102247 Received 27 January 2019; Received in revised form 17 June 2019; Accepted 23 June 2019 Available online 24 June 2019 0895-9811/ © 2019 Elsevier Ltd. All rights reserved.

Journal of South American Earth Sciences 94 (2019) 102247

C. Pavez, et al.

the Tacora magmatic-hydrothermal system the descriptive terms proposed by Nicholson (1993) and suggested in the model of Capaccioni et al. (2011). These correspond to: reservoir equilibrium state, fluid type, reservoir temperature, relieve and heat source. This 3D model allows to explain the current distribution of the physical and chemical surface manifestations observed in the hydrothermal system (Capaccioni et al., 2011; Gaytán, 2012, 2013; Contreras, 2013). Additionally, it provides additional insights into the subsurface geometry and the hydrodynamic of the system related to existing faults. So far, this work corresponds to the unique research completely oriented to observe and understand the behavior of Tacora volcano in depth, using the complement between geological and geophysical results.

velocity anomalies with high Vp values observed in shallow crustal regions represent structural features associated in general with rocks containing fluids (De Natale et al., 2004; Lees, 2007). On the other hand, low Vp/Vs (1.65–1.70) anomalies are interpreted as evidence of melt accumulation, like magma reservoirs and active conduits of liquids and magma (Londoño and Sudo, 2003; Pavez et al., 2016). Indeed, partial melting is generally related with a decrease in Vp (Nakajima et al., 2001; Foulger et al., 2003; De Matteis et al., 2008). However, fluid saturated zones also show low Vp and Vs, producing a high Vp/Vs ratio (Walck, 1988; Lees, 2007). In hydrothermal systems, high Vp/Vs velocity anomalies are ussually associated to steam condensation zones, areas where meteoric water infiltrates, saturation conditions and zones where hydrothermal gases migrate from the reservoir to distal areas due to increase of effective pressure (Ito et al., 1979; Simiyu, 1999; Lees and Wu, 2000; Foulger et al., 2003; De Matteis et al., 2008). In contrast, low Vp/Vs anomalies associated to low Vp values are related with high temperature zones of magmatic volatiles, steam saturated formations and hydrothermally-altered zones (Foulger and Miller, 1995; Boitnott and Kirkpatrick, 1997; Simiyu, 1999; Foulger et al., 2003; Yoshikawa and Sudo, 2004; Vanorio et al., 2005; De Matteis et al., 2008). Particularly, the joint interpretation between ΔVp/Vp (%) and Vp/Vs is widely used to study geothermal areas (Muksin et al., 2013). The combination of these parameters improves the elucidation of the seismic velocity models in terms of lithology and fluid content. The southern segment of the Central Volcanic Zone, in northern Chile, has been the focus of increasing attention due to the large number of hydrothermal systems. In this zone, there are about 90 geothermal areas, of which only a few are under systematic exploration (Tassi et al., 2005, 2010; Lahsen et al., 2010; Aravena et al., 2016). Particularly, Tacora volcano (17°43′S – 69°46′W) (Fig. 1) is under study due to the presence of thermal manifestations, like hot springs and fumarolic fields, and the presence of extensive areas with hydrothermal alteration, which makes it a promising geothermal energy field (García, 2011; Gaytán, 2012, 2013; Contreras, 2013). Through a travel time seismic tomography we attempt to better characterize the shallow hydrothermal system hosted by Tacora volcano, supporting and improving the model proposed by Capaccioni et al. (2011). The information derived from the seismicity distribution (i.e. ΔVp/Vp (%) and Vp/Vs structure) is used to locate in a 3D model of

2. Geological setting Tacora Volcano is a composite stratovolcano located at ∼100 km northwest of Arica, northern Chile region (Figs. 1 and 2). It has a summit elevation of 5980 m a.s.l, and 1700 m from its base (Barrientos, 2013). It lies at the southernmost end of a 10 km-long N-S oriented volcanic lineament that comprise Nevado La Monja, Nevado El Fraile (Peru), and Chupiquiña volcanoes (Fig. 2b). An Andesite sample from the base of Tacora volcano was dated using the K-Ar method showing an age of 0.489 ± 0.015 Ma (Wörner et al., 2000). Meanwhile two lava samples from the volcanic edifice were dated by Barrientos (2013) through the 40Ar/39Ar StepHeating method. These samples exhibit ages of 0.340 ± 0.060 and 0.363 ± 0.007 Ma, all corresponding to the Middle Pleistocene. Two main sub-units have been identified in the volcano area (Contreras, 2013). The first one is composed by andesitic to dacitic lava flows, which form a large part of the volcanic building, and pyroclastic rocks (Montecinos, 1970; Barrientos, 2013). Meanwhile, the second sub-unit is formed by a dacitic dome, exposed in the middle north sector of the volcano whose rocks are composed by phenocrystals of plagioclase, sanidine, biotite, amphibole and quartz in smaller amounts. This sub-unit is slightly altered to clay minerals (Barrientos, 2013; Contreras, 2013). Moraines composed by volcanic fragments with thicknesses up to 200 m (Contreras, 2013) and glacial valleys related to the Holocene glaciation (∼11,000 ka) are recognized in the E, S and SE flanks and at the base of the volcano (De Silva and Francis, 1991; Capaccioni et al., 2011; Contreras, 2013).

Fig. 1. Digital elevation model showing the regional context of Tacora volcano. Inside the red box Tacora volcano and the volcanic lineament are shown. To the lower left corner, the inset is showing the digital elevation model corresponding to an overview of South America. This overview is also displaying the limits for the Northern, Central and Southern Volcanic Zones along the Central Andes. A yellow box indicates the outline corresponding to the regional context of this study. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) 2

Journal of South American Earth Sciences 94 (2019) 102247

C. Pavez, et al.

Fig. 2. (a) Structural geology in the southern border region of Perú. Incapuquio – Challaviento fault system and its inferred extension to northern Chile is highlighted (modified from Contreras, 2013). The study area is indicated with a dotted yellow box. To the lower right corner, a zoom out is showing in orange the complete extension and shape of La Yarada aquifer (Rojas, 2008) (b) Local map of Tacora volcano and the surrounding area showing the seismic network coverage. Green triangles indicate the location of the temporary seismic stations used in this work. Purple squares correspond to the NW and W fumarolic fields registered by Capaccioni et al. (2011), and to the weak summit degassing reported by Contreras (2013). Pink ascending arrows indicate the up-flow zone proposed by Gaytán (2012, 2013). To the lower left corner, UTM zone number and designator are shown in white. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Barrientos (2013) describes superficial pyroclastic flow deposits stratigraphically above Pliocene to Pleistocene volcaniclastic rocks, that could not be dated due to atmospheric related weathering and

According to Capaccioni et al. (2011), the last explosive activity of the Tacora volcanic system is probably related with an explosion crater about 500 m in diameter, located on the NW side of the volcano. 3

Journal of South American Earth Sciences 94 (2019) 102247

C. Pavez, et al.

(Capaccioni et al., 2011; Contreras, 2013). Meanwhile hot springs (Aguas Calientes de Tacora) are present at ∼2 km to the SW of the volcano (Aguilera, 2010). These fluids have temperatures near 40 °C and high content of SO4 in addition to Cl, Na, Ca and Mg. Based on the aforementioned manifestations, several areas with hydrothermal alterations have been identified. These are particularly located on the volcanic building, and show argillic alteration and silicification, with temperatures bordering 100 °C and pH in the range of 1–2 (Contreras, 2013). The study area is close to the southernmost border of La Yarada shallow aquifer (Rojas, 2008), located in Tacna, Perú (Fig. 2a). This aquifer drains toward the Pacific Ocean, and its waters are originated by rainfall in the western part of the Andes (Rojas, 2008). Although the vertical limits of the aquifer are not known exactly, the depth ranges from 100 m to 1 km, reaching its maximum depth near the Chile – Perú border (Rojas, 2008). Beneath Tacora volcano is located the Challaviento reverse fault that belongs to the Incapuquio - Challaviento fault system of Middle Eocene age (Fig. 2a). The Incapuquio – Challaviento fault system (ICFS) is exposed northwestward through Peruvian territory. The southeastern prolongation in Chile (Fig. 2a) is covered by younger volcaniclastic units (Oaxaya and Visviri Formations; Oligocene - Pleistocene) and the previously mentioned volcanic units and deposits (Contreras, 2013; David et al., 2004). This system, of WNW-ESE orientation, is mainly composed by two reverse faults, Incapuquio and Challaviento, and other secondary structures (David et al., 2004). The ICFS includes also subvertical faults and currently is seismically active in Peruvian territory (David et al., 2004; Flores et al., 2005). So far, no research was done regarding its seismic activity in Chile.

Fig. 3. L-curve for the inversion process. (R2: 0.995, RMSE: 2.976).

Table 1 Percentage variation per iteration. Iteration

Percentage Improvement

1 2 3 4 5 6

88.42% 44.83% 17.83% 13.63% 6.33% 3.91%

3. Methodology alteration. Currently, the caldera is covered by glaciers above ∼5500 m. Hantke (1939) describes two eruptive events in 1930 and 1937, although there is no recognizable record of these activities. Presently, permanent fumarolic degassing characterizes the upper portion of the NW and W flanks and weak degassing occurs on the volcano summit (Contreras, 2013) (Fig. 2b). These manifestations, developed mainly above the glacial deposits, are mainly composed by vapor water, CO2, N2, H2S, CH4 and magmatic gases SO2, HF and HCl

3.1. Travel time seismic tomography To achieve an accurate characterization of the magmatic - hydrothermal system, a 3D body wave velocity model was determined. A seismic network of 17 short period, three-components, and continuously recording seismic stations was deployed in the area of Tacora volcano from August to December 2014 (Fig. 2b). The coverage area was 20 × 15 km2 including both, the volcanic edifice and the NE

Fig. 4. (a) Local map of Tacora volcano and the surrounding area. Nevado La Monja and Chupiquiña volcanoes has been marked with a light blue arrow. Locally recorded shallow seismicity is observable from the surface. The inferred Challaviento fault has been drawn in the area according to Contreras (2013). Continuous lines correspond to the vertical sections shown in Fig. 5b and d (b) Final hypocenters of local seismicity beneath Tacora volcano are shown in a 3D view. A shallow cluster is boxed in both figures. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) 4

Journal of South American Earth Sciences 94 (2019) 102247

C. Pavez, et al.

Fig. 5. Iso-surfaces with 3D projection showing the Vp/Vs ratios of the final velocity model of Tacora volcano and the surrounding area. All sections use the same color palette for velocity ratios. (a) Low Vp/Vs rates in the range of 1.65–1.68. (b) Low Vp/Vs rates: north – south cross section showing the velocity ratio anomalies in the range of 1.65–1.68. Profile location is shown in Fig. 4a (c) High Vp/Vs rates in the range of 1.74–1.80. (d) High Vp/Vs rates: north - south cross section showing the velocity ratio anomalies in the range of 1.74–1.80. Profile location is shown in Fig. 4a. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 6. Cross section with a 2D projection of the final velocity model for Vp/Vs. An inferred fault is shown with a dotted line. At the surface, a fumarolic field is represented with a light orange symbol. Blue and red zones are showing high and low Vp/Vs values (under 1.68 and above to 1.76). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

P wave arrival times are between 0.1 and 0.5 s while those for S waves typically are about twice that amount. The hypocenters were preliminarily determined using HIPOINVERSE (Klein, 1978) based on a 1D velocity model for P- waves (Thierer et al., 2005), and a Vp/Vs ratio of 1.75. The preliminarily hypocenter locations were used as starting values for the tomography inversion. The inversion codes correspond to SPHYPIT90/SPHREL3D90, which

geothermal system area identified by Gaytán (2012, 2013). The body wave data set consisted of P- and S- arrival times from about 1000 events. Earthquakes that were recorded by at least 10 stations and with predicted arrival times using the 1D model within an initial outlier residual threshold of larger of 2 s and 10 per cent of the total travel time were selected. Application of these criteria resulted in a reference data set of 7413 P- and 7269 S- wave arrival times. Nominal uncertainties for 5

Journal of South American Earth Sciences 94 (2019) 102247

C. Pavez, et al.

Fig. 7. Ray density profiles with elevation crossing the study area. To the left, the location of the profiles is shown. The image was generated using Google Earth Pro®, with map data from DigitalGlobe. To the right, AA′ profile, with origin and angle. Green triangles displays the projection of the seismic stations. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) Fig. 8. Cross section with 2D projections of the final velocity model for ΔVp/Vp (%). The proposed clay level and degasification zones are in dotted line. Above them, the fumarolic field is showed with a light orange symbol. The Challaviento fault is displayed with a dotted line. Blue and red zones are showing high and low ΔVp/Vp (%) variations (between 10% − 15% and −10% − −15%, respectively). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

inversion to be iterative, thus possible biases related with the nonlinearity and the reference starting model are reduced. First, all the earthquakes in the original data set are relocated in a given structure, that later is perturbed to minimize the residuals. These steps are performed iteratively until has no significant reduction in the residual variance (Roecker et al., 1987). In this work, the 3D body wave velocity model was obtained in a grid of 594 double blocks, covering a volume of 20 × 15 × 24.5 km3 that includes the volcanic building and the surrounding area. The region was laterally parameterized by 3 × 3 km2 grid spacing, and 1–1.5 km depth separation along 11 layers, ranging from 4.5 km above

were developed by Steve Roecker (Roecker, 1982; Roecker et al., 1987, 1993). The three-dimensional velocity model is determined through minimizing travel time residuals of P- and S- waves from locally recorded earthquakes. The three dimensional structure is parameterized by a set of layers with horizontal interfaces, each of which is divided into a regular orthogonal mesh. The P- and S- wave velocities in each block of the grid are specified independently, and perturbations and adjustments to these velocities and hypocentral coordinates are determined simultaneously (Roecker et al., 1987). Contributions to travel times from station elevations are calculated directly including station corrections. The ray tracing technique allows to the three dimensional 6

Journal of South American Earth Sciences 94 (2019) 102247

C. Pavez, et al.

Fig. 9. (a) Clustering of Vp/Vs ΔVp/Vp (%) using k = 3 and showing the interpretation in the structure proposed in Fig. 8 (b) Silhouette technique for clusters 1, 2 and 3. All the clusters reach a high value (> 0.6).

sea level to −20 km below. The iteration process stopped when the maximum change per iteration and the total residuals were less than 5% and 2%, respectively. After 6 iterations, we obtained 4% of variation of the model with respect to the previous iteration. For the last iteration, 742 blocks were used during the inversion process. Fig. 3 and Table 1 shows the percentage of improvement per iteration. In Fig. 4, the final hypocenters location for shallow seismicity in the vicinity of Tacora volcano can be observed.

4. Final velocity model Based on the thermal manifestations of the study area (Capaccioni et al., 2011, Fig. 2b) and the shallow seismicity distribution (Fig. 4), four 3D iso-surfaces for the Vp/Vs ratio have been generated (Leapfrog Viewer software®, v.4.4.0, 2014). Iso-surfaces 5a & 5c image the low and high Vp/Vs values beneath the study area. These models were cut into vertical sections to better observe the geometry and orientation of the Vp/Vs anomalies (Fig. 5b and d). Profile lines can be observed in Fig. 4a. In the eastern part of the study area, a low velocity anomaly, with ratios under 1.68 (Fig. 5a and b) predominates. The geometry of this low velocity anomaly (Fig. 5b) suggests the presence of a high angle structure (∼60°) with north vergency and NE strike. Meanwhile, going to the west and directly below Tacora volcano, high Vp/Vs values have been identified (Fig. 5c and d). These values, above 1.76, predominate to a depth between ∼4 km a.s.l and 4 km b.s.l, extending from north to south and even beneath the nearby volcanic lineament (Fig. 4). A vertical cross section directly beneath the NW fumarolic field can be observed in Fig. 6. This cross section shows several variations in the Vp/Vs value. The profile shows a low Vp/Vs anomaly between the sea level and 2 km below, that exhibits the same geometry and vergence that the low ratio anomaly observed in Fig. 5b, to which overlies a small intermediate Vp/Vs anomaly. A high Vp/Vs zone can also be observed to the south of the study area. To complement this model, more cross sections can be found in Supplementary Material 1. With the aim to observe the well resolved areas in the velocity model, a ray density profile -AA′- was generated (Fig. 7). In this profile, N(P + S) indicates the total sum between P wave and S wave rays used in the final inversion. The profile has a width of ± 1 km and is covering the entire area in EW direction. AA′ shows a high ray density directly under Tacora volcano and until −4 km depth (N(P + S) ≥400). This number (at least 400 P and S rays per element of the grid), supports a good resolution for the entire area, particularly under the volcanic

3.2. K-means technique for cluster analysis K-means clustering corresponds to a partitioning method which creates k mutually exclusive clusters. The number of clusters k, cluster initialization and distance metric are parameters specified by the user. Different initializations can finish in a different final clustering because K-means only converges to a local minimum. To solve this, the K-means algorithm is runned with multiple different initial partitions, choosing the partition with the smallest squared error. K-means typically use Euclidean distance to compute the distance between points and cluster centers, obtaining spherical or ball-shaped clusters. The centroid for each cluster is defined by the point to which the sum of distances from all objects in that cluster is minimized. The summation process is iterative, until the minimum value is reached (Jain, 2010). To measure the validity and consistency of the clusterization, a Silouehette method was used. The Silouhette values indicate how similar is an object to its own cluster, and how diferent from the others. The values range between −1 and 1. If the majority of the objects have high and positive values, then the configuration is valid (Rousseeuw, 1987). The cluster modeling was determined using the k-means internal function available in Matlab. Vp/Vs ratio and ΔVp/Vp (%) were clustered. An Euclidean distance was selected, in order to be in agreement with the spatial location of the parameters. The number of partitions was k = 3, and the minimum value for the total sum of distances was reached after 5 iterations. 7

Journal of South American Earth Sciences 94 (2019) 102247

C. Pavez, et al.

Fig. 10. Checkerboard test results for Vp and Vp/Vs parameters. a) & b) are showing the original checkerboard pattern, c) & d) are showing the recovered models in EW and NS direction for Vp. Finally, e) & f) correspond to the recovered models in EW and NS direction for Vp/Vs. The same color palette is used for both recovered models. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

and the contact zones between low and high velocity ratios (Fig. 6), we infer that this structure corresponds to the Challaviento fault, that belongs to the Incapuquio – Challaviento fault system (Figs. 2a and 4a; Acosta et al., 2011). The ICFS exerts an important structural control in the southern zone of Perú, and it is seismically active to the Peruvian territory (David et al., 2004; Flores et al., 2005). However, this behavior is not reflected in the proximity of Tacora volcano, and there is no observable relation when the local seismicity is projected to the surface (Fig. 4a). Regarding this issue, some of the local seismic activity is likely related to brittle fracturing generated by fluid mobility, or with the emplacement and location of a magmatic reservoir associated with Tacora volcano. Both, seismicity cluster (Fig. 4a and b) and the proposed magmatic chamber are spatially correlated (Fig. 6). To differentiate and better understand some of the previous features observed in the Vp/Vs model, a ΔVp/Vp (%) 2D cross section was developed (Fig. 8). This final velocity model indicates that the

building.

5. Discussion and conceptual model High Vp/Vs values (1.80–1.85) located at the western sector of the area are interpreted as fluid-saturated rocks (Fig. 5c and d). The origin of these fluids is attributed to glaciers situated above ∼5000 m a.s.l, local precipitation and the proximity of a shallow aquifer of andesitic water (La Yarada) (Taran et al., 1989; Giggenbach, 1992; Rojas, 2008; Capaccioni et al., 2011). Above these high Vp/Vs values, the low temperatures recorded in the W fumarolic field - with temperatures near to the water boiling point (82° - 93 °C) (Capaccioni et al., 2011); -corroborate the presence of liquid water. Orientation and vergence of the main eastern anomalies suggest the presence of a high angle and locally blind structure below the volcanic cover (Figs. 5b and 6). Using the geometry of the low Vp/Vs anomalies 8

Journal of South American Earth Sciences 94 (2019) 102247

C. Pavez, et al.

Fig. 11. Final velocity model at surface level (> 5000 m a.s.l). Section set to the local topography. The inferred projection of Challaviento fault is shown with dotted lines. The fluid migration is indicated with light blue arrows. Purple squares show the thermal manifestations described and located in the area by Aguilera (2010) & Capaccioni et al. (2011). Blue and red zones indicate high and low Vp/Vs rates (between 1.76 and 1.80 and 1.65–1.68, respectively). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Challaviento fault represents the path of magmatic hydrothermal fluids from depth to surface, similar to it was observed in other volcanic and geothermal systems (Foulger, 1982; Muksin et al., 2013; Pavez, 2016; Pavez et al., 2016). This assumption is supported by:

In the same way, this clay level matches at the surface with the fumarolic field registered by Capaccioni et al. (2011). The existence of clay levels is also strongly supported at the surface by the presence of hydrothermal alteration zones of argillic and advanced argillic type in the volcanic edifice and the surrounding area (García, 2011; Contreras, 2013). The forming conditions of these alterations corroborate acid hydrothermal fluids near to 100 °C, which is coincident with the recorded temperatures in the NW fumarolic field (Capaccioni et al., 2011). Combined with the presence of fumarolic fields, surface alterations also support the existence of a shallow magmatic-hydrothermal reservoir (Figs. 2b, 4a and 6), because they form from the interaction between sulfated fluids with meteoric shallow waters having magmatic source compounds (Nicholson, 1993; Contreras, 2013). Distinct rock types or fluid content can be identified by the application of cross-plots and statistical cluster analysis techniques (Muksin et al., 2013). To correlate empirically velocity ratio ranges of seismic waves with potential lithological domains, a cluster analysis was then carried out. The clustering was done using the ΔVp/Vp (%) parameter, mainly due to the identifiable changes in the P- wave velocity. The clearly identifiable variations in ΔVp/Vp (%) (Fig. 8), represents stronger variations of κ (bulk modulus) over μ (shear modulus), so gas or liquid phases can be identified. Based on this assumption, a k-means clustering function was used (Fig. 9a). The number of clusters was selected by analyzing the Silhouette value, which results in k = 3 (Fig. 9b). In the first cluster, κ presents a percentage decline, indicating high temperatures and liquid or gaseous phases. Meanwhile, the third cluster (Fig. 9) groups a percentage increase suggesting low temperatures and solid phases. Finally, the second cluster, with a data average of m = 1.712, could be related to the location of clay caps, if and when it is followed by a decrease in the ΔVp/Vp (%) value, that is linked with a gas phase (Fig. 8; O'Connell and Budiansky, 1974; Toksoz et al., 1976). In our case, the decreasing in the ΔVp/Vp (%) value is observed (Fig. 8) and is directly related with a degassing zone, -the NW fumarolic field-. Therefore, the structure would allow us to identify clay levels inside the magmatic-hydrothermal systems is the next one, from depth

1. The inferred fault (Challaviento) runs through the main hydrothermal manifestation located at the NW flank of the Tacora volcano (Capaccioni et al., 2011) (Figs. 4a, 6 and 8). 2. The small seismicity cluster, recorded between 2 km depth and sea level (Fig. 4b), is also coinciding with a low value of ΔVp/Vp (%) (Fig. 8). The location of this cluster, observed below the E flank of the Tacora volcano (Fig. 4a and b), is in addition spatially correlated with the fluid up-flow zone proposed by Gaytán (2012, 2013; Fig. 2b). This could indicate that the cluster is related to the circulation of a large amount of hot fluids, which in turn are producing the brittle fracturing in the area. The location and orientation of the low Vp/Vs anomaly (Vp/ Vs < 1.68) and low ΔVp/Vp (%) next to the Challaviento fault (Fig. 5a and b, 6 & 8) suggests the presence of a magmatic-hydrothermal reservoir that governs the temperature range (270–310 °C) of the geothermal reservoir proposed by Capaccioni et al. (2011). In addition to the proposed reservoir, the orange zones observed in Fig. 5a and b could correspond to a network of active channels. These channels responsible for the circulation of magmatic – hydrothermal fluids - can be supported considering the link between anomalies, surface manifestations, and the interpretative fluid circulation model presented by Capaccioni et al. (2011). The presence of this magmatic-hydrothermal reservoir gives rise to clay levels (Gaytán, 2012, 2013, Contreras, 2013). These are identified with intermediate Vp/Vs values (Fig. 6), low ΔVp/Vp (%) values (Fig. 8) and low resistivity (Soyer and Hallinan, 2012), due to the interaction of soluble magmatic fluids (HCl, HF and SO2) that condense in the shallow aquifer. Under this analysis, the NW area described as a clay level by Capaccioni et al. (2011), is now supported and spatially located by the seismic tomography (Figs. 6 and 8). 9

Journal of South American Earth Sciences 94 (2019) 102247

C. Pavez, et al.

Fig. 12. Schematic representation of the conceptual model of the magmatic-hydrothermal system hosted by Tacora Volcano (modified from Capaccioni et al., 2011; Contreras, 2013), based on geochemistry of fluid circulation and seismic tomography. The results presented by Capaccioni et al. (2011) are marked with two asterisks (**). The results obtained in this paper are shown with one asterisk (*). Common interpretations are marked with both symbols. As in most igneous systems in active volcanoes, the magmatic-hydrothermal system can be represented with different Vp/Vs values: a deep magmatic reservoir with the lowest Vp/Vs values; a hydrothermal reservoir with hot magmatic volatiles (low Vp/Vs); and a shallow aquifer with highest Vp/Vs values at the upper part of the system, given by the interaction between meteoric waters and magma volatiles. Clay minerals, formed in relation with the magmatic process, appears beneath the fumarolic field (intermediate Vp/Vs).

to surface: a low ΔVp/Vp (%) value, linked with the hydrothermalmagmatic reservoir, an intermediate ΔVp/Vp (%) value linked with the clay level, followed by a low and well located low ΔVp/Vp (%) related with a degasification zone and fumarolic field type thermal manifestations on the surface (Figs. 8 and 9). To check the reliability of our model and in order to investigate which kind of anomalies can be reconstructed, we developed a checkerboard test. For this, the models for Vp and Vp/Vs were perturbed by anomalies alternating with ± 5% of variation, creating a checkerboard pattern. NS and EW profiles were selected in order to show the results. Fig. 10a and b shows the original pattern. Meanwhile, Fig. 10c and d and 10e & 10f exhibit the recovered checkerboard models for Vp and Vp/Vs, respectively. The reconstructed images for Vp indicate a good resolution until −3.5 km depth at local scale for both NS and EW profiles. This is supporting the interpretation of the model mainly under Tacora volcano, reinforcing the magmatic chamber location and the degasification zone related to the NW fumarolic field. This high-resolution area also coincides with the obtained using the ray density analysis (Fig. 7). On the other hand, the anomalies for the Vp/Vs model are well resolved until −1.0 km b.s.l, and with degradation in recovery quality under this depth. However, considering the obtained results for the Vp checkerboard test and the ray coverage, this should not affect our interpretation below that depth. It is also relevant to point out that our interpretation of the results obtained through the seismic tomography includes only features that have sizes of the same order of

magnitudes and locations of the recovered models. To generate a spatially constrained analysis of the interpretative model of fluid circulation proposed by Capaccioni et al. (2011), a surface horizontal projection of the magmatic-hydrothermal system was developed (Fig. 11). This model consists of an overlap between topography, fumarolic fields location, and the first and shallow layer corresponding to the Vp/Vs velocity model (between 3.0 and 4.5 km a.s.l). High Vp/Vs values indicate the location of the recharge zone of the magmatic – hydrothermal system. Topography of the area, descending eastward (from ∼5.500 to ∼4.000 m.a.s.l) favors fluid migration, which results in an area with intermediate Vp/Vs values (1.68–1.74). These fluids are part of the hydrodynamic activity of the system generating convection cells around low velocity anomalies (Vp/Vs < 1.68). Regarding the surface alterations, the low value anomaly located to the north-east (Fig. 11), would be related to the argillic alteration previously observed in the zone (Contreras, 2013). Through the final models obtained for Vp/Vs and ΔVp/Vp (%) (Figs. 8 and 11) an schematic model that enhance and support the results presented by Capaccioni et al. (2011) was generated (Fig. 12). With them, is possible to locate spatially the fault system under the volcanic building, the magmatic source, a NW clay level and the source of gases proposed by Capaccioni et al. (2011). Finally, both models allowed us to classify the system according to the parameters proposed by Nicholson (1993). The magmatic – hydrothermal system hosted in Tacora volcano, would correspond to a high enthalpy dynamic system, 10

Journal of South American Earth Sciences 94 (2019) 102247

C. Pavez, et al.

of andesitic – magmatic fluid type, vapor dominated and high relief, which lays within category of many others Andean magmatic complexes (Aravena et al., 2016).

References Acosta, H., Alván, A., Mamani, M., Oviedo, M., 2011. Geología de los Cuadrángulos de Pachía y Palca. Hojas 36-v y 36-x. Bol.139 Carta. Geol.Nac. INGEMMET, Lima. Aguilera, F., 2010. Origin and Evolution of Fluid in Volcanoes, Geothermal Fields and Thermal Discharges of Central Volcanic Zone in Northern Chile (17°43’S-25°10’S). Bicentury Ph. D Thesis. 2008. Aravena, D., Muñoz, M., Morata, D., Lahsen, A., Parada, M.A., Dobson, P., 2016. Assessment of high enthalpy geothermal resources and promising areas of Chile. Geothermics 59, 1–13. Barrientos, J., 2013. Evaluación y zonificación preliminar del peligro volcánico del Volcán Tacora, XV región de Arica y Parinacota, Andes Centrales del norte de Chile. Geologist Thesis. Universidad de Chile. Boitnott, G., Kirkpatrick, A., 1997. Interpretation of field seismic tomography at the geysers geothermal field, California. In: 22th Workshop on Geothermal Reservoir Engineering, Standford University, Stanford, California, January 27-29, 1997. Capaccioni, B., Aguilera, F., Tassi, F., Darrah, T., Poreda, R.J., Vaselli, O., 2011. Geochemical and isotopic evidences of magmatic inputs in the hydrothermal reservoir feeding the fumarolic discharges of Tacora volcano (northern Chile). J. Volcanol. Geotherm. Res. 208, 77–85. Contreras, A., 2013. Caracterización de la mineralogía de alteración hidrotermal en superficie del volcán Tacora y sus alrededores, Región de Arica y Parinacota. Geology Thesis. Universidad de Chile. David, A., Audin, L., Tavera, H., Hérail, G., 2004. Sismicidad cortical y fallas recientes en el sur del Perú. Resúmenes extendidos del XII Congreso Peruano de Geología. Public. Esp. N°6 290–293 SGP. De Matteis, R., Vanorio, T., Zollo, A., Ciuffi, S., Fiordelisi, A., Spinelli, E., 2008. Threedimensional tomography and rock properties of the Larderello-Travale geothermal area. Italy. Phys. of the Earth and Planetary Interiors 168, 37–48. https://doi.org/10. 1016/j.pepi.2008.04.019. De Natale, G., Troise, C., Trigila, R., Dolfi, D., Chiarabba, C., 2004. Seismicity and 3D substructure at Somma- Vesuvius Volcano: evidence for magma quenching. Earth Planet. Sci. Lett. 221, 181–196. De Silva, S., Francis, P., 1991. Volcanoes of Central Andes. Springer Verlag, Berlin, pp. 216. Flores, A., Acosta, J., Bedoya, C., Sempere, T., 2005. Oligocene – neogene tectonics and sedimentation in the forearc of southern Peru, Tacna area (17. 5° - 18. 5° S). In: 6th International Symposium of Andean Geodynamics. Extended Abstract, pp. 269–272. Foulger, G., 1982. Geothermal exploration and reservoir monitoring using earthquakes and the passive seismic method. Geothermics 11 (4), 259–268. Foulger, G.R., Miller, A.D., 1995. Three-dimensional vp and vp/vs structure of the Hengill Triple Junction and geothermal area, Iceland, and the repeatability of tomographic inversion. Geophys. Res. Lett. 22 (10), 1309–1312. https://doi.org/10.1029/ 94GL03387. Foulger, G.R., Julian, B.R., Pitt, A.M., Hill, D.P., Malin, P.E., Shalev, E., 2003. Threedimensional Crustal Structure of Long Valley Caldera , California , and Evidence for the Migration of CO 2 under Mammoth Mountain, vol. 108. pp. 1–16. 564. http:// doi.org/10.1029/2000JB000041. García, M., 2011. Reconocimiento geológico preliminar del área del Volcán Tacora, Región de Arica y Parinacota. Advanced Mining Technology Center. Universidad de Chile. Gaytán, D., 2012. Concesiones Volcán Tacora y Licancura III. Informe y presentación a Ministerio de Energía, 08 de Agosto de 2012. Gaytán, D., 2013. Estado actual de las actividades de exploración comprometidas In Informe de avance de actividades Concesión de Exploración “Volcán Tacora”, 29 de Marzo de 2013. Giggenbach, W.F., 1992. Relative importance of thermodynamic and kinetic processes in governing the chemical and isotopic composition of carbon gases in highheatflow sedimentary basins. Geochem. Cosmochim. Acta 61, 3763–3785. Gritto, R., Jarpe, S.P., 2014. Temporal variations of Vp/Vs ratio at the Geysers geothermal field, USA. Geothermics 52, 112–199. https://doi.org/10.1016/j.geothermics.2014. 01.012. Hantke, G., 1939. Übersicht über die Vulkanische Tätigkeit vom April bis Dezember 1938. Z. Dtsch. Geol. Ges. 91, 757–765. Ito, H., DeVilbiss, J., Nur, A., 1979. Compressional and shear waves in saturated rock during water-steam transition. J. Geophys. Res. 84, 4731–4735. Jain, A.K., 2010. Data clustering: 50 years beyond K-means. Pattern Recogn. Lett. 31 (I.8), 651–666. Jousset, P., Haberland, C., Bauer, K., Arnason, K., 2011. Hengill geothermal volcanic complex (Iceland) characterized by integrated geophysical observations. Geothermics 40 (1), 1–24. http://doi.org/10.1016/j.geothermics.2010.12.008. Klein, F.W., 1978. Hypocenter Location Program HYPOINVERSE. U.S. Geol. Surv. OpenFile Repp. 78–694 598. Koulakov, I., 2013. Studying deep sources of volcanism using multiscale seismic tomography. J. Volcanol. Geotherm. Res. 257, 205–226. http://doi.org/10.1016/j. jvolgeores.2013.03.012. Kuhn, M., 2004. Reactive Flow Modelling of Hydrothermal Systems. Springer Verlag. Lahsen, S., Muñoz, N., Parada, M.A., 2010. Geothermal development in Chile. In: Proceedings World Geothermal Congress. Lees, J.M., 2007. Seismic tomography of magmatic systems. J. Volcanol. Geotherm. Res. 167, 37–56. http://doi.org/10.1016/j.jvolgeores.2007.06.008. Lees, J.M., Wu, H., 2000. Poisson's ratio and porosity at coso geothermal area, California. J. Volcanol. Geotherm. Res. 95, 157–173. Londoño, J.M., Sudo, Y., 2003. Velocity structure and a seismic model for Nevado del Ruiz Volcano (Colombia). J. Volcanol. Geotherm. Res. 119 (1–4), 61–87.

6. Conclusions - The main fault inferred in the study region (Challaviento fault) corresponds to the primary influence on the system's fluid mobility. This can be corroborated analyzing its spatial relation with the fumarolic fields located on the NW flank of the volcanic edifice (Figs. 6 and 8). The horizontal projection (Fig. 11) indicates that this fault is also related to the hydrothermal alteration observed in the area. In the same way, this fault could also control the emplacement and location of the magmatic-hydrothermal reservoir. - We determined intermediate ranges of Vp/Vs ratios and percentage variations of Vp in relation to the existence and location of clay levels. According to the observed clay level in the region, Vp/Vs values are in the range of 1.70–1.74 exclusively for areas located over the magmatic-hydrothermal reservoir and with an average of m = 1.712. The percentage variations of Vp indicate that those values are in the range of −6.19 – 0.44%. - The groups observed in the cluster analysis (Fig. 9) allow to determine, using the Vp percentage variations, features that are not reached analyzing Vp/Vs only (Figs. 6 and 8). To establish lithological differences, it is imperative to analyze the elastic parameters independently, thus complementing the delimitation of the domains observed in the seismic tomography. Our analysis allows to define clay levels in volcanic systems, even in absence of surface manifestations. A similar study regarding clay levels has been documented in Tinguiririca Volcano (Pavez et al., 2016). In that case, the clay level was located directly below Los Humos fumarolic field and above of a lobe shape magmatic reservoir. This is in accordance with the results in the present research. - All the previous evidence suggests the existence of a degassing magmatic-hydrothermal reservoir located between sea level and 2 km depth. This reservoir is consistent with the magmatic-hydrothermal signature proposed by Capaccioni et al. (2011). - Travel time seismic tomography, which consider the local seismic activity, allows to constrain the different geometric domains of the conceptual model of the magmatic-hydrothermal system hosted in the Tacora volcano and the surrounding area (Capaccioni et al., 2011). The information derived from the seismicity distribution (i.e. ΔVp/Vp (%) and Vp/Vs structure) was used to define the spatial location of the recharge areas, magmatic reservoirs, degasification zones and clay levels in a 3D model of Tacora magmatic-hydrothermal system. Finally, the obtained model is coherent with previous research developed by Gaytán (2012, 2013) and with the INFINERGEO SpA's magnetotelluric study which can be found in Soyer and Hallinan (2012).

Acknowledgments The development of this article was funded by the doctoral fellowship Mecesup, project UCH0708 and supported by Advanced Mining Technology Center (AMTC). We appreciate the collaboration and support of INFINERGEO SpA. We thank to the free version of Leapfrog®: 3D & 2D images were generated by using the software Leapfrog Viewer®.

Appendix A. Supplementary data Supplementary data to this article can be found online at https:// doi.org/10.1016/j.jsames.2019.102247. 11

Journal of South American Earth Sciences 94 (2019) 102247

C. Pavez, et al.

Workshop on Geothermal Reservoir Engineering. Standford University, Stanford, California, January, pp. 25–27. Soyer, W., Hallinan, S., 2012. Magnetotelluric (MT) Survey, 3D Modelling and Gravity Study Tacora, Chile for Infinergeo SpA by Fugro Electro Magnetics Italy. Milan with Fugro Interra S.A., Santiago. Tanaka, S., Hamaguchi, H., Nishimura, T., Yamawaki, T., Ueki, S., Nakamichi, H., Tsutsui, T., Miyamachi, H., Matsuwo, N., Oikawa, J., Ohminato, T., Miyaoka, K., Onizawa, S., Mori, T., Aizawa, K., 2002. Three‐dimensional P‐wave velocity structure of Iwate volcano, Japan from active seismic survey. Geophys. Res. Lett. 29 (10). https://doi. org/10.1029/2002GL014983. Tassi, F., Martínez, C., Vaselli, O., Capaccioni, B., Viramonte, J., 2005. Light hydrocarbons as redox and temperature indicators in the geothermal field of El Tatio (northern Chile). Appl. Geochem. 20 (11), 2049–2062. Taran, Y.A., Pokrovsky, B.G., Esikov, A.D., 1989. Deuterium and oxygen-18 in fumarolic steam and amphiboles from some Kamchatka volcanoes: “andesitic waters”. Dokl. Akad. Nauk SSSR 304, 440–443. Tassi, F., Aguilera, F., Darrah, T., Vaselli, O., Capaccioni, B., Poreda, R.J., Delgado Huertas, A., 2010. Fluid geochemistry of hydrothermal systems in the AricaParinacota, Tarapacá and Antofagasta regions (northern Chile). J. Volcanol. Geotherm. Res. 192 (1–2), 1–15. Thierer, P.O., Flueh, E.R., Kopp, H., Tilmann, F., Comte, D., Contreras, S., 2005. Local earthquake monitoring offshore Valparaiso Chile. N. Jb. Geol. Paläont. (Abh.) 236, 173–183. Toksoz, M.N., Cheng, C.H., Timur, A., 1976. Velocities of seismic waves in porous rocks. Geophysics 41, 621–645. Tryggvason, A., Rögnvaldsson, S., Flóvenz, O., 2002. Three-dimensional imaging of the P and S wave velocity structure and earthquake locations beneath Southwest Iceland. Geophys. J. Int. 151, 848–866. Vanorio, T., Virieux, J., Capuano, P., Russo, G., 2005. Three dimensional seismic tomography from P waves and S wave microearthquake travel times and rock physics characterization of the Campi Flegrei Caldera. J. Geophys. Res. 110 (B3, 14.). Walck, M.C., 1988. Three-dimensional variations in shear structure and Vp/Vs for the Coso region, California. J. Geophys. Res. 93, 2047–2052. Wörner, G., Hammerschmidt, K., Henjes-Kunst, F., Lezaun, J., Wilke, H., 2000. Geochronology (40Ar/39Ar, K-Ar and He-exposure ages) of Cenozoic magmatic rocks from Northen Chile (18°-22°S): implications for magmatism and tectonic evolution of the Central Andes. Rev. Geol. Chile 205–240. Yoshikawa, M., Sudo, Y., 2004. Three Dimensional Seismic Velocity Structure beneath the Otake-Hatchobaru Geotermal Area at Kuju Volcano in Central Kyushu. Annuals of Disas. Prev. Res. Inst., Kyoto Univ., Japan No 47 B.

Montecinos, F., 1970. Geología del yacimiento de azufre del volcán Tacora, Departamento de Arica, vol. 25. Instituto de Investigaciones Geológicas, Santiago, pp. 41. Muksin, U., Bauer, K., Haberland, C., 2013. Seismic Vp and Vp/Vs structure of the geothermal area around Tarutung (North Sumatra, Indonesia) derived from local earthquake tomography. J. Volcanol. Geotherm. Res. 260, 27–42. Nakajima, J., Matsuzawa, T., Hasegawa, A., Zhao, D., 2001. Three dimensional structure of Vp, Vs, and Vp/Vs beneath northeastern Japan: implications for arc magmatism and fluids. J. Geophys. Res. 106, 21843–21857 N B10. Nicholson, K., 1993. Geothermal Fluids: Chemistry and Exploration Techniques. Springer Verlag. O'Connell, R.J., Budiansky, B., 1974. Seismic velocities in dry and saturated cracked solids. J. Geophys. Res. 79, 5412–5426. Pavez, C., Tapia, F., Comte, D., Gutierrez, F., Lira, E., Charrier, R., Benavente, O., 2016. Structural characterization of the hydrothermal system of the Tinguiririca volcanic complex, Central Chile, using passive seismic tomography. J. Volcanol. Geotherm. Res. 310, 107–117. Pavez, C., 2016. Modelamiento conceptual y numérico termal de sistemas hidrotermales basados en tomografía sísmica en los volcanes Tacora y Tinguiririca, Chile. PhD Thesis. Universidad de Chile In Spanish. Poletto, F., Farina, B., Carcione, J.M., 2018. Sensivity of seismic properties to temperature variations in a geothermal reservoir. Geothermics 76, 149–163. Roecker, S.W., 1982. Velocity structure of the Pamir-Hindu Kush region: possible evidence of subducted crust. J. Geophys. Res. 87, 945–959. https://doi.org/10.1029/ JB087iB02p00945. 649. Roecker, S.W., Yen, Y.H., Tsai, Y., 1987. Three dimensional P and S wave velocity structure beneath Taiwan: deep structure beneath an arc- continent collision. J. Geophys. Res. 92, 547–570. https://doi.org/10.1029/JB092iB10p10547. Roecker, S.W., Sabitova, T.M., Vinnick, L.P., Burmakov, Y.A., Golvanov, M.I., Mamatkanova, R., Murinova, L., 1993. Three dimensional elastic wave velocity structure of the western and central Tien Shan. J. Geophys. Res. 98 (B9), 779–795. https://doi.org/10.1029/93JB01560. Rojas, A., 2008. Modelamiento para el análisis de la variación hidrogeológica espacial del acuífero La Yarada- Tacna. Master Tesis. Universidad Nacional de Ingeniería. Rousseeuw, P.J., 1987. Silouhettes: a graphical aid to the interpretation and validation of cluster analysis. J. Comput. Appl. Math. 20, 53–65. Saccorotti, G., Piccinini, D., Zupo, M., Mazzarini, F., Chiarabba, C., Agostinetti, N.P., Licciardi, A., Bagagli, M., 2014. The deep structure of the Larderello – travale geothermal field (Italy) from integrated, passive seismic investigations. Energy Procedia 59, 227–234. https://doi.org/10.1016/j.egypro.2014.10.371. Simiyu, S.M., 1999. Seismic Velocity Analysis in the Olkaria Geothermal Field. 24th

12