# Geofluid Systems of Koryaksky-Avachinsky Volcanoes (Kamchatka, Russia)

Authors: A. Kiryukhin; V. Lavrushin; P. Kiryukhin; P. Voronin
Journal: Geofluids
 (2017)

---

## Abstract

The Koryaksky-Avachinsky volcanogenic basin, which has an area of 2530 km2, is located 25 km from Petropavlovsk-Kamchatsky City and includes five Quaternary volcanoes (two of which, Avachinsky (2750 masl) and Koryaksky (3456 masl), are active), and is located within a depression that has formed atop Cretaceous basement rocks. Magma injection zones (dikes and chamber-like shapes) are defined by plane-oriented clusters of local earthquakes that occur during volcanic activity (mostly in 2008–2011) below Koryaksky and Avachinsky volcanoes at depths ranging from −4.0 to −2.0 km and +1.0 to +2.0 km, respectively. Water isotopic (δD, δ18O) data indicate that these volcanoes act as recharge areas for their adjacent thermal mineral springs (Koryaksky Narzans, Isotovsky, and Pinachevsky) and the wells of the Bystrinsky and Elizovo aquifers. Carbon δ13С data in СО2 from CO2 springs in the northern foothills of Koryaksky Volcano reflect the magmatic origin of CO2. Carbon δ13С data in methane CH4 reservoirs penetrated by wells in the Neogene-Quaternary layer around Koryaksky and Avachinsky volcanoes indicate the thermobiogenic origin of methane. Thermal-hydrodynamic TOUGH2 conceptual modeling is used to determine what types of hydrogeologic boundaries and heat and mass sources are required to create the temperature, pressure, phase, and CO2 distributions observed within the given geological conditions of the Koryaksky-Avachinsky volcanic geofluid system.


---

## Body

## 1. Introduction

Avachinsky and Koryaksky volcanoes are located 25–30 km from the city of Petropavlovsk-Kamchatsky, which has a population of approximately 250 thousand people, contains in its interior significant resources of underground heat and groundwater, and represents a potential danger (see AppendixA with Figures). The cone of Avachinsky Volcano was formed 3500 years ago and has produced 15 historical eruptions since 1737. Koryaksky Volcano experienced a significant increase in fumarolic activity in 2008-2009 [1]. The Koryaksky-Avachinsky Basin contains two fresh groundwater reservoirs (Elizovo and Bystrinsky), which are used as the water supply of Petropavlovsk-Kamchatsky and its neighboring cities. The Ketkino low-temperature geothermal field is used for balneological purposes; it also features projects of the EGS utilization of the magma chamber of Avachinsky Volcano. Nevertheless, the geofluid system of the Koryaksky-Avachinsky volcanoes is not fully understood in terms of its heat source distributions and the types of hydrologic boundaries on the surfaces of the volcanoes, which prevents more extensive use of its potential resources.Magma injections below active volcanoes have been a significant focus of recent studies [2–4], which revealed their hydrofracturing nature and the plane-oriented swarms of earthquakes with which they are associated. Thus, in this paper, we used plane-oriented earthquake clusters (retrieved from the Kamchatka Branch of the Geophysical Survey Russia Academy of Sciences (KB GS RAS) catalogs 2001–2016) to track dikes injected below the above-mentioned volcanoes [5, 6]. The magma injection zones (i.e., dikes and chamber-like shapes) defined using this method were defined as heat sources for their adjacent hydrogeological reservoirs and surface features (hot springs and fumaroles) in terms of thermal-hydrodynamic models.Thus, this paper focuses on the Koryaksky-Avachinsky volcanic cluster and aims to achieve the following objectives: (1) using seismic data to identify dike swarms and magma chambers (i.e., potential heat sources); (2) using deep wells, thermal features, and the geochemistry of cold springs (i.e., water and carbon isotopes and water and gas chemistry) and temperature data to estimate the pressure/temperature/phase parameters of reservoirs and their mass/heat recharge conditions; (3) using TOUGH2-modeling, based on the above-mentioned data, to verify and estimate the role of Koryaksky Volcano as an injector of magma and cold water into adjacent structures and the creation of hydrothermal circulation and geothermal and gas reservoirs. It is also noteworthy that the rate and extent of water injection into the plumbing systems of active volcanoes are crucial to prevent or trigger catastrophic eruptions (i.e., hundreds of km3 in size) due to heat overbalance between magma and water injection rates, as was shown earlier by Fournier and Pitt [8] using Yellowstone caldera as an example.


## 2. Brief Review of Studies of Magma Plumbing Systems of Volcanoes

Volcanoes and crustal magma chambers result from the ascent of magma from primary magma chambers at depths of 150–200 km [11]. Crustal magma chambers are fed by magma from primary magma reservoirs (zones of magma accumulation) and supply magma needed for the formation of inclined dikes, sills, and dikes, including magma-feeding dikes that produce the eruption of magma onto the ground surface for the generation of volcanic edifices themselves. The formation of sills and dikes during initial emplacement through the viscous shells of magma chambers expands the chambers. Many dikes are “frozen” in their host rocks and do not reach the ground, whereas others change their orientations during propagation and become sills. A review of paleovolcanological evidence and the existing thermohydrodynamic models of magma chambers shows the following characteristics (Gudmudsson, 2012): (1) many crustal magma chambers are produced from sills; (2) active crustal magma chambers inject magma into associated volcanoes, with most chambers existing in a partially molten (pore-like) state; (3) ellipsoids are thermally stable shapes of magma chambers; (4) any discussion of hydraulic fracturing and the injection of dikes should consider the excess pressure of magma (and the strength of the host rocks needed to resist hydrofracturing, which is 1–6 MPa); (5) the excess pressure during eruptions decreases in an exponential manner, until the dike loses its hydraulic connection to the chamber; (6) the pressure of magma in the chamber during a period of repose must equal the overall lithostatic pressure; (7) the conditions that lead to hydraulic fracture at the top of a magma chamber are reached in two ways, either by increasing the absolute pressure of the magma in the chamber due to injection from the magma reservoir or by decreasing the horizontal stress due to regional extension; and (8) the local stress field around the magma chamber depends on the shape and depth of the chamber, as well as the geomechanical properties and layering of the host rocks.Hydraulic fracturing is more difficult when the fluid is subject to “freezing” because in this case one needs greater discharge rates to maintain thermal balance above the solidus. For this reason, magma injections in the form of dikes should have larger apertures (up to 1–10 m), which can be deduced from both theoretical calculations and geological observations.The formation of new fissures by hydraulic fracturing and the renewal of activity on preexisting faults are accompanied by seismicity due to the generation of shear fissures in the zone immediately adjacent to the main aseismic zone of a hydraulic fracture fissure (i.e., a seismicity trigger). For this reason, it may be hypothesized that the planes that host microearthquake clusters have the same orientation as the hydraulic fracture fissures produced during magma injection (emplacement). The magnitudes of strike-slip earthquakes with slip amplitudes of 0.1 mm to 1 cm and fissure lengths of a few hundred meters are estimated to range from 1 to 2. This is consistent with the sensitivity of the local seismograph networks that are operated in the area of the Koryaksky-Avachinsky volcanic cluster.The results of hydrogeomechanical CFRAC modeling [12, 13] show that it is feasible that plane-oriented clusters of earthquakes beneath active volcanoes may indicate processes of magma fracking or dike formation [6, 14].The approach to this problem is also motivated by the observations reported by Sigmundsson et al. [3], who described the injection of magma from the magma chamber beneath the Bardarbunga central-type volcano, Iceland, which occurred in August 2014 and was accompanied by a dike that was propagated at a distance of 50 km. Bardarbunga volcano (which has a caldera that is 8 × 11 km in size and produced 23 eruptions during the last 1100 years) is located in the central part of the Icelandic rift zone (with an extension rate of 19 mm/year). The magma chamber of Bardarbunga volcano is assumed to be located at a depth of 10–15 km. The volume of lava erupted in 2014 is estimated to be 1.4 km3 (thus representing the largest lava eruption in Iceland since 1783). Seismic event data were used to identify the shape of the 50 km long segmented dike, which included 11 plane-oriented earthquake clusters (with the number of earthquakes, which had magnitudes of up to 5, ranging from 57 to 1181). The injection of the dike lasted for 22 days. It is noteworthy that the first dike segment was normal to the caldera boundary (along the axes of least stress), while the next dike that was generated coincided with the strike of the rift zone. The eruption started at the last dike segment and lasted for 6 months. GPS and InSAR data yielded an estimate of the volume of the dike system of 0.6 km3 and that of a dike wall opening of 320 cm; the decrease in the volume of the magma chamber was estimated to be 0.3 km3.Ground deformation that occurred during the 2012-2013 eruption of Tolbachik Volcano (Kamchatka) (9 months, 0.54 km3) was estimated based on multiple satellite-based radar observations (InSAR) [15]. Climatic conditions (snow cover) made comparisons of interferograms only feasible for the surveys conducted in August-September 2012 and 2013. The results of 3D geomechanical modeling indicated that the deformation observed during the time span of interest (1 year) was the result of the emplacement of a radial dike into Tolbachik Volcano (with a dip angle of 80° to the west-northwest). Later, a more detailed analysis of seismic data preceding the November 27, 2012, Tolbachik lava eruption [5] revealed that the injection of magma resulted in a series of dikes trending west-northwestward at absolute depths ranging from −4 to +3 km in a zone situated to the southeast of the Ploskii Tolbachik Volcano edifice. The dikes penetrated into a nearly horizontal permeable zone at an absolute depth of approximately zero, producing sills and emplacing a magma-conducting dike along the top of the zone of cinder cones (with a dip angle of 50° toward the azimuth at 300°), which was located 5.5 km from the epicenter of the initial magma injection.Another example was presented by Dumont et al. [4], who identified numerous dike injections (in a 60 km × 5 km area) in the Afar rift system produced during 2005–2010 using InSAR and seismic data. A mid-segment magma chamber at a depth of ~4 km and a crustal chamber at a depth of more than 15 km are considered to be sources of the dike injections, which occurred in different directions above the source magma chambers along rift-zone axes.The above discussion shows that the injection of magma beneath volcanoes during the generation of dikes and sills due to magma-driven fracturing is largely analogous to the injection of fluids into wells, which is associated with the subsequent hydraulic fracturing and fissure generation in their host formations. The data from observations of fluids injected into wells in oil, gas, and geothermal fields have been widely used to evaluate the geometry of fracture systems and the state of stress in the reservoirs using seismic data [16]. This approach may also be applied in order to understand magmatic fracturing occurring beneath active volcanoes.Finally, after magmatic activity ends or shifts to another volcanic funnel, water comes there. An example of the high-rate cold springs (0.7–6.8 m3/s) located in the eastern foothills of the Cascades (Oregon) was given by James et al. [17]. Their water isotopic (δD, δ18O) data indicated that the summits of the Cascade volcanoes serve as a recharge area for these springs, which are located 30–60 km apart. Thus, the magmatic plumbing systems of dormant volcanoes switch to inject cold water into adjacent hydrological reservoirs.


## 3. Geological Setting

### 3.1. Hydrogeological Stratification

