Contention between supply of hydrothermal fluid and conduit obstruction: inferences from numerical simulations

Abstract

We investigate a volcanic hydrothermal system using numerical simulations, focusing on change in crater temperature. Both increases and decreases in crater temperature have been observed before phreatic eruptions. We follow the system’s response for up to a decade after hydrothermal fluid flux from the deep part of the system is increased and permeability is reduced at a certain depth in a conduit. Our numerical simulations demonstrate that: (1) changes in crater temperature are controlled by the magnitude of the increase in hydrothermal fluid flux and the degree of permeability reduction; (2) significant increases in hydrothermal flux with decreases in permeability induce substantial pressure changes in shallow depths in the edifice and decreases in crater temperature; (3) the location of maximum pressure change differs between the mechanisms. The results of this study imply that it is difficult to predict eruptions by crater temperature change alone. One should be as wary of large eruptions when crater temperature decreases as when crater temperature increases. It is possible to clarify the implications of changes in crater temperature with simultaneous observation of ground deformation.

Introduction

In recent years, precursors of phreatic eruptions have been captured by observations near the crater (Barberi et al. 1992; Rouwet et al. 2014) and are diverse, including ground deformation and changes in gravity, total magnetic field, gas composition, heat flux, and crater temperature (Barberi et al. 1992; Rouwet et al. 2014). Measurement of crater temperature is relatively easy and has been done for a long time. Changes in crater temperature preceding phreatic eruption have been documented for several eruptions (Minakami 1939; Ossaka et al. 1997; Ehara 2007; Christenson et al. 2010; Strehlow et al. 2017). Thermal observation using drones and remote sensing has become feasible in recent years (Harvey et al. 2016; Mori et al. 2016), and these methods make acquisition of crater temperature easier. In this research, we focus on change in crater temperature as a precursor of a phreatic eruption.

An increase in crater temperature (often the temperature of a crater lake) has been observed as a precursor of phreatic eruption in many cases, including almost all events at Ruapehu during 1992–2012 (Christenson et al. 2010; Strehlow et al. 2017), Meakan in 1988 and 2008 (Japan Meteorological Agency 2012a), Azuma in 1977 (Japan Meteorological Agency 2012b), and most events at Kusatsu-Shirane in 1938, 1982, and 1983 (Minakami 1939; Ossaka et al. 1997). In such cases, an increase in the supply of hydrothermal fluid (or heat) from depth is often proposed as the primary mechanism causing phreatic eruption (Rouwet et al. 2014). However, phreatic eruptions after decreases in crater temperature have also been reported in some cases, such as Ruapehu in 1988, 2006, and 2007 (Christenson et al. 2010), Meakan in 1996 (Japan Meteorological Agency 2012a), some events of Kusatsu-Shirane in 1938, 1982, and 1983 (Minakami 1939; Ossaka et al. 1997), Tokachidake in 2004 (Takahashi et al. 2017). The mechanism leading to phreatic eruption after decreases in crater temperature has been proposed to be pressure increase in the shallow part of the edifice due to permeability decrease in the shallow part of the edifice, perhaps within a conduit that feeds hydrothermal discharge (Christenson et al. 2010; Rouwet et al. 2014; Strehlow et al. 2017). A permeability decrease in the shallow part of such a conduit could restrict flow of heat and mass to the crater and cause a decrease in crater temperature (Christenson et al. 2010; Tanaka et al. 2017). Either physical or chemical processes may cause obstruction of the conduit. Candidate physical processes include mechanical reworking of surface deposits (e.g., D’Oriano et al. 2014) and collapse of the conduit wall (e.g., Calvari et al. 2016). Candidate chemical processes include deposition of native sulfur and hydrothermal minerals in the conduit.

Many previous studies have compared numerical modeling of the hydrothermal system to the volcanic monitoring record. Most of the numerical studies that have treated this particular problem invoked an increase in hydrothermal fluid (or heat) supply from the deep part of the system (Todesco et al. 2010; Fournier and Chardot 2012; Currenti et al. 2017). However, relatively few studies have examined the influence of permeability changes in the shallow part of the conduit (Tanaka et al. 2017). No numerical modeling studies have discussed how both the crater temperature and the surrounding edifice react to a competition between increased hydrothermal flux from the depth and reduced permeability in the shallow part of the conduit.

The goal of the study is to examine how crater temperature and other observations at an active volcano can be used for prediction of a phreatic eruption.

Numerical simulation

