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.
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 Appendix A 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: () using seismic data to identify dike swarms and magma chambers (i.e., potential heat sources); () 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; () 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): () many crustal magma chambers are produced from sills; () active crustal magma chambers inject magma into associated volcanoes, with most chambers existing in a partially molten (pore-like) state; () ellipsoids are thermally stable shapes of magma chambers; () 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); () the excess pressure during eruptions decreases in an exponential manner, until the dike loses its hydraulic connection to the chamber; () the pressure of magma in the chamber during a period of repose must equal the overall lithostatic pressure; () 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 () 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.
The basin basement comprises Upper Cretaceous , 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 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 (Figure 1). 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: () a time difference ( = 1 day); () a distance difference in the horizontal plane ( < 6 km); and () a requirement of a nearly planar orientation (i.e., a distance from the event to the plane that is less than 200 m). When the resulting cluster contains more than five elements , 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 element contains the coordinates (, , ). For (the number of elements in the cluster) points with coordinates (, , ), we can find the equation of the fitting plane using the least-squares method. The solution thus reduces to solving a set of linear equations as follows:
These equations are then solved using Cramer’s rule. In this way, we obtain the coefficients , , for the equation of a plane, which is defined as . The next step is to find the unit vector that is normal to the fitting plane , where is the determinant of the equations that result in the following geological parameters: dip angle and the azimuth of dip . 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 (Figure 1) revealed the following geomechanical features (Figures 1–4): () 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; () 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); () 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); () 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).
4. Koryaksky-Avachinsky Basin Discharge/Recharge Conditions Based on Gas, Chemical, and Isotopic Data
Intervals of deep sampling wells (Table 2(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.
Fluids 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, Table 3(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.
Thermal 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)).
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 , 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С = ) and R3 (Radyginskaya area δ13C = ) 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 (Nizhne-Kvakchikskoye) and (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 ) 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, Figure 5) are characterized by a range of isotopic composition reflecting the mixture of meteoric and magmatic (, , ) 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 Figure 2. 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 = 500 m () and = 200 m (), respectively; thus, the total number of model elements was 1276 (excluding elements above the topographic surface). A model width of = 1 km () was also assigned.
The following geological units were represented in the model as domains with different petrophysical properties: () 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 ) on the top surface of Koryaksky Volcano (domain “Fixe1”); and (6) fixed-state boundary conditions (single-phase P, T, , 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).
The sources and sinks in the model were assigned using the following method (Table 4(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 is , 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)).
Figure 7 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).
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 (Figure 9); 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.
The issue of 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.
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.
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 (Figure 7) 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 ], 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 (Figure 12). 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)).
The temperature distribution in this case was reshaped from a symmetrical (Figure 7) 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.
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).
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 (Figure 1) 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 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 increase of 6–12°С is compared to the annual maximum monthly temperature . 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 Section 4) 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).
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, Figure 1). 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, Figure 1). 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
() 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).
() 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: () 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.
() 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: () 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.
Appendix
A. Photos of the Thermal Features of the Koryaksky-Avachinsky Volcanic Cluster (Kamchatka)
B. Analysis of the Sensitivity of the Frac-Digger Algorithm to the Criteria of Temporal and Spatial Proximity
This section demonstrates the feasibility of using the Frac-Digger program to extract plane-oriented clusters of seismic events (dikes, fracs) from seismic data catalogs within a reasonable range of a base set of temporal and spatial proximity parameters. The input data used for this test include 5160 seismic events collected by KB GS RAS during the time period from Jan. 2000 to July 2016 using the seismic station network around the Koryaksky-Avachinsky volcanic group. We ran Frac-Digger by varying a base set of parameters, which are responsible for the selection of the plane-oriented clusters: () , cluster size, , 8, 10; (2) the distance () from the event to the plane, = 150, 200, 300 m; and (3) the time interval () allocated to the cluster, = 1, 7, 30 days.
Figures 31, 32, and 33 show that if the cluster size varies from 6 to 10, we obtain () estimated average depths of dike injections remain at elevations ranging from −4000 to −2500 masl (Koryaksky) and from 1000 to 2000 masl (Avachinsky) (Figure 31); (2) estimated average dip angles of dikes remain in the range of 59.18 to 59.44° (Figure 32); and (3) the estimated dip azimuths of dikes show two dominant directions at 90° (east) and 270° (west) (Figure 33).
Figures 34, 35, and 36 show that if the distance () from the event to the plane varies from 150 to 300 m, we obtain the following: () estimated average depths of dike injections remain at elevations ranging from −4000 to −3000 masl (Koryaksky) and from 1000 to 2000 masl (Avachinsky) (Figure 34); (2) estimated average dip angles of dikes remain in the range of 58.12 to 59.35° (Figure 35); and (3) the estimated dip azimuths of dikes show two dominant directions at 90° (east) and 270–290° (west) (Figure 36).
Figures 37, 38, and 39 show that if the time interval () allocated to the cluster varies from 1 to 30 days, we obtain the following: () estimated average depths of dike injections remain at elevations from −4000 to −3000 masl (Koryaksky) and from 500 to 2000 masl (Avachinsky) (Figure 37); (2) estimated average dip angles of dikes remain in the range of 59.35 to 63.64° (Figure 38); and (3) the estimated dip azimuths of dikes show two dominant directions at 70–90° (east) and 270–290° (west) (Figure 39).
Finally, Figure 40 shows the Frac-Digger sensitivity of estimated dike traces at a depth of −3000 masl below Koryaksky Volcano determined using a range of the above-mentioned base sets of temporal and spatial proximity parameters. It can be visually seen that the main trends of dike orientation are maintained, despite the fact that the parameters of cluster selection have been significantly changed. This indicates the reasonably low sensitivity of the Frac-Digger algorithm to spatial criteria (, proximity or distance from the event to the plane; , the number of events included in a cluster) and temporal proximity criteria (, the time interval allocated for earthquake events included in a cluster), as is demonstrated above.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was funded by the Russia Science Foundation Project no. 16-17-10008.