The Koryaksky-Avachinsky volcanogenic basin, which has an area of 2530 km2, contains five Quaternary volcanoes (two of which, Avachinsky (2750 masl) and Koryaksky (3456 masl), are active) and subbasins of volcanogenic and sedimentary Neogene-Quaternary deposits that are up to 1.4 km thick (Figure 1). The basin is located in a depression that has formed atop Cretaceous basement rocks and is generally characterized by a low-temperature gradient of 24°C/km.Figure 1
Geological map of the Koryaksky-Avachinsky volcanogenic basin. (1) The summits of the Avachinsky, Koryaksky, Kozelsky, Arik, and Aag volcanoes; (2) Avachinsky, Koryaksky, Kozelsky volcanoes, and their eruptive products; (3) Pinachevsky extrusions Q2-3; (4) thermal features (for details, see Table 1): FA, fumaroles on Avacha Volcano; FK, fumaroles on Koryaksky Volcano; K1, K2, K3, K7, and K8, thermal mineral springs of Koryaksky Narzan; IS, Izotovsky; VD, Vodopadny; CH, Chistinsky; Va, Vakinsky; (5) deep hydrogeological wells; (6) KB GS RAS seismograph stations; (7) dikes traced at −3000 masl below Koryaksky Volcano and 1500 masl below Avachinsky Volcano; (8) glaciers. Notes: (1) the isolines show the topographic surface, and the ticks along the axes represent intervals of 5 km. (2) For the extended map, including positions of all deep wells, see Figure 1 in Kiryukhin et al. [7].Table 1
Hot and cold springs of the Koryaksky-Avachinsky volcanogenic basin.
Name
##
E
N
Z masl
Cold Spring C5
CS
158.8070
53.3243
910
Cold spring C7
C7
158.7975
53.3434
899
Chistinsky Narzan
CH
158.7287
53.4367
602
Avacha Fumarole Left
FA
158.8309
53.2559
2720
Avacha FumaroleNZ frac
FA
158.8301
53.2563
2672
Avacha Fumarole Central
FA2
158.8339
53.2546
2720
Koryaksky Fumarole
FK
158.6960
53.3300
2900
Isotovsky S
IS
158.7322
53.3920
786
Koryaksky Narzan 1
KI
158.7429
53.3882
920
Koryaksky Narzan 2
K2
138.7552
33.3924
800
Koryaksky Narzan 3
K3
158.7439
53.3852
922
Koryaksky Narzan 7
K7
158.7535
53.3812
899
Koryaksky Narzan 8
K8
158.7557
53.3733
967
Vodopadny
VD
158.7400
53.4215
569The basin basement comprises Upper CretaceousК2, which is represented by metamorphic rocks of metasandstone, metasiltstone, and phyllites with interbedded shales and microquartzites. The porosity values of these rocks are a few tenths of a percent. The matrix of these rocks has a low permeability of 0.001–0.01 mD. However, there are fracture zones that can be used to test the flows of reservoir water, which have flow rates ranging from 0.05 to 6.1 kg/s at a depth interval of 1438–1490 m (well E1, Figure 1). The average thermal conductivity of the Cretaceous deposits is 2.8 W/m С°.The Miocene (N1)-Quaternary aquifer complex is composed of pyroclastic and volcanogenic-sedimentary formations. Their porosity values are quite high and range from 0.36 to 0.48; their well productivity indexes range from 0.004 kg/s/bar (well P2) to 0.1 kg/s/bar (well GK1 Pinachevskaya). The average thermal conductivity of the Paleogene-Quaternary deposits is 1.5 W/m С°. The aquifer system of the Pinachevsky Q1-2 extrusions is composed of andesite and rhyolite extrusions and includes vent andesite, dacite, and rhyolite formations (with thicknesses of more than 200–500 m). According to laboratory studies, their porosity is 0.12 and their permeability is 24 mD.The artesian volcanogenic basin (AVB) includes the aquifers of water-glacial formations, which comprise Holocene alluvial deposits, Upper Pleistocene-Holocene marine and alluvial-marine horizons, Upper Pleistocene glacial and fluvioglacial complexes, Holocene aquifers, and a proluvial and diluvial-proluvial complex. Its permeability ranges from 10 to 3000 mD, according to the results of the well-testing of the Bystrinsky freshwater reservoir.


### 3.2. Magmatic Activity beneath Koryaksky-Avachinsky Volcanoes from 2000–2016

#### 3.2.1. Seismic Data and Method of Plane-Oriented Clusters Identification