We used the ‘STAR’ (Pritchett 1995) as the numerical code to describe hydrothermal fluid circulation by finite difference method, with ‘BRNGAS’ as the equation of state (Pritchett 1995). These tools enabled us to calculate heat and mass flow in a system with three pore components: water (H₂O), an incondensable gas (air in the present study), and a soluble salt (NaCl in the present study), and three pore phases: liquid, gas, and solid precipitate. The ‘BRNGAS’ is valid for pressures to 200 MPa and over the temperature range 0–350 °C. In this study, we consider the system in which degassing magma is supposed to be located far below the water table at around sea level (more exactly, it depends on the permeability of the edifice) and hydrothermal fluid ascends intensively through the central conduit. The calculation region was axisymmetric 2D to represent a simplified conical edifice (Fig. 1). Temperature and pressure were held constant at the ground surface (at 20 °C and 1.013 × 10⁵ Pa with incoming air) and along the vertical boundary on the downstream side (hydrostatic conditions and 30 °C km⁻¹). Thermally insulating and hydraulically impermeable conditions were imposed at the bottom of the model. Meteoric recharge at the land surface was injected at a constant rate (17 cm yr⁻¹ at 20 °C), and a constant heat flow of 6.9 × 10⁻² W m⁻² (c.a. the mean continental heat flux) was supplied at the base of the model. The modeled edifice had uniform permeability (5.0 × 10⁻¹³ m²), with the exception of a high permeability conduit (1.0 × 10⁻¹⁰ m²) near the symmetry axis. The permeability of volcanic rocks ranges widely (10⁻¹⁸–10⁻¹⁰ m², Saar and Manga 1999; Sruoga et al. 2004). We adopted a single value in this research, but the influence of differences in host rock permeability should be considered in future work. The anisotropy of permeability (difference between horizontal permeability and vertical permeability) varies greatly, but in this research, for simplicity, it was assumed that there is no anisotropy of permeability and that the porosity is uniform (0.1). The rock properties are listed in Table 1.

Fig. 1

Two-dimensional, axisymmetric computational domain. The domain extends radially to 5225 m and reaches a maximum depth of 2275 m. The red column is the conduit. A permeability reduction at PCB (permeability changeable block) is imposed at the time of conduit obstruction

Table 1 Rock properties

To support the fumarolic activity at the crater, hydrothermal fluid was injected at the bottom of the permeable central conduit at fixed enthalpy (~1345 kJ/kg corresponding to liquid H₂O at ~300 °C and 10 MPa) and a rate of 1000 kg/s. This injection rate was chosen to reproduce the initial heat discharge rate through the block corresponding to the vent (radius of 50 m; the smallest cell size in the model) being roughly 100 MW, which is a representative value for active fumarolic activities at volcano summit. The rest of the convective heat input (~1250 MW) spreads to the edifice and discharges on the flanks of the volcano and through the open distal boundary.

In subsequent simulations, the initial condition was obtained by simulating a long-lasting (thousand year) injection of hydrothermal fluids, representing near-steady-state conditions. The initial condition showed a wide high-temperature region (100–300 °C) and a thick unsaturated zone (Fig. 2), with a water table located at 1100-m model height. The temperature at the crater (the center of the shallowest node (25 m depth) of the conduit near the symmetry axis) was about 135 °C, which indicates superheated steam discharge at the crater.

Fig. 2

Initial near-steady distribution of a temperature, b pressure, and c gas volume fraction

We investigated changes in conditions (temperature and pore pressure) caused by an increase in hydrothermal fluid flux from depth and/or a reduction in permeability at shallow depth in the conduit (Table 2). The hydrothermal flux was increased stepwise at t = 0 by 0.5 and 1.0 times the initial value. Permeability was reduced at the permeability changeable block (PCB) in the conduit. Christenson et al. (2010) performed simulations for a simple one-dimensional case to estimate the permeability reduction due to deposition of hydrothermal minerals. Their simulations, in which volcanic gas was injected into a porous zone saturated with lake water, revealed that significant deposition of quartz, sulfur, and clay minerals could occur over a short timescale of about 10 days. They also pointed out that deposition of sulfur reduced permeability from 10⁻¹² to 10⁻¹⁷ m². Thus, it is conceivable that the permeability in the conduit could decrease by several orders of magnitude on a short timescale. We decreased the permeability of PCB by 1–3 orders of magnitude at t = 0. After introducing such changes, we observed the condition of the hydrothermal system for 10 years of model time.

Table 2 Hydrothermal flux and permeability of PCB in each run

Results

