Magma Fracking Beneath Active Volcanoes Based on Seismic Data and Hydrothermal Activity Observations

: Active volcanoes are associated with microearthquake (MEQ) hypocenters that form plane-oriented cluster distributions. These are faults delineating a magma injection system of dykes and sills. The Frac-Digger program was used to track fracking faults in the Kamchatka active volcanic belt and fore-arc region of Russia. In the case of magma laterally injected from volcanoes into adjacent structures, high-temperature hydrothermal systems arise, for example at Mutnovsky and Koryaksky volcanoes. Thermal features adjacent to these active volcanoes respond to magma injection events by degassing CO 2 and by transient temperature changes. Geysers created by CO 2 gaslift activity in silicic volcanism areas also flag magma and CO 2 recharge and redistributions, for example at the Uzon-Geyserny, Kamchatka, Russia and Yellowstone, USA magma hydrothermal systems. Seismogenic faults in the Kamchatka fore-arc region are indicators of geofluid fracking; those faults can be traced down to 250 km depth, which is within the subduction slab below primary magma sources.


Introduction
Volcanoes and crustal magma chambers are the products of magma injection from primary magma chambers at depths of 150-200 km below the Kuril-Kamchatka volcanic zone ( Figure 1). These magma injections create new hydraulic fractures and reactivate existing faults. The formation of shear cracks in adjacent to the main aseismic opening triggers microseismicity [1]. Earthquake magnitudes (M) for shear deformations with amplitudes of 0.1 mm-1 cm and fracture lengths of a few hundred meters are estimated to range from 1 to 2 (Ks from 3.5 to 5.5), corresponding to the sensitivity of local seismic networks operating in the areas of the Klyuchevskoy, Koryaksky-Avachinsky, and Mutnovsky-Gorely volcano groups. We propose that the planes defined by clusters of microearthquakes correspond to magmatic intrusions. For a search of plane-oriented clusters of earthquakes, we use our original programs Frac-Digger (RU reg. ## 2016616880), and for production feed-zones, we use Frac-Digger2 (RU reg. ## 2017618050) (these applications discrete fracture systems using seismic and geological data). In these programs, the following criteria are used to select clusters of earthquake hypocenters or production feed zones: (1) proximity in time δt (Frac-Digger only); (2) proximity within the horizontal plane δR; (3) proximity to plane orientation δZ (distance between the object and plane); and (4) minimum number of elements in cluster N. This method has already been tested on the example of the 2012 eruption of Tolbachik volcano [2], the 2000-2016 MEQ data of Koryaksky-Avachinsky volcanic cluster [3][4][5], Paratunsky graben (adjacent to Vilyuchinsky volcano) [6], the Mutnovsky-Gorely volcanic cluster [7], and the Klyuchevkoy group of volcanoes [8]. Note that the methodology of dyke tracking using the Frac-Digger program is presented in detail in Section 3.2.1, "Seismic Data and Method of Plane-Oriented Clusters Identification" [5].
By coupling magma fracking observations beneath volcanoes with the monitoring of thermal features in adjacent hydrothermal systems, we can develop appropriate geofiltrational and thermalhydrodynamic TOUGH2-models that allow an assessment of recoverable geothermal resources and explain the mechanisms of the anomalous hydrothermal perturbations. A recently revealed mechanism [9] of geysers activity due to cyclic magmatic CO2-gaslift pointed to the possibility of understanding the mechanisms of relationship between anomalous hydrothermal perturbations, volcano activity, and strong earthquakes. Continuous measurements of CO2 partial pressures in the Koryaksky Narzan thermal mineral springs (northern sector of Koryaksky volcano, 2017-2019) and in the Mutnovsky geothermal Power Plant condenser (2019) provide evidence of this process.

Koryaksky-Avachinsky Volcanoes
Plane-oriented clusters of 5160 seismic events from data of the Kamchatka Branch Federal Research Center United Geophysical Survey Russia Academy of Sciences (KB FRC UGS RAS) of 01.2000-07.2016) below Koryaksky and Avachinsky volcanoes were used to identify magma injection zones [4,5]. Cluster identification was carried out using our Frac-Digger program. The previously defined criteria were used to include a new event in the cluster: δt ≤ 1 day, δR ≤ 6 km, δZ ≤ 0.2 km, N ≥ 6.
It was found that 30% of the total number of 5160 earthquakes formed 204 plane-oriented clusters. A careful analysis of the orientations of dykes and sills (using stereograms and histograms) showed the following: (1) the dykes beneath Koryaksky Volcano were mostly emplaced in the depth range from −5000 to 0 m masl; they mostly strike north-south (75% of dykes have strike azimuths between 320° and 40°) and have dip angles over 50° (70% of the dykes), and (2) the dykes beneath Avacha Volcano were mostly in the depth range 1000-2000 m masl. Their strikes are distributed over all possible directions, but there is a local maximum of dykes striking nearly north-south with 25% of dykes having azimuths between 350° and 10°. Steeply dipping dykes are not obviously dominant (61% have dip angles over 45°).
Additionally: (1) A shallow crustal magma chamber, likely as a combination of sills and dykes, lies below the southwestern base of Koryaksky Volcano at depths ranging from -2 to -5 km masl and a width of 2.5 km; (2) dyke accumulation in a nearly north-south zone (7.5 by 2.5 km, with absolute depths ranging from -2 to -5 km occurs beneath the northern sector of Koryaksky Volcano; and (3) a shallow magma chamber in the cone of Avachinsky Volcano lies at absolute depths ranging from 1 to 2 km and has a width of 1 km.
Thus, we propose the geomechanical state according to [10]: (1) extension is recorded beneath Koryaksky Volcano producing normal faulting (NF conditions), the vertical stress Sv is the maximum stress, the horizontal principal stresses are coincident with the north-south (SHmax) and the eastwest (Shmin) directions, (2) the cone of Avacha Volcano is found to be under conditions of radial extension.
The relationship M = 0.67 log (V) − 0.82 between the upper bound of the induced seismic event M and magmatic injection volume V, suggested by [11], may be used for recorded dykes volume estimates beneath Koryaksky in a range from 220 m 3 to 1.988 × 10 6 m 3 , and the total recorded volume of the injected magma is estimated at 5.862 × 10 6 m 3 during 2008-2011. Avachinsky volcano characterized by a range from 111 m 3 to 0.064 × 10 6 m 3 , and the total recorded volume of the injected magma is estimated at 0.169 × 10 6 m 3 during 2008-2011.
It is also worth noting that seismogenic planes and microearthquake (MEQ) mechanisms planes (19 mechanisms estimated were used) are intersected with the angle of 51° (on average) [12]. This may have the following geomechanical interpretation: seismogenic planes mark opening mode fractures (hydrofracking-type fractures) formed orthogonal to the least effective stress during magma injection, whereas MEQ mechanisms planes are mark shear faults, which initiated from the main opening mode fault (dyke) using pre-existing fractures systems.
There were six magma injections below the northern foothills of Koryaksky volcano detected in 2011-2019 based on MEQ analysis data using Frac-Digger criteria (δt = 1 day, δR = 6 km, δZ = 0.2 km, N = 6): one dyke on 02.08.2011, one dyke on 28.02.2016, one dyke on 22.08.2017, two dykes on 14.11.2017, and one dyke on 22.04.2019. All of these dykes were injected just below the Koryaksky Narzan and Isotovsky thermal mineral springs area. The depth of magma injections was estimated from seismic data from −7000 to +1000 masl, with M ranges from 0.85 to 2.95 ( Figure 2, Table 1).  Table 1 below. Note: Dykes patches in 3D view are convex polygons, with vertices in the points of projections of plane-oriented clusters of earthquakes hypocenters on an approximation plane. If a less severe Frac-Digger criterion in time, δt ≤ 30 day, is used, the number of detected dykes significantly increases. By using this approach, it was found that as a result, 91% of the total number of 5160 earthquakes 452 plane-oriented clusters were formed. Figure 3 shows that 68 of these features were interpreted as dykes from 2011 to 2019 beneath Koryaksky volcano. Most of them were injected in a northeast sector of the volcano structure, where thermal springs activity occurs.

Mutnovsky and Gorely Volcanoes
We already used plane-oriented earthquake clusters (KB FRC UGS RAS, catalogs 2009-2017) to track the dykes injected beneath and around the Mutnovsky volcano by using the Frac-Digger program [5]. Thus, magmatic injection zones (shape-like dykes and sills, [13,14]) defined in this manner were treated as pathways and heat sources for ascending high-temperature upflows that feed adjacent hydrogeological reservoirs and surface feature (hot springs and fumaroles) discharge.
Four seismic stations record seismicity in the Mutnovsky-Gorely volcanic cluster. A total of 1336 earthquakes were recorded by the KB FRC UGS RAS in the edifices and basement of the Mutnovsky and Gorely volcanoes between January 2009 and February 2017. Cluster identification was carried out using our Frac-Digger program. The same criteria as above were used to include a new event in the cluster. Our treatment of these data used the method previously outlined, which yielded 973 earthquakes (72.8%) that make up 74 plane-oriented clusters between November 2013 and October 2016 (most of them associated with the Mutnovsky volcano) ( Figure 4).
The analysis of seismic activity at the Mutnovsky volcano revealed the following geomechanical features.  Thus, dyke geometry indicated reverse fault (RF) geomechanic conditions [10] in the vicinity of the Mutnovsky volcano (SHmax > Shmin > Sv, SHmax striking to NE), while there is a trend to normal fault (NF) conditions 10 km away in the region where the production geothermal reservoir formed. The relationship M = 0.67 log (V) − 0.82 between the upper bound of the induced seismic event M and magmatic injection volume V, suggested by [11], may be used for recorded dyke volume estimates in a range from 520 m 3 to 0.252 × 10 6 m 3 , and the total recorded volume of the injected magma is estimated at 0.67 × 10 6 m 3 from 2009-2016.
There were additional 22 magma injections identified during 02.2017 to 10.2019 in the Mutnovsky-Gorely geothermal area. Most of the dykes beneath the Mutnovsky volcano have a dip angle ranging from 20° to 40° and were injected at the depth ranging from −4.0 to −2.0 km masl and dipping in the SEE-SE (120°-130°) direction. Most of the dykes beneath the Gorely volcano have a dip angle ranging from 40° to 70° and were injected at the depth ranging from −2.0 to −1.5 km masl and dipping in the NWW (280°-320°) direction. It is worth noting two sub-parallel dykes (dip angles from 39° to 43°, dip azimuth from 149° to 153°, depth range from −2.6 to −1.9 km masl, M from 1.25 to 1.7) injected in the vicinity of well 035 in 2019 ( Figure 4).

Northern Group of Volcanoes
To recover the sequence of magma injections in the area of the Klyuchevskoy Volcanic Cluster, we used a catalog containing data of seismic monitoring at 19 seismic stations operated by the KB FRC UGS RAS for the 2000-2017 period (the total number of earthquakes that have been recorded between January 1, 2000 and August 23, 2017 is 122,451). The method that we used to identify planeoriented earthquake clusters and fitting planes using the Frac-Digger software was described in [2][3][4][5]. The relevant calculated plane patches are interpreted as zones of magma emplacement in the form of dykes and sills. The calculations assumed the following parameters for the identification of planeoriented clusters: δt ≤ 1 month, δR ≤ 6 km, δZ ≤ 1 km, N ≥ 6. The calculations yielded 1788 clusters that contained 117,412 earthquakes (96% of the total number of recorded events).
A vertical profile along the axis of the Klyuchevskoy Volcanic Cluster that extends from the Tolbachik volcanoes to Shiveluch ( Figure 5) shows that when extrapolated downward, the dykes beneath Tolbachik and Shiveluch intersect at a point that lies between absolute depths of -165 and -205 km beneath Klyuchevskoy; it is at this location that the region of primary magma melting is thought to reside, or there is a primary magma chamber to provide magma for all active volcanoes in the Klyuchevskoy Cluster (Klyuchevskoy, Plosky Tolbachik, Bezymyanny, and Shiveluch). The 2000-2017 magma injections directly beneath Klyuchevskoy were concentrated in the depth ranges between −31 and −28 km masl (35%) and between −1 and +2 km masl (20%), where we hypothesize the existence of a crustal chamber (C2) and a peripheral chamber (C1), respectively (we mean chamber as a plexus of sills and dykes here). The crustal magma chamber (C2) was found to contain dykes that dip at angles of 70-80° (30%) and sills that dip at 15-20° (8%); one also finds an increasing number of nearly east-west striking dykes that supply magma to the nearby Bezymyanny and Krestovsky volcanoes. Nine episodes of intensive dyke emplacement and seven episodes of intensive sill emplacement were recorded in the C2 chamber during the period of interest (2000-2017); the dyke episodes preceded the sill episodes by an advance time of one year (two cases), the episodes occurred simultaneously (five cases), or they did not terminate in sill emplacement (one case), which provides evidence of a change in the geomechanical condition around the C2 chamber in the range from horizontal extension to radial compression.
The dip angles and azimuths of dykes and sills in the peripheral chamber (C1) are distributed rather uniformly; one notes some increase in the number of dykes that dip at angles of 35° to 50° and 65° to 80°. There is also a set of dykes that dip at northeast azimuths (i.e., striking along the Klyuchevskoy-Krestovsky line).
Tolbachik Volcano is characterized by injections of magma in the depth range between −8 and −1 km masl, with most dykes dipping more steeply than 60° (67%). In addition, we identified a set of dykes that dip at azimuths of 180-220° (26%) and strike along the Ostry Tolbachik to Plosky Tolbachik to the Bol'shaya Udina line.
Shiveluch is characterized by magma injections in the depth ranges between −4 and −2 km masl (37%) and between 0 and +2 km (47%), where two peripheral magma chambers are hypothesized to reside. The injections mostly occurred in the form of sills that dip at angles below 5°, with some trend of prevailing western dip azimuths.
Bezymyanny Volcano receives magma from the crustal magma chamber beneath Klyuchevskoy Volcano (C2), which lies at a depth between −31 and −28 km masl. The injections beneath Bezymyanny occurred in the depth range between −2 and +2 km (93%), where a peripheral magma chamber is hypothesized to exist. The injections mostly occurred in the form of sills that dip at angles below 15° (27%) and low-angle dykes that dip at 50° or less (58%); the dip azimuths are distributed uniformly over all directions.
The volumes of magma injection from magma chambers can be estimated using empirical relationships that relate the maximum magnitude of the triggered seismicity to the volume of injected magma [11]: M = 0.67 log(V) − 0.82, where M is the maximum magnitude of triggered seismic events and V is the injection volume. The results from this assessment of dyke and sill injections volumes for the volcanoes in the Klyuchevskoy Cluster during 2000-2017 are Klyuchevskoy (peripheral chamber C1) − 0.6 × 10 6 m 3 , Klyuchevskoy (crustal chamber C2) − 2.4 × 10 6 m 3 , Bezymyanny − 4.7 × 10 6 m 3 , Tolbachik -13.2 × 10 6 m 3 , Shiveluch -22 × 10 6 m 3 . Assuming the enthalpy of magma at a temperature of 1200 °С to be 1000 kJ/kg and the density of magma to be 2800 kg/m 3 , one can estimate the mass discharge and the corresponding thermal power due to magma injections for the time elapsed from January 2000 to August 2017.
The thermal power values for magma injections in the volcanic plumbing systems considered here and the output estimates for the respective volcanoes from [15] yields the following ratios between the discharge of magmas stored beneath volcanoes and that of magma ejected onto the ground surface (volume of intruded magma)/(volume of erupted magma): 0.8% for Klyuchevskoy, 14.9% for Bezymyannyy, 23.2% for Tolbachik, and 72.9% for Shiveluch.
Thus, Kluchevskoy volcano magma injections (dykes) took place in "normal fault" (NF) conditions, forming a permeable reservoir down to −35 km msl and two magma chambers within it; Shiveluch volcano magma injections (sills) occurs in "reverse fault" (RF) conditions at shallow levels from −4 km masl to −2 km masl, forming sills in an area 15 km across.

Kamchatka East Volcanic Belt and Adjacent Shelf Area
Critically stressed faults are key players in creation of the productive reservoirs and generation of strong earthquakes [10]. We performed a Frac-Digger analysis of the Kamchatka's regional seismicity to identify the network of such active faults, which are identified as plane-oriented clusters of hypocenters in Frac-Digger terms [16].
We used the regional catalog of earthquakes from KB FRC UGS RAS, which includes 5972 earthquakes with M values above 4.25 during the time period from 01.1980 to 02.2016. The following criteria for seismogenic planes selection were used: N ≥ 6, δt ≤ ∞, δZ ≤ 10 km, δR ≤ 100 km. These criteria correspond to the assumption of the existing network of continuously active regional faults.
Based on the above, we found 156 plane-oriented clusters of earthquakes, which are interpreted as active seismogenic faults. Most of them are located beneath Kamchatka's eastern shelf, between the shore line and ocean trench ( Figure 6). Seventeen faults are found to be the most active (with more than 100 earthquakes each). Most of the active faults are characterized by dip angles from 50° to 70° and a dip azimuth NWW from 300° to 310°, striking subparallel to the ocean trench. This points to the extension in the NWW direction and NF geomechanical conditions as a whole.
Seismogenic fault #110 includes the hypocenter of earthquake M = 7.3 on 24 Nov 1971, which was the strongest felt in Petropavlovsk-Kamchatsky since 1959. Faults with other directions (especially in shallow conditions at elevations above −10 km masl) are found too, which reflects local geomechanical features. One of such faults is #4 (175 earthquakes included, dip angle 53°, dip azimuth 217°) striking in the direction of Petropavlovsk-Kamchatsky.
Seismogenic fault #35 includes the hypocenter of earthquake M = 6.9 on 17 Aug. 1983, which was accompanied with the ground deformations described by [17]. Seismogenic fault #35 has a dip angle of 60.3°, dip azimuth of 281.5°, and includes 77 hypocenters of earthquakes. Seismogenic faults distributions point to a lower limit of the hydrofrack propagation into the subduction plate (down to −250 km masl), where this fluid can act as the key ingredient for magma melting to recharge the primary chambers of active volcanoes (Figure 7).
It is also worth noting that 98% of all earthquakes formed a plane-oriented clusters network, which characterizes the geomechanical conditions of collided plates and suggests hydrofracking mechanisms due to fluids (possible phases are: water, oil, gas) generation there.
In case of water, the empirical relationship M = 0.67 log (V) + 1.42 between the upper bound of the induced seismic event M and water injection volume V, suggested by [11], may be used for water volumes recharge estimates into the active faults network of the Kamchatka shelf. This yields in a range from 0.00007 to 0.3 km 3 , and the total volume of the recharged water is estimated at 1.87 km 3 during the time interval from 1980 to 2016 (over approximately 750 km of arc length).
In case of gas, the empirical relationship is slightly modified (supercritical CO2 injection) to M = 0.67 log (V) − 0.30 [11]; then, this may be used for rough estimates of gas volumes generated from deep hydrocarbon sources to be injected into the active faults network of the Kamchatka shelf. This yields in a range from 0.024 to 111 km 3 , and the total volume of the injected gas is estimated at 691 km 3 during the time interval from 1980 to 2016. In relation to this, it is worth noting that methane hydrates and methane submarine discharges are widely distributed along the east coast of Kamchatka.
Then, we perform an additional analysis of the 3D distribution of the 200 strongest earthquakes (М > 5.65) from the above-mentioned catalog of seismic events of KB FRC UGS RAS. The following selection criteria were used in Frac-Digger2: δZ = 4 km, δR = 100 km, N = 5. The output results show that 102 of the 200 strongest earthquakes form 11 plane-oriented clusters of hypocenters. A comparative analysis of seismogenic faults orientation with the mechanisms of corresponding earthquakes (estimated at http://www.globalcmt.org/CMTsearch.html) shows an average intersection angle of 33°. This may have the following geomechanical interpretation: seismogenic planes are marked as opening mode fractures (hydrofractures), which formed orthogonal to the least effective stress during magma injection, while the EQ's mechanisms planes are marked as shear faults that initiated from the main opening mode fault (dyke) using pre-existing fractures systems.  Figure 1), showing traces of active seismogenic faults. The blue color covers the area where plane-oriented seismogenic faults are possible, which points to hydrofrack or brittle rock conditions. Numbers correspond to the fault numbers shown in Figure 6. The star denotes suggested active volcanoes' primary magma chamber positions.

Uzon-Geysernaya Caldera Geysers
Geysers are examples of cyclically erupted boiling hot springs. It has been shown that the driving mechanism of geysers cycling is gas-lift assisted eruptions, where non-condensable gases (mainly CO2) are important players [9,18] due to CO2 significantly drops boiling temperatures.
There is evidence from gas sampling at the Velikan Geyser that at the time of full activity before Jan. 3, 2014, gas composition was dominated by CO2. Gas sampling from September 2013 showed that the gas component of the hydrothermal reservoir that fed the Velikan Geyser was dominated by carbon dioxide (СО2, 61.5%) and nitrogen (N2, 32.1%), along with а significant amount of methane (CH4, 5.8%) and hydrogen (Н2, 0.45%) [9]. Gas composition in 2014-2019 in the Velikan and Bolshoy geysers reveals nitrogen becoming a dominant gas, while CO2 declined to less than 1% [19]. In contrast, a new geyser, Shaman (Mutny), developed in 2008 in a channel of the former hot springs of Uzon caldera, just 12 km apart [20]. The gas composition of the newly formed Shaman Geyser was characterized by CO2 domination according to our gas sampling in 2018 and 2019.
We interpret this as a redistribution of the magmatic gas recharge from magma plumbing systems of the Uzon-Geysernaya caldera in the following way: the Valley of Geysers' hydrothermal system magmatic CO2 recharge was reduced, while that of the Uzon geothermal reservoir CO2 recharge increased.

Koryaksky Narzan Thermal Springs
The Koryaksky Narzan CO2-reach thermal springs (12-14 °C) and Isotovsky thermal springs   (Figure 9). The most significant drop observed is 0.4 °C, which may be caused by a H2O-CO2 boiling temperature drop in a spring pool due to the magmatic CO2 release associated with dykes injections in the adjacent area.

Mutnovsky Production Reservoir
The Mutnovsky (Dachny) geothermal reservoir includes at least two production faults geometrically and hydraulically connected to the magma fracking system of Mutnovsky volcano [7] (see also Figure 4).
We performed 30-day continuous monitoring of the non-condensable gases (most of gas content is CO2) partial pressure in the turbine condenser of the geothermal power plant. To estimate PCO2, we performed simultaneous measurements of the steam condensate pressure Pc and temperature Tc; then, PCO2 was calculated as a difference between Pc and saturation pressure corresponding to temperature Tc. Figure 10 shows the transient PCO2 change during the observational period of time from 25.08.2019 to 25.09.2019. It is clearly seen that at least 14 maxima of PCO2 synchronized with 14 minimums of Tc, which detect non-condensable gas arrivals into the turbine from the production geothermal reservoir. Some of these PCO2 peaks may be related to magmatic gas recharge impulses, followed by the magma fracking processes described above (see Section 2.2 of this paper).

Discussion
The key question in this paper is: How can we verify that plane-nested earthquakes are tracking magma injections in the shapes of dykes or sills beneath active volcanoes? The answer is in the form of a question, too: what other fluid except magma can demonstrate such hydrofracking capabilities? Superheated water and non-condensable gases (possibly CO2) are other candidates for such working fluid duties. However, superheated water in shallow permeable fracture reservoirs conditions is very sensitive to host rock temperatures and is easily converted into high-compressible two-phase conditions, forming geothermal production fields. Another fluid is CO2, especially magmatic CO2 having less compressibility as compared to superheated steam, and it may act as a working fluid coupling with magma. If so, then we should extend our term of magma fracking beneath active volcanoes to the term magma+CO2 fracking beneath active volcanoes. CO2 impact is also useful to explain traces of fracks in a slope of Koryaksky volcano (Figure 3), with no associated magma discharge on the surface at the same time. In this connection, it is worth noting that magmatic CO2 redistribution appears to be the key reason to explain the transfer of geysers activity from Geysers Valley to Uzon in Kamchatka and recent (2018) reactivation of Steamboat Geyser in Yellowstone.
Another important point is the accuracy of seismic hypocenter data and uniqueness of the Frac-Digger method for plane-oriented shapes definition. There is no uniqueness, since we found more fracks using less severe criteria of selections in Frac-Digger (see Section 2.1). Thus, we suggest that 3D distributions of magma fracks be considered a plausible scenario of magma fracking beneath active volcanoes, but these may include more or less magma+CO2 dykes, depending on the Frac-Digger selection parameters. Nevertheless, some Frac-Digger selections pointed to 90%-95% of earthquake hypocenters belonging to plane-oriented clusters, meaning that fracking is a dominant process there.
High-temperature (HT) geothermal system formation due to magma fracking is also well explained in terms of the magma thermal-hydrodynamic modeling.
Conceptual 2D iTOUGH2-EOS1sc thermal hydrodynamic modeling of the Mutnovsky magmatichydrothermal system [7] reasonably explains its evolution over the most recent 1500-5000 years in terms of heat recharge (supplied by injected dykes from the active funnel Mutnovsky-4) and mass recharge (water injected through the dormant volcanic funnels Mutnovsky-3 and possibly Mutnovsky-2) conditions. We emphasize that the magmatic injection rate is approximately equivalent to the heat discharge rate (455 MWt). This is equivalent to approximately 455 kg/s of magma, which is in turn equivalent to from 15.3 to 25.6 km 3 of magma over 3000-5000 years. This volume value is comparable to the available space for dyke accommodation under regional horizontal extension conditions [23].
Conceptual TOUGH2 modeling was used to understand and explain the mechanism of the formation of the hydrothermal system beneath northern slope of Koryaksky Volcano [4,5]. For this purpose, the following terms were found to be crucial in this model: (1) heat sources of 20 MW/km 3 (340 MW total in 17 km 3 of reservoir rocks) and gas (CO2) sources of 10 g/s/km 3 acting during 7000 years in the zones of magma injections; and (2) cold water recharge of 580 kg/s through the volcanic funnel to the deep dyke injection area. The modeling results reasonably match the Na-K geotemperature estimates of geothermal reservoirs (300 °C), the isotopic values (δD, δ 18 O) of highelevation meteoric water recharge, the concentrations of magmatic CO2 (up to 4 g/kg) in the hot springs on the northern slope of Koryaksky Volcano, and the thermal reactions to the dyke injections recorded in Isotovsky and Koryaksky Narzan hot springs. This modeling also indicates that a hidden high-temperature geothermal reservoir is present also beneath the southern slope of Koryaksky Volcano (at an elevation of −1 km), which may become a subject of future drilling explorations.
Another interesting related issue is where magma came from to recharge volcanoes for magma fracking and consequent eruptions. The possible answer is that water release from the subduction plate (due to hydro-fracking in the plate itself) forms upflows into the mantle edge, which creates melting conditions for primary magma chambers at the depth of 150 km (Figure 7). This closed loop of water-magma-water circulation in the subduction zone may be responsible for strong earthquakes, when hydrofracking drives the opening of shear faults. A mini-loop of such watermagma interplay may be Karymsky Volcano, which has been working almost continuously since 1996 due to water recharge fed from Karymsky Lake into dip seismogenic faults #120 and #122 and from there into the magma chambers of the volcano ( Figure 6).

Conclusions
(1) Active volcanoes are injectors of magma and water into adjacent structures, which creates high-temperature production reservoirs.
(2) Seismic data reveal magma hydrofracking and stress conditions around active volcanoes.
(4) Seismogenic faults on the Kamchatka shelf are indicators of geofluid generation and water propagation to 250 km depth.
(5) Geysers result from CO2-gaslift in active silicic volcanic areas. Continuously performed measurements of CO2 partial pressures in hydrothermal features are a possible key to understanding the mechanisms of the relationship between anomalous hydrothermal perturbations [24], volcanic activity, and strong earthquakes.