In this study, we assume that the emplacement of magma in a fractured medium beneath active volcanoes is analogous to the injection of fluids into wells with subsequent hydraulic fracturing occurring in their host formations. Six stations can record seismicity in the Koryaksky-Avachinsky volcanic cluster (Figure1). The absolute uncertainty in the locations of the hypocenters and epicenters of microearthquakes in this area is estimated to be 1 km [18, 19]. A total of 5160 earthquakes have been recorded by the Kamchatka Branch Geophysical Survey Russia Academy of Sciences (KB GS RAS) as having occurred in the edifices and basement of the Koryaksky-Avachinsky volcanoes during the period between January 2000 and July 2016. Our treatment of this dataset using the method outlined below yielded 1540 earthquakes comprising 204 plane-oriented clusters for the period between January 2008 and February 2016.Cluster identification was carried out using our Frac-Digger program (RU #Reg. 2016616880). The following is a brief explanation of the algorithm used in this program. The first element of the cluster is removed from the initial list during each iteration. The following criteria are used to include a new event in a cluster: (1) a time difference (δt = 1 day); (2) a distance difference in the horizontal plane (δR < 6 km); and (3) a requirement of a nearly planar orientation (i.e., a distance from the event to the plane (δZ) that is less than 200 m). When the resulting cluster contains more than five elements (N>5), that cluster is treated as completed and is added to the list of plane-oriented clusters. All elements of a resulting cluster are removed from the initial list of elements (in cases when the cluster size > 5). This procedure is then reiterated until the initial list of elements is exhausted.The calculation of the parameters of a plane-oriented cluster is based on a list of cluster elements. Each elementi contains the coordinates (xi, yi, zi). For N (the number of elements in the cluster) points with coordinates (xi, yi, zi), we can find the equation of the fitting plane z=ax+by+c using the least-squares method. The solution thus reduces to solving a set of linear equations as follows:(1)∑xi2∑xiyi∑xi∑xiyi∑yi2∑yi∑xi∑yinabc=∑xizi∑yizi∑zi.These equations are then solved using Cramer’s rule. In this way, we obtain the coefficientsa, b, c for the equation of a plane, which is defined as z=ax+by+c. The next step is to find the unit vector that is normal to the fitting plane n=(a/Δ,b/Δ,-1/Δ), where Δ is the determinant of the equations that result in the following geological parameters: dip angle β=arccos⁡(1/Δ)∗180/π and the azimuth of dip α=arctan⁡(a/b)∗180/π. The analysis of the sensitivity of this algorithm, where plane-oriented clusters are selected according to the criteria of temporal and spatial proximity, indicates that the above criteria produce selection results that are both physically and geologically reasonable (Appendix B). The results of the hydrogeomechanical CFRAC modeling [12, 13] performed indicate that it is feasible that plane-oriented clusters of earthquakes beneath active volcanoes may indicate magma fracking or dike formation processes [6, 14].


#### 3.2.2. Distributions of Dikes, Sills, and Magma Chambers

The analysis of the seismic activity at Koryaksky Volcano (Figure1) revealed the following geomechanical features (Figures 1–4): (1) the 2008–2009 summit steam-gas eruption of Koryaksky Volcano was accompanied by 153 plane-oriented earthquake clusters that are interpreted here to represent zones where dikes and sills were emplaced during magma injection; (2) the precursory period of this eruption began with magma filling the crustal chamber (the top is at an absolute depth of −3 km and the chamber is 2.5 km wide) near the southwestern base of Koryaksky Volcano (from July 2008 to January 2009); (3) magma was injected into a nearly meridional zone (7.5 by 2.5 km, with absolute depths ranging from −2 to −5 km) in the northern sector of Koryaksky Volcano simultaneously during the most intense period of the summit steam-gas eruption (from February 2009 to March 2010); (4) following the saturation of the plumbing system beneath Koryaksky Volcano, magma was injected into the cone of Avachinsky Volcano (2010–2016, at absolute depths ranging from 1 to 2 km).Figure 2
Geological cross-sectional AB (see Figure1). (1) Basement with metamorphic Upper Cretaceous K2 rocks in its upper part; (2) Miocene (N1)-Quaternary aquifer complex; (3) Koryaksky Volcano (or/and Avachinsky Volcano in Figures 3 and 4); (4) traces of dikes and sills, defined from plane-oriented MEQ clusters; (5) assumed main feeding magma channels of Koryaksky Volcano (or/and Avachinsky Volcano in Figures 3 and 4); (6) earthquake hypocenters, in which their sizes are proportional to their magnitudes (KB GS RAS data 2001–2016). The ticks along the axes represent intervals of 1 km, and axes are labeled in meters.Figure 3
Geological cross-sectional CD (see Figure1). Legend is the same as in Figure 2, and axes are labeled in meters.Figure 4
Geological cross-sectional EF (see Figure1). Legend is the same as in Figure 2, and axes are labeled in meters.


## 3.1. Hydrogeological Stratification

The Koryaksky-Avachinsky volcanogenic basin, which has an area of 2530 km2, contains five Quaternary volcanoes (two of which, Avachinsky (2750 masl) and Koryaksky (3456 masl), are active) and subbasins of volcanogenic and sedimentary Neogene-Quaternary deposits that are up to 1.4 km thick (Figure 1). The basin is located in a depression that has formed atop Cretaceous basement rocks and is generally characterized by a low-temperature gradient of 24°C/km.Figure 1
Geological map of the Koryaksky-Avachinsky volcanogenic basin. (1) The summits of the Avachinsky, Koryaksky, Kozelsky, Arik, and Aag volcanoes; (2) Avachinsky, Koryaksky, Kozelsky volcanoes, and their eruptive products; (3) Pinachevsky extrusions Q2-3; (4) thermal features (for details, see Table 1): FA, fumaroles on Avacha Volcano; FK, fumaroles on Koryaksky Volcano; K1, K2, K3, K7, and K8, thermal mineral springs of Koryaksky Narzan; IS, Izotovsky; VD, Vodopadny; CH, Chistinsky; Va, Vakinsky; (5) deep hydrogeological wells; (6) KB GS RAS seismograph stations; (7) dikes traced at −3000 masl below Koryaksky Volcano and 1500 masl below Avachinsky Volcano; (8) glaciers. Notes: (1) the isolines show the topographic surface, and the ticks along the axes represent intervals of 5 km. (2) For the extended map, including positions of all deep wells, see Figure 1 in Kiryukhin et al. [7].Table 1
Hot and cold springs of the Koryaksky-Avachinsky volcanogenic basin.
Name
##
E
N
Z masl
Cold Spring C5
CS
158.8070
53.3243
910
Cold spring C7
C7
158.7975
53.3434
899
Chistinsky Narzan
CH
158.7287
53.4367
602
Avacha Fumarole Left
FA
158.8309
53.2559
2720
Avacha FumaroleNZ frac
FA
158.8301
53.2563
2672
Avacha Fumarole Central
FA2
158.8339
53.2546
2720
Koryaksky Fumarole
FK
158.6960
53.3300
2900
Isotovsky S
IS
158.7322
53.3920
786
Koryaksky Narzan 1
KI
158.7429
53.3882
920
Koryaksky Narzan 2
K2
138.7552
33.3924
800
Koryaksky Narzan 3
K3
158.7439
53.3852
922
Koryaksky Narzan 7
K7
158.7535
53.3812
899
Koryaksky Narzan 8
K8
158.7557
53.3733
967
Vodopadny
VD
158.7400
53.4215
569The basin basement comprises Upper CretaceousК2, which is represented by metamorphic rocks of metasandstone, metasiltstone, and phyllites with interbedded shales and microquartzites. The porosity values of these rocks are a few tenths of a percent. The matrix of these rocks has a low permeability of 0.001–0.01 mD. However, there are fracture zones that can be used to test the flows of reservoir water, which have flow rates ranging from 0.05 to 6.1 kg/s at a depth interval of 1438–1490 m (well E1, Figure 1). The average thermal conductivity of the Cretaceous deposits is 2.8 W/m С°.The Miocene (N1)-Quaternary aquifer complex is composed of pyroclastic and volcanogenic-sedimentary formations. Their porosity values are quite high and range from 0.36 to 0.48; their well productivity indexes range from 0.004 kg/s/bar (well P2) to 0.1 kg/s/bar (well GK1 Pinachevskaya). The average thermal conductivity of the Paleogene-Quaternary deposits is 1.5 W/m С°. The aquifer system of the Pinachevsky Q1-2 extrusions is composed of andesite and rhyolite extrusions and includes vent andesite, dacite, and rhyolite formations (with thicknesses of more than 200–500 m). According to laboratory studies, their porosity is 0.12 and their permeability is 24 mD.The artesian volcanogenic basin (AVB) includes the aquifers of water-glacial formations, which comprise Holocene alluvial deposits, Upper Pleistocene-Holocene marine and alluvial-marine horizons, Upper Pleistocene glacial and fluvioglacial complexes, Holocene aquifers, and a proluvial and diluvial-proluvial complex. Its permeability ranges from 10 to 3000 mD, according to the results of the well-testing of the Bystrinsky freshwater reservoir.


## 3.2. Magmatic Activity beneath Koryaksky-Avachinsky Volcanoes from 2000–2016

### 3.2.1. Seismic Data and Method of Plane-Oriented Clusters Identification

In this study, we assume that the emplacement of magma in a fractured medium beneath active volcanoes is analogous to the injection of fluids into wells with subsequent hydraulic fracturing occurring in their host formations. Six stations can record seismicity in the Koryaksky-Avachinsky volcanic cluster (Figure1). The absolute uncertainty in the locations of the hypocenters and epicenters of microearthquakes in this area is estimated to be 1 km [18, 19]. A total of 5160 earthquakes have been recorded by the Kamchatka Branch Geophysical Survey Russia Academy of Sciences (KB GS RAS) as having occurred in the edifices and basement of the Koryaksky-Avachinsky volcanoes during the period between January 2000 and July 2016. Our treatment of this dataset using the method outlined below yielded 1540 earthquakes comprising 204 plane-oriented clusters for the period between January 2008 and February 2016.Cluster identification was carried out using our Frac-Digger program (RU #Reg. 2016616880). The following is a brief explanation of the algorithm used in this program. The first element of the cluster is removed from the initial list during each iteration. The following criteria are used to include a new event in a cluster: (1) a time difference (δt = 1 day); (2) a distance difference in the horizontal plane (δR < 6 km); and (3) a requirement of a nearly planar orientation (i.e., a distance from the event to the plane (δZ) that is less than 200 m). When the resulting cluster contains more than five elements (N>5), that cluster is treated as completed and is added to the list of plane-oriented clusters. All elements of a resulting cluster are removed from the initial list of elements (in cases when the cluster size > 5). This procedure is then reiterated until the initial list of elements is exhausted.The calculation of the parameters of a plane-oriented cluster is based on a list of cluster elements. Each elementi contains the coordinates (xi, yi, zi). For N (the number of elements in the cluster) points with coordinates (xi, yi, zi), we can find the equation of the fitting plane z=ax+by+c using the least-squares method. The solution thus reduces to solving a set of linear equations as follows:(1)∑xi2∑xiyi∑xi∑xiyi∑yi2∑yi∑xi∑yinabc=∑xizi∑yizi∑zi.These equations are then solved using Cramer’s rule. In this way, we obtain the coefficientsa, b, c for the equation of a plane, which is defined as z=ax+by+c. The next step is to find the unit vector that is normal to the fitting plane n=(a/Δ,b/Δ,-1/Δ), where Δ is the determinant of the equations that result in the following geological parameters: dip angle β=arccos⁡(1/Δ)∗180/π and the azimuth of dip α=arctan⁡(a/b)∗180/π. The analysis of the sensitivity of this algorithm, where plane-oriented clusters are selected according to the criteria of temporal and spatial proximity, indicates that the above criteria produce selection results that are both physically and geologically reasonable (Appendix B). The results of the hydrogeomechanical CFRAC modeling [12, 13] performed indicate that it is feasible that plane-oriented clusters of earthquakes beneath active volcanoes may indicate magma fracking or dike formation processes [6, 14].


### 3.2.2. Distributions of Dikes, Sills, and Magma Chambers

The analysis of the seismic activity at Koryaksky Volcano (Figure1) revealed the following geomechanical features (Figures 1–4): (1) the 2008–2009 summit steam-gas eruption of Koryaksky Volcano was accompanied by 153 plane-oriented earthquake clusters that are interpreted here to represent zones where dikes and sills were emplaced during magma injection; (2) the precursory period of this eruption began with magma filling the crustal chamber (the top is at an absolute depth of −3 km and the chamber is 2.5 km wide) near the southwestern base of Koryaksky Volcano (from July 2008 to January 2009); (3) magma was injected into a nearly meridional zone (7.5 by 2.5 km, with absolute depths ranging from −2 to −5 km) in the northern sector of Koryaksky Volcano simultaneously during the most intense period of the summit steam-gas eruption (from February 2009 to March 2010); (4) following the saturation of the plumbing system beneath Koryaksky Volcano, magma was injected into the cone of Avachinsky Volcano (2010–2016, at absolute depths ranging from 1 to 2 km).Figure 2
Geological cross-sectional AB (see Figure1). (1) Basement with metamorphic Upper Cretaceous K2 rocks in its upper part; (2) Miocene (N1)-Quaternary aquifer complex; (3) Koryaksky Volcano (or/and Avachinsky Volcano in Figures 3 and 4); (4) traces of dikes and sills, defined from plane-oriented MEQ clusters; (5) assumed main feeding magma channels of Koryaksky Volcano (or/and Avachinsky Volcano in Figures 3 and 4); (6) earthquake hypocenters, in which their sizes are proportional to their magnitudes (KB GS RAS data 2001–2016). The ticks along the axes represent intervals of 1 km, and axes are labeled in meters.Figure 3
Geological cross-sectional CD (see Figure1). Legend is the same as in Figure 2, and axes are labeled in meters.Figure 4
Geological cross-sectional EF (see Figure1). Legend is the same as in Figure 2, and axes are labeled in meters.


## 3.2.1. Seismic Data and Method of Plane-Oriented Clusters Identification

In this study, we assume that the emplacement of magma in a fractured medium beneath active volcanoes is analogous to the injection of fluids into wells with subsequent hydraulic fracturing occurring in their host formations. Six stations can record seismicity in the Koryaksky-Avachinsky volcanic cluster (Figure1). The absolute uncertainty in the locations of the hypocenters and epicenters of microearthquakes in this area is estimated to be 1 km [18, 19]. A total of 5160 earthquakes have been recorded by the Kamchatka Branch Geophysical Survey Russia Academy of Sciences (KB GS RAS) as having occurred in the edifices and basement of the Koryaksky-Avachinsky volcanoes during the period between January 2000 and July 2016. Our treatment of this dataset using the method outlined below yielded 1540 earthquakes comprising 204 plane-oriented clusters for the period between January 2008 and February 2016.Cluster identification was carried out using our Frac-Digger program (RU #Reg. 2016616880). The following is a brief explanation of the algorithm used in this program. The first element of the cluster is removed from the initial list during each iteration. The following criteria are used to include a new event in a cluster: (1) a time difference (δt = 1 day); (2) a distance difference in the horizontal plane (δR < 6 km); and (3) a requirement of a nearly planar orientation (i.e., a distance from the event to the plane (δZ) that is less than 200 m). When the resulting cluster contains more than five elements (N>5), that cluster is treated as completed and is added to the list of plane-oriented clusters. All elements of a resulting cluster are removed from the initial list of elements (in cases when the cluster size > 5). This procedure is then reiterated until the initial list of elements is exhausted.The calculation of the parameters of a plane-oriented cluster is based on a list of cluster elements. Each elementi contains the coordinates (xi, yi, zi). For N (the number of elements in the cluster) points with coordinates (xi, yi, zi), we can find the equation of the fitting plane z=ax+by+c using the least-squares method. The solution thus reduces to solving a set of linear equations as follows:(1)∑xi2∑xiyi∑xi∑xiyi∑yi2∑yi∑xi∑yinabc=∑xizi∑yizi∑zi.These equations are then solved using Cramer’s rule. In this way, we obtain the coefficientsa, b, c for the equation of a plane, which is defined as z=ax+by+c. The next step is to find the unit vector that is normal to the fitting plane n=(a/Δ,b/Δ,-1/Δ), where Δ is the determinant of the equations that result in the following geological parameters: dip angle β=arccos⁡(1/Δ)∗180/π and the azimuth of dip α=arctan⁡(a/b)∗180/π. The analysis of the sensitivity of this algorithm, where plane-oriented clusters are selected according to the criteria of temporal and spatial proximity, indicates that the above criteria produce selection results that are both physically and geologically reasonable (Appendix B). The results of the hydrogeomechanical CFRAC modeling [12, 13] performed indicate that it is feasible that plane-oriented clusters of earthquakes beneath active volcanoes may indicate magma fracking or dike formation processes [6, 14].


## 3.2.2. Distributions of Dikes, Sills, and Magma Chambers

The analysis of the seismic activity at Koryaksky Volcano (Figure1) revealed the following geomechanical features (Figures 1–4): (1) the 2008–2009 summit steam-gas eruption of Koryaksky Volcano was accompanied by 153 plane-oriented earthquake clusters that are interpreted here to represent zones where dikes and sills were emplaced during magma injection; (2) the precursory period of this eruption began with magma filling the crustal chamber (the top is at an absolute depth of −3 km and the chamber is 2.5 km wide) near the southwestern base of Koryaksky Volcano (from July 2008 to January 2009); (3) magma was injected into a nearly meridional zone (7.5 by 2.5 km, with absolute depths ranging from −2 to −5 km) in the northern sector of Koryaksky Volcano simultaneously during the most intense period of the summit steam-gas eruption (from February 2009 to March 2010); (4) following the saturation of the plumbing system beneath Koryaksky Volcano, magma was injected into the cone of Avachinsky Volcano (2010–2016, at absolute depths ranging from 1 to 2 km).Figure 2
Geological cross-sectional AB (see Figure1). (1) Basement with metamorphic Upper Cretaceous K2 rocks in its upper part; (2) Miocene (N1)-Quaternary aquifer complex; (3) Koryaksky Volcano (or/and Avachinsky Volcano in Figures 3 and 4); (4) traces of dikes and sills, defined from plane-oriented MEQ clusters; (5) assumed main feeding magma channels of Koryaksky Volcano (or/and Avachinsky Volcano in Figures 3 and 4); (6) earthquake hypocenters, in which their sizes are proportional to their magnitudes (KB GS RAS data 2001–2016). The ticks along the axes represent intervals of 1 km, and axes are labeled in meters.Figure 3
Geological cross-sectional CD (see Figure1). Legend is the same as in Figure 2, and axes are labeled in meters.Figure 4
Geological cross-sectional EF (see Figure1). Legend is the same as in Figure 2, and axes are labeled in meters.


## 4. Koryaksky-Avachinsky Basin Discharge/Recharge Conditions Based on Gas, Chemical, and Isotopic Data

Intervals of deep sampling wells (Table2(a)) that encountered the basement of the Koryaksky-Avachinsky Basin are characterized by their chemical composition, which corresponds to the dilution of oceanic sediment water due to the infiltration of meteoric waters. Therefore, these mostly represent dominantly Cl-Na waters with minor concentrations of НСО3, Са, Mg, and SO4. There is no indication of significant water metamorphism, and their Cl/Na ratios are close to seawater values. Increases in the concentrations of Ca and Mg may be related to CO2-leaching processes from the host rocks.Table 2
(a)
Water chemical compositions of deep wells of the Koryaksky-Avachinsky volcanogenic basin, according to sampling data collected from 1970–2014. Wells: GK1, G1, K6, R3, R2, 53, E1, KO1, K5, and K6 (well positions shown in Figure1). GK1, G1: data from G. Ryabinin (pers. com., 2011); R2, R3 (A. Evtukhov et al., 1995), and E1: Kopylova G. (pers. com., 2011). The remaining samples were collected by A. Kiryukhin and T. Rychkova, and analyses were performed in the Central Chemistry Lab of the Institute of Volcanology FEB RAS. TSiO2 and TNa-K: hydrogeochemical geothermometers (TSiO2=1309/(5.19-log⁡10(I3))-273, Fournier, 1977; and TNa-K=855.6/(log⁡10(Na/K)+0.8573)-273, White, 1970).
##
Data
рН
Н
С
О
3
-
C
l
-
S
O
4
-
Na+
K+
Ca2+
Mg2+
SiO2
Na/K
T
N
a
-
K
T
S
i
O
2
GK1
04.10.1970
7.5
207
2199
1050
70
241
15
148
G1
14.04.1988
8.8
502
4929
100
3450
48
24
24
72
42
K6
17.11.1990
8.5
2124
1175
25
119
24
46
47
65
98
K8
17.11.1990
8.3
14
6358
13
3125
91
598
234
44
34
85
96
R3
05.06.1992
7.6
144
12869
16
6667
167
679
271
61
40
75
111
R2
15.03.1993
8.2
518
6901
4308
0
110
122
11
43
E1
30.08.2011
12.1
68
344
10
285
47
36
2
2
6
53
05.07.2011
8.2
51
24
154
28
3
28
22
35
9
195
86
GK1
28.07.2014
6.9
60
6594
173
3387
132
778
0
59
26
105
110
K01
09.06.2014
8.2
26
5460
53
2759
75
689
0
26
37
80
74
K5
30.06.2014
8.9
16
2411
5
1336
37
180
0
39
36
81
91
K6
09.06.2014
7.2
139
2375
2
1337
53
208
0
52
25
105
104
(b)
Water chemical compositions (ppm) of hot and cold springs of the Koryaksky-Avachinsky volcanogenic basin, according to sampling data collected from 2011 to 2016, ##: denoted names, see Table1 for explanation; gas, l/s; water, kg/s. The samples were collected by A. V. Kiryukhin, and analyses were performed in the Central Chemistry Lab of the Institute of Volcanology FEB RAS. TSiO2 and TNa-K, hydrogeochemical geothermometers (TSiO2 = 1309/(5.19 − log⁡10(I3)) − 273, Fournier, 1977; and TNa-K = 855.6/(log⁡10(Na/K) + 0.8573) − 273, White, 1970)
##
Data
T°C
gas
w
рН
Н
С
О
3
-
C
l
-
S
O
4
-
Na+
K+
Ca2+
Mg2+
SiO2
Na/K
T
N
a
-
K
T
S
i
O
2
K8
09.07.2010
12
2
6.3
1072
21
40
56
11
86
137
104
4.9
281
139
K1
09.07.2010
10
10
6.4
687
40
78
50
9
60
101
84
5.7
257
128
K2
09.07.2010
11
15
6.2
877
45
91
62
12
67
131
94
5.0
278
134
K3
09.07.2010
11
6.5
667
36
77
47
9
62
97
82
5.4
265
126
IS
25.09.2010
48
0.10
7.4
787
33
1167
150
29
415
136
99
5.1
273
137
K7
25.09.2010
9
7
788
31
128
55
13
70
117
90
4.4
296
132
K2
25.09.2010
7.4
872
43
102
61
13
51
142
81
4.9
281
126
IS
28.07.2011
39
0.07
6.9
800
28
1114
323
103
401
2
66
3.1
359
115
K1
29.07.2011
6.6
678
32
115
64
14
76
92
40
4.6
291
92
K2
29.07.2011
14
6.3
866
41
134
80
21
84
120
51
3.8
322
103
VD
29.07.2011
42
6.5
429
<0.3
67
47
22
38
47
111
2.1
448
143
C5
30.07.2011
~1000
6.5
9
1
1
1
1
1
1
3
1.0
K1
14.08.2012
11
6.4
661
28
81
55
8
68
83
169
6.5
240
169
K8
14.08.2012
12
6.4
774
25
202
85
16
98
103
80
5.4
266
125
K2
14.08.2012
13
6.3
960
21
9
74
9
212
36
98
8.0
213
136
IS
14.08.2012
46
6.6
813
50
1536
162
91
501
137
54
1.8
499
106
K1
06.08.2013
12
15–20
6.5
696
33
98
67
10
63
90
94
6.9
231
134
K2
06.08.2013
13
10
6.2
857
43
53
73
13
69
109
103
5.5
263
139
IS
06.08.2013
46
0.20
6.5
852
35
1003
153
37
386
106
115
4.1
309
145
K3
07.08.2013
13
6.5
657
34
58
66
10
66
90
95
6.9
232
134
K1
08.08.2014
13
6.5
664
34
154
50
8
80
101
87
6.3
243
130
K2
08.08.2014
14
6.2
798
44
138
63
12
62
120
91
5.4
264
132
IS
08.08.2014
6.6
805
30
769
157
29
419
24
101
5.3
267
138
C7
09.08.2014
2
~1000
6.8
52
2
3
4
1
12
2
60
3.3
C5
09.08.2014
4
~1000
6.7
5
1
2
1
0
2
0
13
7.3
AV2166
10.08.2014
5.6
6
2
3
1
0
3
0
12
7.1
228
45
K2
05.08.2015
14
6.4
6.40
6.3
812
46
93
56
12
82
109
156
4.5
293
164
IS
05.08.2015
6.7
868
138
1287
129
39
481
173
156
3.3
350
164
K1
05.08.2015
6.2
562
37
106
58
9
74
81
155
6.5
238
163
K1
03.08.2016
6.2
555
36
86
50
10
70
79
72
5.3
269
120
K2
03.08.2016
5.9
748
46
99
61
14
88
105
77
4.4
295
123
IS
03.08.2016
6.3
658
28
1316
162
32
445
102
93
5.1
275
133
K7
04.08.2016
5.9
682
32
73
42
14
86
90
67
3.1
360
116
K8
04.08.2016
6
1053
14
48
60
12
104
119
86
5.0
277
129
VD
03.08.2016
6
236
1
15
23
9
43
1
125
2.7
393
150
CHD
03.08.2016
4.6
26
0
7
5
2
5
0
54
3.4
343
106
CHF
04.08.2016
4.7
1
1
12
3
0
4
0
24
6.3
245
71
SH
04.08.2016
8.5
18
1
1
1
1
15
4
86
1.3
C5
05.08.2016
6.7
7
1
3
1
0
3
0
7
6.0
AV2200
06.08.2016
5.3
0
1
4
1
0
2
0
2Fluids of volcanic basement rocks (outside the magma injection zones) are characterized by gas compositions containing a broad distribution of methane (~70 vol.%, wells R3, K1, GK1, Table3(a)). The most striking example is well R3, which reveals a gas reservoir at the depth interval of 366–455 m, which has a pressure of 24.2 bar. Testing this well showed the production of 3.15 kg/s of water and 4.02 l/s of gas at a discharge temperature of 18°C and depths ranging within 366–1503 m.Table 3
(a)
Gas chemistry (vol%) and carbon isotopic data of deep wells of the Koryaksky-Avachinsky volcanogenic basin, according to sampling data collected from 1988 to 2013. Wells: GK1, K1, and R3 (well positions are shown in Figure1). GK1, G1, data from G. Ryabinin (pers. com., 2011); the remaining samples were collected by A. Kiryukhin and T. Rychkova. Chemical analyses were performed in the Central Chemistry Lab (Guseva V. I.). Isotopic analyses (δ13C) were performed by B. G. Pokrovsky
##
Data
T°C
H2
Ar
O2
N2
CO2
CH4
Sum
δ
13С (СН4)
δ
13С (СО2)
GK1
25.05.1988
0
0.1
0.0
20.9
0.2
78.7
99.9
K1
07.08.2013
0.005
0.1
0.1
22.6
0.7
71.3
94.8
-
34
-
27.7
R3
23.10.2013
0.017
0.1
37.5
0.1
60.7
98.4
-
61
-
49.7
(b)
Gas chemistry (vol%) and carbon isotopic data of hot and cold springs of the Koryaksky-Avachinsky volcanogenic basin, according to sampling data collected from 2011–2016; ##, denoted names, see Table1 for explanation. Sampling was performed by A. Kiryukhin, T. Rychkova, and V. Lavrushin. Chemical analyses were performed in the Central Chemistry Lab (Guseva V. I.). Isotopic analyses (δ13C) were performed by B. G. Pokrovsky
##
Data
T°C
H2
Ar
O2
N2
CO2
CH4
Sum
δ
13С (СН4)
δ
13С (СО2)
K2
09.07.2010
0
0.1
0.2
11.0
88.6
0.04
99.9
K2
25.09.2010
0.001
0.3
2.9
28.8
68.1
0.03
100.1
K2
06.08.2013
0.105
0.2
0.1
23.7
75.6
0.04
99.7
K2
06.08.2013
0.1
27.1
72.4
0.05
99.7
-
40.3
-
8.8
K3
06.08.2013
0.0
62.2
37.7
0.05
100.0
-
23.1
-
9
CH
29.07.2011
0
0.0
0.0
2.5
96.7
0.08
99.3
IS
06.08.2013
46
<0.005
0.0
0.3
2.9
95.2
0.00
98.4
-
36.8
-
6.9
FA
10.08.2014
166
0.686
2.4
0.3
3.4
92.0
0.45
99.2
-
6.8
IS
08.08.2014
0.001
0.1
0.3
4.7
94.7
0.01
99.8
K2
08.08.2014
14
0.005
0.2
0.2
32.9
65.0
0.05
98.5
K2
05.08.2015
14
0.008
0.3
0.7
28.2
70.7
0.03
99.9
FA
07.08.2015
94
0.370
0.4
7.1
30.3
60.8
0.05
99.0
FA
07.08.2015
94
0.400
0.4
6.4
27.4
64.5
0.06
99.0
FA
06.08.2016
94
0.31
0.5
4.54
17.62
76.1
0.03
99.1
-
10.1
K1
03.08.2016
<0.005
0.93
2.09
76.38
19.7
0.05
99.1
-
11.6
K2
03.08.2016
<0.005
0.42
1.54
70.67
27.0
0.03
99.6
-
21
K8
04.08.2016
<0.005
0.55
1.02
20.44
76.8
0.01
98.8
-
50.7
-
18.7
VD
03.08.2016
<0.005
0.24
2.36
12.52
84.5
0.04
99.7
-
22.1
CHD
03.08.2016
<0.005
0.27
1.87
26.57
70.5
0.35
99.5
-
22.3
-
8.1
CHF
04.08.2016
<0.005
0.68
9.5
35.56
48.0
0.35
94.0
-
34
-
6Thermal mineral waters discharged by springs on the foothills of volcanoes differ significantly from waters circulated in basement rocks. Spring waters are characterized by the presence of СО2 in the gas phase (Table 3(b)), lower pH values (Table 2(b)) and lower salinity values (1–~4 g/l); these are related to the significant mixing of these waters with meteoric waters at near-surface conditions. This represents НСО3-Na water with low Cl concentrations that are comparable to the Cl concentrations of adjacent rivers and creeks. Increases in their SO4 concentrations are also observed, especially in the Isotovsky spring (IS). All of the above-mentioned features may be explained as a result of the interaction of meteoric waters enriched by volcanic CO2 gas upflows with volcanogenic rocks. СО2 leaching may produce enrichments in Na and sometimes Са. Sources of SO4 may include either volcanic gases (with CO2) or leaching by meteoric waters from hydrothermally altered rocks in the subcrater zones of volcanoes. The formation of such waters in kinetically active zones is reflected by a wide range of geothermometry values (200–300°C) (Table 2(b)), while basement waters show a narrower range of geothermometry values.The conditions of water recharge in the Koryaksky-Avachinsky volcanogenic basin and thermal mineral springs can be easily identified using the isotopic compositions (δD, δ18O) of water (Figure 5). The light isotopic compositions of the Koryaksky Narzan, Izotovsky, and Pinachevsky-2 springs, as well as the wells of the Bystrinsky groundwater reservoir, correspond to recharge areas at elevations of 2000–2500 masl (which can only represent the slopes of Koryaksky and Avachinsky volcanoes, especially their glacier sites). The water recharge areas of the Chistinsky Narzany group are located at elevations of 800–1500 masl, which correspond to the position of the central part of the Pinachevsky extrusions, with the Arik and Aag volcanoes. The chemical and gas compositions of the fluids of thermal mineral springs in the zone of magmatic injections differ sharply from those of the basement groundwater, as their gas compositions predominantly comprise СО2 (78–97.5%) and record high-temperature Na-K geothermometry values (more than 290–320°С) (Table 2(b)).Figure 5
Isotopic compositions of water (δD and δ18O). (1) Cluster of hot mineral springs (Isotovsky, Koryaksky Narzan, Pinachevsky, total of 36 samples, 2010–2016) and Bystrinsky cold groundwater reservoir (7 samples from 2010–2016); (2) one cold spring at a foothill of Avachinsky Volcano (2200 masl) and Elisovo cold groundwater reservoir (7 samples, 2014–2016); (3) cluster of hot mineral springs (Chistinsky, Vodopadny, total of 4 samples, 2014–2016) and cold springs at a north foothill of Koryaksky Volcano (total of 7 samples, 2010–2016); (4) water phase from Western Kamchatka CH4-water reservoirs (Kshuksky, Nizhne-Kvakchiksky, total of 12 samples, 2010–2016); (5) water condensate from Avachinsky and Koryaksky Volcano fumaroles (14 samples, 1997, 2010–2016); (6) Ketkinsky geothermal field line (5 samples, 2014); (7) Kamchatka representative sites of meteoric water samples (with elevations, masl) and meteoric line. Sampling was performed by A. Kiryukhin (1–7), T. Rychkova (6), Y. Taran (5), N. Malik (5), and P. Voronin (3). Analyses of samples were performed by P. O. Voronin and A. Y. Polyakov using an LGR IWA 45EP isotopic analyzer.Theδ13С carbon isotopic compositions of the СО2 in the free gas samples of the carbonated springs on the north slope of Koryaksky Volcano (K2, K3 and IS) (Table 3(b)) range from −9.0 to -6.9‰, which is within the range recorded by the Koryaksky and Avachinsky fumaroles (see below); they also record a significant magmatic fraction in the CO2 component of the above-mentioned springs [20]. The carbon δ13С isotopic compositions of the СН4 in the free gas samples from the methane wells of K1 (Ketkinsky geothermal field, δ13С = -34‰) and R3 (Radyginskaya area δ13C = -61‰) differ significantly: gases discharged from well K1 formed at higher temperature conditions than the gases penetrated by well R3, where a significant fraction of methane is “marsh methane” of a microbial origin. Interestingly, the wells in the gas condensate fields in western Kamchatka fall within this range from −36.9 to -38‰ (Nizhne-Kvakchikskoye) and -43.1‰ (Kshukskoye) as a result of the magma-hydrothermal gas generation induced in its host volcanogenic-sedimentary rocks. Nevertheless, the low values of δ13С (СО2) (from −49.7 to -16.8‰) observed in methane wells indicate that some CO2 here has a nonmagmatic origin (Table 3(a)) and is paragenetically related to methane sources.According to Taran, gases from the fumaroles of Avachinsky and Koryaksky volcanoes (FA, FK, Figure5) are characterized by a range of isotopic composition reflecting the mixture of meteoric and magmatic (-58‰<δD<-30‰, -0.1‰<δ18O<+7‰, -11.8<δ13С<-5.2) end members. The presence of high methane concentrations (up to 0.3–0.6 mmole/mole) suggests that a considerable amount of regional methane reaches thermal water from the basements of volcanoes. The analyses of our samples of the Avachinsky fumaroles that were collected during 2014–2016 and those of the Ketkinsky geothermal field that were collected in 2014 generally support the above-mentioned mixing model (Figure 5).


## 5. Conceptual Thermal-Hydrodynamic Modeling of the Formation of the Hydrothermal System beneath Koryaksky Volcano

### 5.1. TOUGH2 Model Setup

Thermal-hydrodynamic modeling was used here as a tool to verify a conceptual model of the formation of a hydrothermal reservoir beneath the active Koryaksky Volcano. For this purpose, we applied the TOUGH2 family of codes, which was designed for the modeling of multiphase nonisothermal flows in porous/fractured reservoirs [21]. Here, we used TOUGH2 with different fluid equations of state (EOS) modules, corresponding to the conditions of the geofluids of Koryaksky Volcano: EOS1, water in two phases; EOS1 + tracer, water in two phases plus tracer, which here is considered to be a chloride; and EOS2, water plus CO2 in two phases.The model geometry includes the upper part of the volume of the dike injection beneath Koryaksky Volcano and the surface thermal features on its northern slope (i.e., Koryaksky Narzan and Isotovsky thermal mineral springs). This model corresponds to the vertical cross-sectional AB shown in Figure2. The bottom of the model was restricted to an absolute depth of −4 km, while the top of the model coincides with the top surface of Koryaksky Volcano. A 2D rectangular grid was generated using regular subdivisions in its horizontal and vertical directions, of ΔX = 500 m (NX=49) and ΔZ = 200 m (NZ=35), respectively; thus, the total number of model elements was 1276 (excluding elements above the topographic surface). A model width of ΔY = 1 km (NY=1) was also assigned.The following geological units were represented in the model as domains with different petrophysical properties: (1) basement metamorphic rocks (domain “Base”), (2) Neogene-Pleistocene volcanogenic-sedimentary rocks (domain “NQ1”), (3) Koryaksky volcanic rocks (domain “Volca”), and (4) magma injection units (domain “Dike”) (Figure 6 and Table 4(a)). We also defined model domains to assign different types of near-surface boundary conditions: (5) fixed-state boundary conditions (two-phase P, Sg, and PCO2) on the top surface of Koryaksky Volcano (domain “Fixe1”); and (6) fixed-state boundary conditions (single-phase P, T, PCO2, and Cl) on the foothills of Koryaksky Volcano at an elevation of 300 masl (domain “Fixe2”) (Table 4(b)). Fixed-state boundary conditions reflect unsaturated conditions in Koryaksky Volcano (“Fixe1”) and conditions of adjacent local groundwater aquifer systems (with water level elevations of up to 300 masl) (“Fixe2”). Other initial conditions assigned to the model are shown in Table 4(b).Table 4
(a)
Rock material properties used in the model. Notes: Corey relative permeability values of Slr = 0.3 and Sgr = 0.05 and a specific heat of 1000 J/kg°C were assigned to all model domains. Rock density, porosity, and heat conductivity data are from Chernyak et al. [9], whereas permeability data are from Kiryukhin et al. [10], using analogous reservoirs
Domain parameters
Model domains
Base
NQ1
Volca
Dike
Fixe 1 (Fixed State)
Fixe 2 (Fixed State)
Permeability, 10−15 m2
1
10–100
1000
1000
1000
10–100
Porosity
0.02
0.20
0.2
0.02
0.2
0.20
Density, kg/m3
2700
2200
2200
2700
2200
2200
Heat conductivity, W/m°C
2.0
1.4
1.4
2.0
1.4
1.4
(b)
Initial and boundary conditions assigned in the model
Primary variables
EOS module
Model domains
Base
NQ1
Volca
Dike
Fixe 1 (fixed state)
Fixe 2 (fixed state)
Pressure Pa
EOS1/EOS2
HydrostaticΔP/ΔZ=10 MPa/km
1
e
5
Temperature °C
EOS1/EOS2
Δ
T
Δ
Z
=
24
°
C
/
k
m
Saturation
EOS1/EOS2
0.9
Chloride mass fraction
EOS1 + tracer
19
e
-
3
1
e
-
6
1
e
-
6
CO2 partial pressure, Pa
EOS2
40
(c)
Sources/sinks assigned in the model
Sources/Sinks
Model domains
Base
NQ1
Dike
Heat flux/rate
0.05 W/m2 (lower boundary)
1.0–2.5 MW
Mass rate (CO2)
1
e
-
3 kg/s
Wells on deliverabilityPI=1e-11 m3, Pb=1 bar
Thermal features: FK, K1, IS, VD, CHFigure 6
Model grid, zonation, and boundary conditions (Base, NQ1, Volca, Dike, Fixe 1, Fixe 2, see Table4), and thermal features assigned as “wells on deliverability” (FK, K1, IS, VD, CH, see Table 1). This model position corresponds to the vertical cross-sectional AB shown in Figure 2. Horizontal and vertical axes are labeled in meters.The sources and sinks in the model were assigned using the following method (Table4(c)). Conductive heat flux was defined in the bottom of the model, and constant heat rates were assigned to the model elements corresponding to dike injection zones with a possible grid resolution (Figure 6, Table 4(c)). Constant mass rates of CO2 were assigned to the same model elements (if using module EOS2). These sources were considered to be unknown parameters during the calibration of the model (to match geothermometers for hot mineral springs and CO2 content data). The discharge thermal features were modeled in terms of “wells” using the “well on deliverability” option (Table 4(c)); each “well” (FK, K1, IS, VD, CH) represented in the model corresponds to a thermal feature according to Table 1.The time of the modeling was defined as 7000 years. We believe that the shape of Koryaksky Volcano and the conditions of magma injection did not change significantly during this period. This time also corresponds to the duration of the postglacial period in Kamchatka, when relatively stable meteoric water recharge conditions of Koryaksky Volcano were maintained.


### 5.2. Conceptual TOUGH2 Modeling Results

We performed a number of direct TOUGH2 runs to estimate the rates of the heat sources assigned in the zones of dike injections (“Dike” domain in the model). In the range of 0.1 to 10 MW per model element (the volume of each element is1·108m3, and a total of 85 elements of the “Dike” domain occur in the model), the value of 1 MW was found to be reasonable to explain the temperature evolution of the hydrothermal system below Koryaksky Volcano to match the estimated geothermometry temperatures (253–333°C) in the thermal mineral springs discharged on the northern foothills of Koryaksky Volcano (i.e., Koryaksky Narzan, Isotovsky, see Table 2(b)).Figure7 shows that, after 7000 years of the evolution of the hydrothermal system, the volume of the geothermal reservoir at temperatures above 250°C reached 12.0 km3, with a maximum temperature of 323°C. This reservoir is located at an absolute depth below −0.3 km. The volume of the two-phase subreservoir at the same time is estimated to be 5.2 km3, with its saturation ranging from 0.05 to 0.19 (Figure 8).Figure 7
Distributions of modeled temperatures (°C) (TOUGH2-EOS1) (6A). Horizontal and vertical axes are labeled in meters.Figure 8
Distributions of modeled gas saturation values (Sg) (TOUGH2-EOS1) (6A). Horizontal and vertical axes are labeled in meters.The distributions of fluid pressures after 7000 years show a significant relative pressure minimum of −78 bar at an elevation of −3.9 km at the base of a high-temperature upflow zone (Figure9); that pressure difference pumps down water from the shallow meteoric groundwater systems on the foothills of Koryaksky Volcano (where the fixed-state pressure/temperature boundary is assigned in the “Fixe2” domain of the model, see Figure 6) into the deep geothermal reservoir.Figure 9
Distributions of modeled pressures (bar) (TOUGH2-EOS1) (6A). Horizontal and vertical axes are labeled in meters.The issue of
C
O
2
distribution in the hydrothermal system was analyzed using the TOUGH2 model with a CO2 (EOS2) fluid module [21]. In this model, sources of CO2 were assigned to a magma injection region (Table 4(c)), assuming the magmatic origin of CO2 in the hydrothermal system below Koryaksky Volcano (see Section 3 and Table 3(b)).The TOUGH2-EOS2 output modeling results of the evolution of the hydrothermal system are similar to those obtained above (without CO2) in terms of their temperature, pressure, and saturation distributions, while the CO2 partial pressures (which reach up to 2.4 bar) cause some maximum temperatures to decrease (−7°C) and saturations to increase (+0.08) in the geothermal reservoir. Figure 10 shows the distributions of partial CO2 pressures in this model, which reach values of above 1 bar in regions of dike intrusions (with a maximum value of 2.5 bar), follow high-temperature circulation patterns, and diffuse into cold water inflow regions based on the nonlinear solubility properties of CO2.Figure 10
Distributions of modeled partial CO2 pressures (bar) (TOUGH2-EOS2) (CO2-8). Horizontal and vertical axes are labeled in meters.The distribution of chloride brines in a hydrothermal system was analyzed using the TOUGH2-EOS1+tracer option, where a tracer is assigned as a chloride mass fraction in a fluid. Chloride brines (of up to 13 g/kg) were found in the deep wells of the basement and the Neogene rock units of the Koryaksky-Avachinsky volcanogenic basin. Therefore, in this model, we consider the scenario of an initially seawater-saturated basin (with chloride concentrations of 19 g/kg, or a mass fraction of 0.019) as the initial conditions before magma injection started (Table 4(b)). Figure 11 shows the distribution of the mass fraction of chloride in the model; it also clearly shows a diluted meteoric water circulation pattern (along the 0.0135 isoline) that traces a streamline from a shallow meteoric groundwater system (representing the recharge area for the hydrothermal system) into a region of deep dike intrusions (where heat sources also act as pumps) and then passes to a discharge area zone, where the Koryaksky Narzan and Isotovsky hot springs, K1 and IS, are defined in this model.Figure 11
Distributions of modeled chloride Cl− mass fractions (TOUGH2-EOS1 + tracer) (6B + CL). Horizontal and vertical axes are labeled in meters.


## 5.1. TOUGH2 Model Setup

Thermal-hydrodynamic modeling was used here as a tool to verify a conceptual model of the formation of a hydrothermal reservoir beneath the active Koryaksky Volcano. For this purpose, we applied the TOUGH2 family of codes, which was designed for the modeling of multiphase nonisothermal flows in porous/fractured reservoirs [21]. Here, we used TOUGH2 with different fluid equations of state (EOS) modules, corresponding to the conditions of the geofluids of Koryaksky Volcano: EOS1, water in two phases; EOS1 + tracer, water in two phases plus tracer, which here is considered to be a chloride; and EOS2, water plus CO2 in two phases.The model geometry includes the upper part of the volume of the dike injection beneath Koryaksky Volcano and the surface thermal features on its northern slope (i.e., Koryaksky Narzan and Isotovsky thermal mineral springs). This model corresponds to the vertical cross-sectional AB shown in Figure2. The bottom of the model was restricted to an absolute depth of −4 km, while the top of the model coincides with the top surface of Koryaksky Volcano. A 2D rectangular grid was generated using regular subdivisions in its horizontal and vertical directions, of ΔX = 500 m (NX=49) and ΔZ = 200 m (NZ=35), respectively; thus, the total number of model elements was 1276 (excluding elements above the topographic surface). A model width of ΔY = 1 km (NY=1) was also assigned.The following geological units were represented in the model as domains with different petrophysical properties: (1) basement metamorphic rocks (domain “Base”), (2) Neogene-Pleistocene volcanogenic-sedimentary rocks (domain “NQ1”), (3) Koryaksky volcanic rocks (domain “Volca”), and (4) magma injection units (domain “Dike”) (Figure 6 and Table 4(a)). We also defined model domains to assign different types of near-surface boundary conditions: (5) fixed-state boundary conditions (two-phase P, Sg, and PCO2) on the top surface of Koryaksky Volcano (domain “Fixe1”); and (6) fixed-state boundary conditions (single-phase P, T, PCO2, and Cl) on the foothills of Koryaksky Volcano at an elevation of 300 masl (domain “Fixe2”) (Table 4(b)). Fixed-state boundary conditions reflect unsaturated conditions in Koryaksky Volcano (“Fixe1”) and conditions of adjacent local groundwater aquifer systems (with water level elevations of up to 300 masl) (“Fixe2”). Other initial conditions assigned to the model are shown in Table 4(b).Table 4
(a)
Rock material properties used in the model. Notes: Corey relative permeability values of Slr = 0.3 and Sgr = 0.05 and a specific heat of 1000 J/kg°C were assigned to all model domains. Rock density, porosity, and heat conductivity data are from Chernyak et al. [9], whereas permeability data are from Kiryukhin et al. [10], using analogous reservoirs
Domain parameters
Model domains
Base
NQ1
Volca
Dike
Fixe 1 (Fixed State)
Fixe 2 (Fixed State)
Permeability, 10−15 m2
1
10–100
1000
1000
1000
10–100
Porosity
0.02
0.20
0.2
0.02
0.2
0.20
Density, kg/m3
2700
2200
2200
2700
2200
2200
Heat conductivity, W/m°C
2.0
1.4
1.4
2.0
1.4
1.4
(b)
Initial and boundary conditions assigned in the model
Primary variables
EOS module
Model domains
Base
NQ1
Volca
Dike
Fixe 1 (fixed state)
Fixe 2 (fixed state)
Pressure Pa
EOS1/EOS2
HydrostaticΔP/ΔZ=10 MPa/km
1
e
5
Temperature °C
EOS1/EOS2
Δ
T
Δ
Z
=
24
°
C
/
k
m
Saturation
EOS1/EOS2
0.9
Chloride mass fraction
EOS1 + tracer
19
e
-
3
1
e
-
6
1
e
-
6
CO2 partial pressure, Pa
EOS2
40
(c)
Sources/sinks assigned in the model
Sources/Sinks
Model domains
Base
NQ1
Dike
Heat flux/rate
0.05 W/m2 (lower boundary)
1.0–2.5 MW
Mass rate (CO2)
1
e
-
3 kg/s
Wells on deliverabilityPI=1e-11 m3, Pb=1 bar
Thermal features: FK, K1, IS, VD, CHFigure 6
Model grid, zonation, and boundary conditions (Base, NQ1, Volca, Dike, Fixe 1, Fixe 2, see Table4), and thermal features assigned as “wells on deliverability” (FK, K1, IS, VD, CH, see Table 1). This model position corresponds to the vertical cross-sectional AB shown in Figure 2. Horizontal and vertical axes are labeled in meters.The sources and sinks in the model were assigned using the following method (Table4(c)). Conductive heat flux was defined in the bottom of the model, and constant heat rates were assigned to the model elements corresponding to dike injection zones with a possible grid resolution (Figure 6, Table 4(c)). Constant mass rates of CO2 were assigned to the same model elements (if using module EOS2). These sources were considered to be unknown parameters during the calibration of the model (to match geothermometers for hot mineral springs and CO2 content data). The discharge thermal features were modeled in terms of “wells” using the “well on deliverability” option (Table 4(c)); each “well” (FK, K1, IS, VD, CH) represented in the model corresponds to a thermal feature according to Table 1.The time of the modeling was defined as 7000 years. We believe that the shape of Koryaksky Volcano and the conditions of magma injection did not change significantly during this period. This time also corresponds to the duration of the postglacial period in Kamchatka, when relatively stable meteoric water recharge conditions of Koryaksky Volcano were maintained.


## 5.2. Conceptual TOUGH2 Modeling Results

We performed a number of direct TOUGH2 runs to estimate the rates of the heat sources assigned in the zones of dike injections (“Dike” domain in the model). In the range of 0.1 to 10 MW per model element (the volume of each element is1·108m3, and a total of 85 elements of the “Dike” domain occur in the model), the value of 1 MW was found to be reasonable to explain the temperature evolution of the hydrothermal system below Koryaksky Volcano to match the estimated geothermometry temperatures (253–333°C) in the thermal mineral springs discharged on the northern foothills of Koryaksky Volcano (i.e., Koryaksky Narzan, Isotovsky, see Table 2(b)).Figure7 shows that, after 7000 years of the evolution of the hydrothermal system, the volume of the geothermal reservoir at temperatures above 250°C reached 12.0 km3, with a maximum temperature of 323°C. This reservoir is located at an absolute depth below −0.3 km. The volume of the two-phase subreservoir at the same time is estimated to be 5.2 km3, with its saturation ranging from 0.05 to 0.19 (Figure 8).Figure 7
Distributions of modeled temperatures (°C) (TOUGH2-EOS1) (6A). Horizontal and vertical axes are labeled in meters.Figure 8
Distributions of modeled gas saturation values (Sg) (TOUGH2-EOS1) (6A). Horizontal and vertical axes are labeled in meters.The distributions of fluid pressures after 7000 years show a significant relative pressure minimum of −78 bar at an elevation of −3.9 km at the base of a high-temperature upflow zone (Figure9); that pressure difference pumps down water from the shallow meteoric groundwater systems on the foothills of Koryaksky Volcano (where the fixed-state pressure/temperature boundary is assigned in the “Fixe2” domain of the model, see Figure 6) into the deep geothermal reservoir.Figure 9
Distributions of modeled pressures (bar) (TOUGH2-EOS1) (6A). Horizontal and vertical axes are labeled in meters.The issue of
C
O
2
distribution in the hydrothermal system was analyzed using the TOUGH2 model with a CO2 (EOS2) fluid module [21]. In this model, sources of CO2 were assigned to a magma injection region (Table 4(c)), assuming the magmatic origin of CO2 in the hydrothermal system below Koryaksky Volcano (see Section 3 and Table 3(b)).The TOUGH2-EOS2 output modeling results of the evolution of the hydrothermal system are similar to those obtained above (without CO2) in terms of their temperature, pressure, and saturation distributions, while the CO2 partial pressures (which reach up to 2.4 bar) cause some maximum temperatures to decrease (−7°C) and saturations to increase (+0.08) in the geothermal reservoir. Figure 10 shows the distributions of partial CO2 pressures in this model, which reach values of above 1 bar in regions of dike intrusions (with a maximum value of 2.5 bar), follow high-temperature circulation patterns, and diffuse into cold water inflow regions based on the nonlinear solubility properties of CO2.Figure 10
Distributions of modeled partial CO2 pressures (bar) (TOUGH2-EOS2) (CO2-8). Horizontal and vertical axes are labeled in meters.The distribution of chloride brines in a hydrothermal system was analyzed using the TOUGH2-EOS1+tracer option, where a tracer is assigned as a chloride mass fraction in a fluid. Chloride brines (of up to 13 g/kg) were found in the deep wells of the basement and the Neogene rock units of the Koryaksky-Avachinsky volcanogenic basin. Therefore, in this model, we consider the scenario of an initially seawater-saturated basin (with chloride concentrations of 19 g/kg, or a mass fraction of 0.019) as the initial conditions before magma injection started (Table 4(b)). Figure 11 shows the distribution of the mass fraction of chloride in the model; it also clearly shows a diluted meteoric water circulation pattern (along the 0.0135 isoline) that traces a streamline from a shallow meteoric groundwater system (representing the recharge area for the hydrothermal system) into a region of deep dike intrusions (where heat sources also act as pumps) and then passes to a discharge area zone, where the Koryaksky Narzan and Isotovsky hot springs, K1 and IS, are defined in this model.Figure 11
Distributions of modeled chloride Cl− mass fractions (TOUGH2-EOS1 + tracer) (6B + CL). Horizontal and vertical axes are labeled in meters.


## 6. Discussion

In this section, the results of the model are matched to observations in order to improve the model; we also discuss the ability of the model to reproduce the temperature response of hot springs to recent dike injections. Finally, the powers of the modeled heat sources are compared to that of Koryaksky Volcano in terms of the heat content of its erupted lava.


### 6.1. Asymmetric Discharge of Hot Springs and Distribution of Centripetal Pressures

The distribution of the temperatures output by the TOUGH2 modeling (Figure7) is rather symmetrical in shape, but all of the hot springs adjacent to Koryaksky Volcano are located in an area that is 7 to 12 km to the north (Figure 1). Another observed misfit between the modeling results and observations is a centripetal trend of deep well levels rising toward the centers of the Koryaksky-Avachinsky volcanic cluster [7, Figure 7], while the distribution of TOUGH2 modeled pressures (Figure 9) shows a decrease in pressure toward the center of the volcano at an elevation of −1000 masl (Figure 8).A reasonable way to adjust the model to the above-mentioned observations is to assume that cold water injection occurs beneath Koryaksky Volcano, which shifts the thermal anomaly in the northern direction and builds up pressure in the bottom of the volcanic structure. To do this, a domain defined as MPipe was added to the model (Figure12). MPipe represents the funnel of Koryaksky Volcano in the model, where cold water may recharge the underlying volcanic structure. One element of the MPipe domain was assigned in fixed-state conditions (at a pressure of 1 bar and a temperature of 10°C) to define the water level position in the volcano funnel at 700 masl. We also increased the heat rates of the sources up to 2.5 MW (Table 4(c)).Figure 12
Revised model grid, zonation, and boundary conditions (compare to Figure6) (6A + I). Domain MPipe (Koryaksky Volcano funnel) has been added, and the blue rectangle inside this domain corresponds to the water level in the volcano funnel. Horizontal and vertical axes are labeled in meters.The temperature distribution in this case was reshaped from a symmetrical (Figure7) to an asymmetrical (Figure 13) form, with explicit thermal outflow occurring in the northern direction (K1, IS, VD, CH thermal features) and a hidden geothermal reservoir beneath the southern foothills of Koryaksky Volcano. The volume of the geothermal reservoir at a temperature above 250°C reached 8.5 km3, with a maximum temperature of 339°C. This reservoir is located at an absolute depth below −1.0 km. Figure 14 shows that pressures now have a centripetal distribution, with a relative pressure maximum of +12.1 bar at an elevation of −3.9 km at the base of the high-temperature upflow zone. Fluid circulation also switched from a free convection mode (Figure 9) to a forced convection mode (Figure 14). This modeling scenario also shows a significant increase of 93 bar in the fluid pressure at the base of the high-temperature upflow zone (Figure 14 compared to Figure 9), while single-phase conditions prevailed beneath the structure of Koryaksky Volcano.Figure 13
Distributions of modeled temperatures (°C) (TOUGH2-EOS1) (6A + I). Horizontal and vertical axes are labeled in meters.Figure 14
Distributions of modeled pressures (bar) (TOUGH2-EOS1) (6A + I). Horizontal and vertical axes are labeled in meters.Assigning to the TOUGH2-EOS2 model recharge rates of CO2 sources in an injection area (Table 4(c)) generates two separate CO2-enriched reservoirs with partial CO2 pressures of up to 1.5–2.0 bar (with CO2 concentrations of 30–40 g/kg) (Figure 15). The larger reservoir is shallower and is aligned beneath the known CO2 hot springs (i.e., K1, IS, VD, CH).Figure 15
Distributions of modeled CO2 partial pressures (bar) (TOUGH2-EOS2) (CO2-8 + 2I). Horizontal and vertical axes are labeled in meters.The analysis of the redistribution of chloride brines due to hydrothermal circulation was performed using the TOUGH2-EOS1 + tracer option, where the tracer is assigned as a chloride mass fraction in the fluid. In the case of volcanic funnel-type water recharge conditions, nearly the entire basin is purged and diluted of chloride within a few thousand years. Thus, it is most feasible that the chloride component observed in the thermal features on the northern slope of Koryaksky Volcano has a magmatic origin.


### 6.2. Temperature Response of the Isotovsky Hot Spring to Dike Injections

Observational temperature data in hot mineral springs located on the northern slope of Koryaksky Volcano represent an additional possible data source for the verification of thermal-hydrodynamic models. The Isotovsky hot spring is located 7 km from the summit of Koryaksky Volcano (Figure1) and is characterized by flowing temperatures of up to 50°C. The temperature monitoring of this spring was performed during the time period of 2010–2016 using the HOBO U12 temperature logger (with a 15 min−1 frequency of records) (Figure 16). During the observational period, four plane-oriented earthquake clusters (produced by dike injections) were identified in this local area (##163, 165, 194, and 204; these numbers are shown in Figure 1).Figure 16
Temperatures recorded in Isotovsky hot spring and suggested times of dike injections in adjacent areas (#194, 2.08.2011 and #204, 28.02.2016).T, observational data; Tav, annual maximum monthly temperature; Tmax, maximum monthly temperature.Figure16 shows that a postmagmatic temperature increase of 6–12°С occurred during winter, from October 2011 to June 2012, after event #194 (which is interpreted to represent a dike injection that occurred on 02.08.2011). At this time, the maximum monthly temperature Tmax increase of 6–12°С is compared to the annual maximum monthly temperature Tav. In other words, the temperature perturbation caused by the suggested dike injection was recorded 2 months after the dike injection occurred; then, this temperature increase was maintained over a 10-month time period. Another anomalous temperature increase in Isotovsky of up to 10°С was associated with a dike injection that occurred on 28.02.2016.Modeling (using the model parameters defined in Section4) shows that a distance of ≈1 km (between the heat source and the observational point) represents the upper distance limit of such thermal responses. This is close to the value of 2 km observed in the field (Figure 1). Nevertheless, it is difficult to establish definite conclusions about this, as a coarse grid is used in the model and the seismic event coordinate data may not be accurate. The fracture porosity (or fluid active circulation volume) may also significantly increase the distance of the propagation of a temperature anomaly from a heat source (dike) to the point of heat discharge (hot spring). Additional modeling shows that the thermal response of an increase of 11°C in a few months may be achieved at a distance of 2 km, if the fluid active circulation volume is defined as 0.001 (volume fraction) in the model and the heat sources in the dike injection region are assigned as 8 MW (Figure 17).Figure 17
Thermal response (temperature rise, °C) in a model element corresponding to Isotovsky hot spring after additional heat sources (4 × 8 MW) were assigned in the model 2 km apart. Note: active circulation volume (volume fraction) used as a model sensitivity parameter here. (TOUGH2-EOS1) (6A + I+).


### 6.3. Geothermal Budget of Koryaksky Volcano

The relationship between volcanic, hydrothermal, and seismic activity is another issue that the analysis of the Koryaksky-Avachinsky volcanic cluster can be used to assess. Using the identified geothermal resource data from the world’s eight largest geothermal-producing countries and the allocation of active volcanoes around the world, Stefansson [22] suggested that there is a linear relationship between the number of active volcanoes and the geothermal potential that can be used for the production of electricity, which statistically means that one active volcano has an average potential close to that of a hydrothermal system, which is capable of producing 158 MWe of geothermal power. Nevertheless, it is still not clear why, how, and where active volcanoes store their energy when they are not erupting magma onto the Earth’s surface.The TOUGH2-modeling of Koryaksky Volcano yields an example of this mechanism. A total of 85 model elements were assigned heat sources with a rate of 2 MW; thus, the total heat power used to run a model hydrothermal system was 170 MW, and this value should be doubled to 340 MW (since we need to at least double the width of the 2D model from 1 km to 2 km to cover the area of dike injections and surface thermal features, Figure1). Koryaksky Volcano was characterized by an average magma rate of 150 kg/s during its lifetime of 45 kY [23]. This corresponds to an average heat power of Koryaksky Volcano of 150 MW (as 1 kg/s of magma is approximately converted to a heat rate of 1 MW). That means that double this heat power (and magma rate) is needed to run the hydrothermal system beneath Koryaksky Volcano.In the model was specified cold water recharge through the volcanic funnel by using a constant fluid pressure in model element #1149 at an elevation of 700 masl, which yields a water recharge of 289 kg/s; thus, the total cold water recharge required to run a real hydrothermal system should be doubled to 578 kg/s (since we need to at least double the width of the 2D model from 1 km to 2 km to cover the area of dike injections and surface thermal features, Figure1). The area of the structure of Koryaksky Volcano at an elevation above +2000 masl is estimated to be 18.5 km2, most of which is covered by glaciers and permanent snow fields. Assuming that the annual atmospheric water precipitation rate at these elevations is more than 5000 mm [24], this may produce a water recharge rate of more than 3 m3/с. If only 20% of this flow is converted into the underground drain, it would be sufficient to provide water recharge to its underlying hydrothermal systems.


## 6.1. Asymmetric Discharge of Hot Springs and Distribution of Centripetal Pressures

The distribution of the temperatures output by the TOUGH2 modeling (Figure7) is rather symmetrical in shape, but all of the hot springs adjacent to Koryaksky Volcano are located in an area that is 7 to 12 km to the north (Figure 1). Another observed misfit between the modeling results and observations is a centripetal trend of deep well levels rising toward the centers of the Koryaksky-Avachinsky volcanic cluster [7, Figure 7], while the distribution of TOUGH2 modeled pressures (Figure 9) shows a decrease in pressure toward the center of the volcano at an elevation of −1000 masl (Figure 8).A reasonable way to adjust the model to the above-mentioned observations is to assume that cold water injection occurs beneath Koryaksky Volcano, which shifts the thermal anomaly in the northern direction and builds up pressure in the bottom of the volcanic structure. To do this, a domain defined as MPipe was added to the model (Figure12). MPipe represents the funnel of Koryaksky Volcano in the model, where cold water may recharge the underlying volcanic structure. One element of the MPipe domain was assigned in fixed-state conditions (at a pressure of 1 bar and a temperature of 10°C) to define the water level position in the volcano funnel at 700 masl. We also increased the heat rates of the sources up to 2.5 MW (Table 4(c)).Figure 12
Revised model grid, zonation, and boundary conditions (compare to Figure6) (6A + I). Domain MPipe (Koryaksky Volcano funnel) has been added, and the blue rectangle inside this domain corresponds to the water level in the volcano funnel. Horizontal and vertical axes are labeled in meters.The temperature distribution in this case was reshaped from a symmetrical (Figure7) to an asymmetrical (Figure 13) form, with explicit thermal outflow occurring in the northern direction (K1, IS, VD, CH thermal features) and a hidden geothermal reservoir beneath the southern foothills of Koryaksky Volcano. The volume of the geothermal reservoir at a temperature above 250°C reached 8.5 km3, with a maximum temperature of 339°C. This reservoir is located at an absolute depth below −1.0 km. Figure 14 shows that pressures now have a centripetal distribution, with a relative pressure maximum of +12.1 bar at an elevation of −3.9 km at the base of the high-temperature upflow zone. Fluid circulation also switched from a free convection mode (Figure 9) to a forced convection mode (Figure 14). This modeling scenario also shows a significant increase of 93 bar in the fluid pressure at the base of the high-temperature upflow zone (Figure 14 compared to Figure 9), while single-phase conditions prevailed beneath the structure of Koryaksky Volcano.Figure 13
Distributions of modeled temperatures (°C) (TOUGH2-EOS1) (6A + I). Horizontal and vertical axes are labeled in meters.Figure 14
Distributions of modeled pressures (bar) (TOUGH2-EOS1) (6A + I). Horizontal and vertical axes are labeled in meters.Assigning to the TOUGH2-EOS2 model recharge rates of CO2 sources in an injection area (Table 4(c)) generates two separate CO2-enriched reservoirs with partial CO2 pressures of up to 1.5–2.0 bar (with CO2 concentrations of 30–40 g/kg) (Figure 15). The larger reservoir is shallower and is aligned beneath the known CO2 hot springs (i.e., K1, IS, VD, CH).Figure 15
Distributions of modeled CO2 partial pressures (bar) (TOUGH2-EOS2) (CO2-8 + 2I). Horizontal and vertical axes are labeled in meters.The analysis of the redistribution of chloride brines due to hydrothermal circulation was performed using the TOUGH2-EOS1 + tracer option, where the tracer is assigned as a chloride mass fraction in the fluid. In the case of volcanic funnel-type water recharge conditions, nearly the entire basin is purged and diluted of chloride within a few thousand years. Thus, it is most feasible that the chloride component observed in the thermal features on the northern slope of Koryaksky Volcano has a magmatic origin.


## 6.2. Temperature Response of the Isotovsky Hot Spring to Dike Injections

Observational temperature data in hot mineral springs located on the northern slope of Koryaksky Volcano represent an additional possible data source for the verification of thermal-hydrodynamic models. The Isotovsky hot spring is located 7 km from the summit of Koryaksky Volcano (Figure1) and is characterized by flowing temperatures of up to 50°C. The temperature monitoring of this spring was performed during the time period of 2010–2016 using the HOBO U12 temperature logger (with a 15 min−1 frequency of records) (Figure 16). During the observational period, four plane-oriented earthquake clusters (produced by dike injections) were identified in this local area (##163, 165, 194, and 204; these numbers are shown in Figure 1).Figure 16
Temperatures recorded in Isotovsky hot spring and suggested times of dike injections in adjacent areas (#194, 2.08.2011 and #204, 28.02.2016).T, observational data; Tav, annual maximum monthly temperature; Tmax, maximum monthly temperature.Figure16 shows that a postmagmatic temperature increase of 6–12°С occurred during winter, from October 2011 to June 2012, after event #194 (which is interpreted to represent a dike injection that occurred on 02.08.2011). At this time, the maximum monthly temperature Tmax increase of 6–12°С is compared to the annual maximum monthly temperature Tav. In other words, the temperature perturbation caused by the suggested dike injection was recorded 2 months after the dike injection occurred; then, this temperature increase was maintained over a 10-month time period. Another anomalous temperature increase in Isotovsky of up to 10°С was associated with a dike injection that occurred on 28.02.2016.Modeling (using the model parameters defined in Section4) shows that a distance of ≈1 km (between the heat source and the observational point) represents the upper distance limit of such thermal responses. This is close to the value of 2 km observed in the field (Figure 1). Nevertheless, it is difficult to establish definite conclusions about this, as a coarse grid is used in the model and the seismic event coordinate data may not be accurate. The fracture porosity (or fluid active circulation volume) may also significantly increase the distance of the propagation of a temperature anomaly from a heat source (dike) to the point of heat discharge (hot spring). Additional modeling shows that the thermal response of an increase of 11°C in a few months may be achieved at a distance of 2 km, if the fluid active circulation volume is defined as 0.001 (volume fraction) in the model and the heat sources in the dike injection region are assigned as 8 MW (Figure 17).Figure 17
Thermal response (temperature rise, °C) in a model element corresponding to Isotovsky hot spring after additional heat sources (4 × 8 MW) were assigned in the model 2 km apart. Note: active circulation volume (volume fraction) used as a model sensitivity parameter here. (TOUGH2-EOS1) (6A + I+).


## 6.3. Geothermal Budget of Koryaksky Volcano

The relationship between volcanic, hydrothermal, and seismic activity is another issue that the analysis of the Koryaksky-Avachinsky volcanic cluster can be used to assess. Using the identified geothermal resource data from the world’s eight largest geothermal-producing countries and the allocation of active volcanoes around the world, Stefansson [22] suggested that there is a linear relationship between the number of active volcanoes and the geothermal potential that can be used for the production of electricity, which statistically means that one active volcano has an average potential close to that of a hydrothermal system, which is capable of producing 158 MWe of geothermal power. Nevertheless, it is still not clear why, how, and where active volcanoes store their energy when they are not erupting magma onto the Earth’s surface.The TOUGH2-modeling of Koryaksky Volcano yields an example of this mechanism. A total of 85 model elements were assigned heat sources with a rate of 2 MW; thus, the total heat power used to run a model hydrothermal system was 170 MW, and this value should be doubled to 340 MW (since we need to at least double the width of the 2D model from 1 km to 2 km to cover the area of dike injections and surface thermal features, Figure1). Koryaksky Volcano was characterized by an average magma rate of 150 kg/s during its lifetime of 45 kY [23]. This corresponds to an average heat power of Koryaksky Volcano of 150 MW (as 1 kg/s of magma is approximately converted to a heat rate of 1 MW). That means that double this heat power (and magma rate) is needed to run the hydrothermal system beneath Koryaksky Volcano.In the model was specified cold water recharge through the volcanic funnel by using a constant fluid pressure in model element #1149 at an elevation of 700 masl, which yields a water recharge of 289 kg/s; thus, the total cold water recharge required to run a real hydrothermal system should be doubled to 578 kg/s (since we need to at least double the width of the 2D model from 1 km to 2 km to cover the area of dike injections and surface thermal features, Figure1). The area of the structure of Koryaksky Volcano at an elevation above +2000 masl is estimated to be 18.5 km2, most of which is covered by glaciers and permanent snow fields. Assuming that the annual atmospheric water precipitation rate at these elevations is more than 5000 mm [24], this may produce a water recharge rate of more than 3 m3/с. If only 20% of this flow is converted into the underground drain, it would be sufficient to provide water recharge to its underlying hydrothermal systems.


## 7. Conclusions

(1) The fluid recharge/discharge conditions of the Koryaksky-Avachinsky volcanogenic basin were studied using water (δD and δ18O) and δ13С (gas) isotopic data. The isotopic compositions (δD, δ18O) of the thermal mineral springs (Koryaksky Narzans, Isotovsky, and Pinachevsky) and wells of the Bystrinsky and Elizovo aquifers indicated that water recharge occurred from +2000 to +2500 masl at glaciers on the summits of Koryaksky and Avachinsky volcanoes. Chistinsky Narzans was fed from the central part of the Pinachevsky extrusions (Aric and Aag volcanoes) at elevations of 800–1500 masl. The δ13С isotopic compositions of СО2 sampled from Koryaksky Narzan and Isotovky revealed its magmatic origin. Thus, CO2 springs in the northwestern foothills of Koryaksky Volcano formed as a result of the mixing of magmatic gases and melting glacial waters. The hydrothermal reservoir beneath the northern slope of Koryksky Volcano records geothermometric temperatures ranging from 253 to 333°C and CO2 contents of up to 3 g/kg (Koryaksky Narzan).(2) Plane-oriented clusters of seismic events (KB GS RAS data of 2000–2016) below Koryaksky and Avachinsky volcanoes were used to identify magma injection zones: (1) the shallow crustal magma chamber (i.e., a combination of sills and dikes) below the southwestern base of Koryaksky Volcano at absolute depths ranging from −2 to −5 km, with a width of 2.5 km; (2) dike accumulation in a nearly north-south zone (7.5 by 2.5 km, with absolute depths ranging from −2 to −5 km) in the northern sector of Koryaksky Volcano; and (3) the shallow magma chamber (i.e., a combination of sills and dikes) in a cone of Avachinsky Volcano at absolute depths ranging from 1 to 2 km, with a width of 1 km.(3) Conceptual TOUGH2 modeling was used to understand and explain the mechanism of the formation of the hydrothermal system beneath Koryaksky Volcano (Figure 18). For this purpose, the following terms were found to be crucial in this model: (1) heat sources of 20 MW/km3 and gas (CO2) sources of 10 g/s/km3 acting during 7000 years in the above-defined zones of magma injections; and (2) cold water recharge of 580 kg/s through the volcanic funnel to the deep dike injection area. The modeling results reasonably match the Na-K geotemperature estimates of geothermal reservoirs (300°C), the isotopic values (δD, δ18O) of high-level meteoric water recharge, the concentrations of magmatic CO2 (up to 4 g/kg) in the hot springs on the northern slope of Koryaksky Volcano, the thermal reaction to the 02.08.2011 dike injection recorded in Isotovsky hot spring, and the fluid dilution of the original seawater basin beneath the volcano due to cold water recharge circulation. This modeling also indicates that a hidden high-temperature geothermal reservoir is present beneath the southern slope of Koryaksky Volcano (at an elevation of −1 km), which may become a subject of future drilling explorations.Figure 18
Conceptual TOUGH2-based model of hydrothermal circulation beneath Koryaksky Volcano. 1: temperatures, °C; 2; pressures, bar; 3; regions with CO2 partial pressures above 1 bar (i.e., CO2 mass fractions more than 0.02); 4: regions with chloride concentrations of more than 5 ppm; 5: fluid flows: a, cold, b, hot, and c, intermediate; 6: dike injection regions. For additional explanation of this figure, see legend in Figure 2.



---

*Source: 4279652-2017-10-22.xml*