Figure 3 shows the changes in temperature and pore pressure 1 year after increasing the hydrothermal fluid flux from the deep part of the system (Run 10). Heat accumulation and pressurization occur throughout the conduit. However, large changes in temperature and pore pressure occur only near the conduit. Although pressurization is observed below the water table in the entire edifice, the amplitude of this pressurization is less than that near the conduit by one order of magnitude. Figure 4 shows the temporal changes in temperature and pore pressure at shallow depth near the conduit (observation point in Fig. 1) and in crater temperature over a 10-year period. Both temperature and pore pressure increase over the first several years and then remain high.

Fig. 3

Distribution of change in a temperature and b pressure 1 year after doubling the flux of hydrothermal fluid from the deep part of the system (Run 10, Mʹ = 2000 kg s⁻¹, kʹ_PCB = 1.0 × 10⁻¹⁰ m²). Arrow indicates the location where maximum pressure change occurs in 10-year period

Fig. 4

Temporal change in a temperature at the crater (red line) and the observation point (black line) and b pressure at the observation point after doubling the flux of hydrothermal fluid (Run 10, Mʹ = 2000 kg s⁻¹, kʹ_PCB = 1.0 × 10⁻¹⁰ m²)

Only reductions in permeability at PCB (Run 4) induce heating and pressurization below the PCB accompanied by cooling and depressurization above it (Fig. 5). The region with observed temperature change is limited to the vicinity of the conduit above the water table. Pressurization is induced near the conduit below the PCB, and maximum pressurization occurs just beneath the PCB. The temperature at the crater decreases monotonically throughout the subsequent decade; meanwhile, temperature and pore pressure below the PCB increase for several months and then remain high (Fig. 6).

Fig. 5

Distribution of change in a temperature and b pressure 1 year after conduit obstruction (Run 4, Mʹ = 1000 kg s⁻¹, kʹ_PCB = 1.0 × 10⁻¹³ m²). Arrow indicates the location where maximum pressure change occurs in 10-year period

Fig. 6

Temporal change in a temperature at the crater (red line) and the observation point (black line) and b pressure at the observation point after conduit obstruction (Run 4, Mʹ = 1000 kg s⁻¹, kʹ_PCB = 1.0 × 10⁻¹³ m²)

When an increase in hydrothermal fluid is combined with conduit obstruction, several possible patterns of change in the crater temperature can be observed (Fig. 7). A small decrease in permeability (to 10⁻¹¹ m²) accompanied by doubled hydrothermal fluid flux shows results similar to the change induced by increase in hydrothermal fluid flux only. The crater temperature increases for several years and then remains high. Figure 8 shows the distributions of change in temperature and pore pressure at 1 year for such a case (Run 11); the distributions of change are also similar to those influenced by only hydrothermal fluid flux increase (Run 10, see Fig. 3). In a second pattern, the crater temperature decreases monotonically through the decade and is similar to those resulting from conduit obstruction only. However, the distributions of temperature and pore pressures changes differ. Unlike Run 4 (Fig. 5), Run 14 (Fig. 9) shows an increase in temperature around the conduit below the PCB and a small increase in temperature in the edifice below the water table owing to the doubling of hydrothermal flux accompanying the large decrease in permeability (to 10⁻¹³ m²) (Run 14). The small increase in pore pressure throughout the volcano below the water table is induced by the increase in hydrothermal fluid flux. In a third pattern, the crater temperature shows complex temporal change, decreasing gradually for several years and then recovering to initial value. The broader distribution of changes is similar to those of the second pattern (Fig. 10, Run 12).

Fig. 7

Temporal change in temperature at the crater: a constant initial condition and Runs 1–4, b Runs 5–9, and c Runs 10–14

Fig. 8

Distribution of change in a temperature and b pressure 1 year after increasing of hydrothermal fluid flux with minor conduit obstruction (Run 11, Mʹ = 2000 kg s⁻¹, kʹ_PCB = 1.0 × 10⁻¹¹ m²). Arrow indicates the location where maximum pressure change occurs in 10-year period

Fig. 9

Distribution of change in a temperature and b pressure 1 year after increasing of hydrothermal fluid flux with strong conduit obstruction (Run 14, Mʹ = 2000 kg s⁻¹, kʹ_PCB = 1.0 × 10⁻¹³ m²). Arrow indicates the location where maximum pressure change occurs in 10-year period

Fig. 10

Distribution of change in a temperature and b pressure 1 year after increasing of hydrothermal fluid flux with moderate conduit obstruction (Run 12, Mʹ = 2000 kg s⁻¹, kʹ_PCB = 1.0 × 10⁻¹² m²). Arrow indicates the location where maximum pressure change occurs in 10-year period

Figure 11 summarizes our numerical simulations, focusing on temperature change at the crater and the maximum change in pressure at shallow depth (above the regional water table) over the simulated 10-year period. The horizontal axis represents the amplitude of the stepwise increase in hydrothermal fluid flux at the bottom of the system (Mʹ/M₀, where Mʹ and M₀ represent the elevated mass flow rate and the initial one, respectively). The vertical axis represents the degree of permeability reduction at PCB (k_h/kʹ_PCB, where k_h and kʹ_PCB represent the permeability of host rock and that of PCB after the conduit obstruction, respectively). When the permeability in the conduit decreases sufficiently, the crater temperature decreases despite the increase in hydrothermal flux from the deep part of the system. Moreover, Fig. 11 indicates that a large increase in hydrothermal flux accompanied by a large reduction in permeability induces significant pressurization in the shallow part of the edifice despite the decrease in the crater temperature. Maximum pressure change occurs in the conduit just beneath the PCB following conduit obstruction. When obstruction is not introduced (e.g., only increase in the hydrothermal flux is imposed), the maximum pressurization occurs just above the water table in the vicinity of the conduit.

Fig. 11

Summary of the numerical simulations. Symbols denote the pattern of the temperature change at the crater during the simulated 10 years (circles, monotonic increase; squares, monotonic decrease; triangles, fluctuated). The color scale represents the maximum amplitude of pressurization in the shallow part of the edifice over the simulated 10-year period. The numerals specify the run number (Table 2)

Both increases in hydrothermal flux from the deep part of the system and reduction in permeability in the shallow part of the conduit cause pressurization and heat accumulation in the shallow part of the edifice. Crater temperature shows complex behaviors depending on the amount of increase in hydrothermal flux and the amount of reduction in permeability.

Discussion

Clues for observation and prediction of phreatic eruption

It is difficult to predict phreatic eruptions by only the crater temperature alone. One should be wary of the potential for large eruptions when crater temperature decreases. Indeed, eruptions of Scale 4 (Scott 2013) occurred frequently after crater temperature decreased at Ruapehu between 1940 and 2012 (Strehlow et al. 2017). Observations of ground deformation, which can indicate pressure change in the shallow part of the edifice, and crater temperature can be combined to predict phreatic eruption and understand the activity of the hydrothermal system.

Coupling of ground deformation and crater temperature could make it possible to quantitatively investigate hydrothermal system activity. The location of maximum pressure change (see Figs. 3, 5) could reveal whether or not conduit permeability has decreased. We note that the apparent depth of the shallow pressure source at Tokachidake, Japan, varies with time (Takahashi et al. 2017). At Tokachidake, ground deformation suggesting pressure increases in the shallow part of the edifice has been continuing for nearly 10 years and, since 2006, declining crater temperature has been recognized (Takahashi et al. 2017). Decreasing permeability in the shallow part of the volcano may control volcanic activity (Tanaka et al. 2017). The varying depth of the pressure source may reflect repeated occurrences of permeability reduction in the shallow part of the conduit and/or an increase in hydrothermal flux from the deep part of the system. Numerical simulation coupled with a ground deformation model constrained by field observations would enable quantitative examination of this hypothesis.

Conclusion

Changes in hydrothermal fluid flux and permeability reduction can induce heating and pressurization beneath volcanic craters below the depth of conduit obstruction. The change in temperature at the crater itself is controlled by the amplitude of the increase in hydrothermal fluid flux from the deeper part of the system and the degree of permeability reduction.

In this study, simulations were carried out using a simple structure. However, the response of the edifice to changes in of the hydrothermal fluid flux will be influenced by the heterogeneous permeability structures (Todesco et al. 2010; Currenti et al. 2017). The distribution of permeability affects the timing and amplitude of changes in temperature and pressure through space and time (Todesco et al. 2010). Pressure- and temperature-dependent permeability will also influence the behavior of the system (Coulon et al. 2017), as will the porosity of the rock, the topography, the composition of the hydrothermal fluid, the presence or absence of a crater lake, and other factors.

We suggest that it is difficult and potentially misleading to predict eruptions from increases in crater temperature. One should be wary of the potential for large eruptions when crater temperature decreases. Coupling observations of crater temperature change with observations of ground deformation can clarify mechanisms and help predict phreatic eruption. We will couple hydrothermal simulation with a model of ground deformation in future research.