Permeability-control on volcanic hydrothermal system: case study for Mt. Tokachidake, Japan, based on numerical simulation and field observation

We investigate a volcanic hydrothermal system by using numerical simulation with three key observables as reference: the magnetic total field, vent temperature, and heat flux. We model the shallow hydrothermal system of Mt. Tokachidake, central Hokkaido, Japan, as a case study. At this volcano, continuous demagnetization has been observed since at least 2008, suggesting heat accumulation beneath the active crater area. The surficial thermal manifestation has been waning since 2000. We perform numerical simulations of heat and mass flow within a modeled edifice at various conditions and calculate associated magnetic total field changes due to the thermomagnetic effect. We focus on the system’s response for up to a decade after permeability is reduced at a certain depth in the modeled conduit. Our numerical simulations reveal that (1) conduit obstruction (i.e., permeability reduction in the conduit) tends to bring about a decrease in vent temperature and heat flux, as well as heat accumulation below the level of the obstruction, (2) the recorded changes cannot be consistently explained by changing heat supply from depth, and (3) caprock structure plays a key role in controlling the location of heating and pressurization. Although conduit obstruction may be caused by either physical or chemical processes in general, the latter seems more likely in the case of Mt. Tokachidake.Graphical abstract . .


Background
Hydrothermal systems are relevant to a variety of hydrovolcanic phenomena, including hydrovolcanic eruptions (Barberi et al. 1992;Germanovich and Lowell 1995;Browne and Lawless 2001). In this study, we focus on phreatic eruptions, defined by Barberi et al. (1992), as the eruption style caused by underground aquifers, whether or not they are phreatic or geothermal, without direct involvement of a magma body. Although most phreatic eruptions are local phenomena affecting a limited area, they can be accompanied by hazards such as ballistic blocks, base surges, and debris avalanches, and precursors are often subtle (Barberi et al. 1992). Germanovich and Lowell (1995) suggested a typical sequence of processes for phreatic eruptions: (1) dike injection, (2) hydrothermal system heating above the dike, and (3) rapid crack propagation due to thermal stress. This implies that the hydrothermal system may experience some change in state prior to phreatic eruptions. However, changes in the hydrothermal system do not always lead to explosive surface manifestations. Thus, improved understanding of changes in state in the hydrothermal system contributes to the assessment of whether phreatic eruption is likely.
Instrumental monitoring provides insight into ongoing changes in hydrothermal systems beneath specific volcanoes, whereas numerical simulations can provide useful clues as to the physical and/or chemical implications of such changes. Such clues can be obtained by sensitivity analysis (varying boundary conditions and parameters) and by forward modeling of the system's temporal evolution. Hydrothermal simulators have most often been used to assess the productivity of geothermal fields, but have also been applied to volcanic systems. Several previous studies compare volcano-monitoring records to numerical simulations to develop process understanding. For example, Christenson et al. (2010) discussed the eruption of Mt. Ruapehu, New Zealand, in 2007, by analyzing the effect of low-permeability seal formation on the CO 2 flux and temperature of the crater lake. Fournier and Chardot (2012) concluded that ground deformation at White Island, New Zealand, during 2002-2006-2009 was caused by a pore pressure increase at depth, rather than thermal expansion. However, few if any studies so far have used magnetic field monitoring as reference data for numerical simulations.
Many previous studies have reported magnetic field changes accompanying volcanic activity (e.g., Dzurisin et al. 1990;Sasai et al. 2002;Hurst et al. 2004). Such magnetic field changes may reflect temperature changes (thermomagnetic effect), stress changes (piezomagnetic effect), fluid-driven telluric current (the electrokinetic effect), or chemical processes such as hydrothermal alteration of magnetic minerals (Sigurdsson 2000). A large thermomagnetic effect implies a substantial difference between heat supply from depth and surface and/or nearsurface heat discharge.
On the basis of observations of the magnetic total field and surface heat discharge, we focus on the heat balance at shallow depths. At Mt. Tokachidake Volcano in central Hokkaido, Japan (Fig. 1), continuous demagnetization beneath the active crater area has been observed . We aim to investigate the mid-to longterm (years to decade) behavior of the hydrothermal system at Mt. Tokachidake based upon the records of the magnetic total field and thermal activity, using numerical simulations to better interpret the recent volcanic activity summarized in the following section.

A brief summary of historical eruptions and recent activity of Mt. Tokachidake
The major historical eruptions of Mt. Tokachidake in the modern era took place in 1926, 1962, and 1988-1989(Katsui et al. 1990  which brought about the collapse of the NW sector of the preexisting central cone followed by hydrothermal surges and snow melting, resulting in lahars (Uesawa 2010). This eruption left the collapsed topography now called Taisho crater. Intermittent small eruptions continued until 1928. In the 1962 eruption, hydrothermal eruptions also preceded the magmatic eruption of the main phase. The column of the sub-Plinian eruption reached 12 km above the crater. New craters named 62-0, 62-1, 62-2, and 62-3 were formed in the course of this eruption. The total erupted volume was estimated as 7.1 × 10 7 m 3 (Katsui et al. 1963). The 1988-1989 eruption started with a hydrothermal eruption from crater 62-2 in December 1988. Phreatomagmatic eruptions with small-scale pyroclastic flows occurred intermittently until March 1989. Most recently a very small-scale eruption with a colored plume was witnessed in 2004 (JMA 2013).
Crater 62-2 continues to host fumarolic activity but, according to reports that are regularly released to the public by the Japan Meteorological Agency (JMA), both plume height and temperature have gradually decreased since 2000. Meanwhile, localized ground inflation around crater 62-2 has been recognized since 2006 (observations performed by Geological Survey of Hokkaido, Hokkaido University and JMA). In addition, annual magnetic total field measurement campaigns by a group including the authors began in 2008 and revealed continuous demagnetization, suggesting heating at a shallow depth beneath the crater. Based on results from the first 2 years, Hashimoto et al. (2010) proposed that the localized demagnetization and deformation occur where hot gas (primarily water vapor) from depth condenses, releasing its latent heat, and then flows downward as liquid. Under this hypothesis, the amount of heat required to bring about the observed demagnetization would match the gradual decrease in heat discharge from crater 62-2 observed since 2000. Demagnetization would not necessarily be a consequence of enhanced heat supply from depth. However, their model is a conceptual one and does not account for the localized release of latent heat at a certain depth, nor does it provide mechanism for reducing heat discharge from the crater. A fully viable explanation might be common to other volcanoes where remarkable magnetic field changes during quiescent periods have been reported (e.g., Kanda et al. 2010;Kuchi-no-Erabujima Volcano).

Magnetic total field
Intermittent repeated magnetic total field measurements performed once or twice a year in collaboration with the JMA and the Geological Survey of Hokkaido reveal continuous demagnetization, suggesting heat accumulation beneath crater 62-2 since 2008 (Fig. 2). The decrease in the magnetic total field in the southern part of the study area and the increase in the northern part are typical of field changes due to demagnetization. Since the spatial extent is limited to crater 62-2 and its vicinity, the causal demagnetization source was likely shallow. Assuming a spherical demagnetization source, Hashimoto et al. (2010) estimated a source depth of 150 m beneath crater 62-2 (star symbol in Fig. 2). The white bars and contour lines in Fig. 2 indicate changes in the total magnetic field that were calculated by assuming a synthetic demagnetized sphere. Hashimoto et al. (2010) proposed that the changes were likely due to the thermomagnetic effect; localized inflation had been observed during the same period (Volcanic Observations and Information Center, Sapporo District Meteorological Observatory, JMA 2016), but stress-induced piezomagnetism would have different polarity.
On the basis of subsequent observations, the electrokinetic effect is also unlikely, as the cumulative amplitude of the magnetic field changes is quite large and still developing; unrealistic continuous enhancement of shallow hydrothermal flow would be necessary to cause such large effect. On the other hand, near-surface geochemical destruction of magnetic minerals cannot be completely excluded, and may contribute to local changes if TKNM TKSM 62-2 Crater Taisho Crater Fig. 2 Comparison of calculated (white bars) and observed (black bars) total magnetic field changes, after Hashimoto et al. (2010). The star symbol represents the estimated source of demagnetization. The direction of magnetization is assumed to be parallel to the present geomagnetic field at Mt. Tokachidake (inclination 57.4°, declination −9.8°). Contour lines are modeled at 1750 m above sea level, and TKNM and TKSM are the continuous observation points installed in 2014 such destruction takes place in the vicinity of measurement sites. However, we mapped local anomalies by double-checking supplemental measurements in the neighborhood of each point and infer that most chemical demagnetization is very localized. We thus agree that the thermomagnetic effect is likely the dominant cause for the overall field changes around the crater. Figure 3 shows the temporal changes plotted as a simple difference of the magnetic total fields between two continuous monitoring sites (TKSM and TKNM) deployed in 2014. Assuming the thermomagnetic effect is dominant, the decreasing trends represent a temperature increase, or expansion of a high-temperature zone, beneath crater 62-2. Thus, the plot indicates continuous heat accumulation since at least 2008. The temporary reversal in mid-2015 reflects contamination due to expansion of steaming ground that approached the southern station (TKSM). The reversal ceased in late 2015 when the area of streaming ground stabilized.

Thermal activity
Since the moderate eruptive events in 1988-1989, surficial thermal activity at crater 62-2 has peaked twice without eruption, in 1995 and 2000 (Fig. 4). Plume height was observed by JMA at Camera 1 ( Fig. 1) until 2014 and at Camera 2 since 2014 and has gradually decreased since 2000. Plume height was 200-300 m in the early 2000s and had decreased to 50-100 m by the early 2010s. Although plume height does not correspond directly to heat discharge, because it is strongly influenced by instantaneous wind speed and direction, taking daily maxima and monthly averaging greatly reduces variance due to wind effects, so that the long-term trend may be used as a proxy of fumarolic heat discharge rate. JMA and the Geological Survey of Hokkaido have observed fumarolic temperature at crater 62-2 using infrared radiation thermocameras. In the early 2000s, temperatures close to 500 °C were recorded, and although the temperature Year 40 nT phreatic eruption in 2004), the overall trend is decreasing and the temperature is now near 100 °C. Both plume height and temperature reflect gradual reduction of heat discharge from the crater.

Numerical simulations
Mineralogical sealing seems to play a key role in phreatic eruption at some volcanoes (Raoul Island, New Zealand: ; Ruapehu, New Zealand: Christenson et al. 2010). We explore this conceptual but broadly cited hypothesis in our numerical model of the hydrothermal system at Mt. Tokachidake, where we invoke conduit obstruction due to permeability reduction to explain the decrease in heat discharge from the crater and consequent heating beneath it. The goal of the numerical simulations are to investigate (1) whether conduit obstruction can cause heat accumulation beneath the active crater, gradual cooling of the crater fumaroles, and waning of the plume height; (2) the effects of a specific structure such as a caprock or sealing layers; and (3) whether such conduit obstruction is likely to cause other detectable changes (e.g., ground deformation). We use STAR (Pritchett 1995) as the numerical code and HOTH2O (Pritchett 1994) as the equation of state over a temperature range of 0-800 °C. These tools enable us to calculate the heat transfer and mass flow rates of H 2 O (gas, liquid, and two phases) in the hydrothermal system. Pressure, temperature, and phase state at the center of each model cell are obtained by solving conservation equations for mass and energy. We performed mainly 2-D modeling in which the flow is confined to a cross section along a line through the summit of Mt. Tokachidake and crater 62-2 (A-B line in Fig. 1).
The modeled region is gridded in a rectangular coordinate system (Fig. 5). Although heat and mass transfer orthogonal to the cross section cannot be taken into account in the 2-D model, the main slope of the topography is consistent with the line-of-section. In the following calculations, the thickness in the orthogonal direction is assumed to be 250 m. Temperature and pressure are maintained constant at the ground surface (15 °C and 1.013 × 10 5 or 1.013 × 10 3 Pa; the latter pressure corresponds to vapor-phase incoming fluid) and on the vertical boundary at the downstream side (i.e., at x = 0 m; hydrostatic conditions and 30 °C km −1 ). Thermally insulating and hydraulically impermeable conditions are imposed at the bottom and upstream vertical boundaries (i.e., x = 5000 m). Meteoric recharge is injected at the land surface at a constant rate (2.8 × 10 −5 kg m −2 s −1 at 15 °C) in accordance with the local precipitation, and a constant heat flow of 6.9 × 10 −2 W m −2 is supplied at the base of the model.
We calculate the time-dependent physical states (temperature, pore pressure, mass flow rate of liquid/vapor, and heat flow) over four steps. First, assuming homogenous permeability for the entire space and boundary conditions of 15 °C and 1.013 × 10 5 Pa at the ground surface, the model is run for 10 5 years until it reaches a quasi-steady state. Second, we reduce pressure at the ground surface from 1.013 × 10 5 to 1.013 × 10 3 Pa and run the model until another quasi-steady state is achieved (10 5 years). This step is used to induce an unsaturated zone within a shallow part of the edifice, since HOTH2O equation of state cannot explicitly deal with air. The unsaturated zone is treated as a two-phase zone of H 2 O over a less-than-atmospheric pressure range. . In Case C, Region 2 has a higher permeability than the host rock (Region 1). A permeability reduction at PCB (permeability changeable block) is imposed at the time of conduit obstruction Third, the high-permeability conduit and the surrounding low-permeability caprock are introduced in the middle of the slope and at the same time hydrothermal fluid (350 °C, 75 kg s −1 ) is injected at the bottom of conduit to reproduce fumarolic heat discharge of about 100 MW, a value estimated by applying the plume-rise method (Kagiyama 1978) to time-lapse photographs of the plume in 2004 (picture data courtesy of JMA), after reaching a quasi-steady condition (10 3 years). Fourth and finally, an abrupt reduction of permeability is imposed at a particular depth in the conduit, and the system response is observed for about 10 years. We explore four cases (A to D) with different combinations of permeability. The high-permeability (5.0 × 10 −10 m 2 ) column described above is common to all four cases, and a particular block in the conduit (PCB) is treated as an instantaneous low-permeability plug, hindering upward migration of hydrothermal fluids. In Cases A and B, we apply a relatively high permeability (5.0 × 10 −14 m 2 ) to the bulk of the host rock (Region 1), so that an unsaturated zone occurs above 900 m elevation (Fig. 6, for a uniform and constant porosity of 0.1 and permeability of 5 × 10 −14 m 2 throughout the model space). We lack concrete field evidence for the thickness of the unsaturated zone other than the existence of natural hot springs at around 600-1200 m asl. Additional geophysical observations and research would be needed to reveal hydrological structure within Mt. Tokachidake and to run more realistic simulations in the future. An unsaturated domain does not occur for the lower-permeability (5.0 × 10 −16 m 2 ) host rock in Cases C and D. We introduce a sealing structure or caprock surrounding the conduit in Case A only. For the shape and location of the caprock, we refer to the magneto-telluric modeling of Yamaya et al. (2010), who invoked hydrothermal fluid to explain the very conductive (~0.5 Ωm) column (150-700 m depth) beneath crater 62-2, and also imaged an overlying less conductive zone (10 Ωm), implying sealing structure. For comparison, a sealing structure is not implemented in Case B. A relatively high-permeability zone (Region 2) is introduced between the conduit and the caprock to represent altered and fractured rocks in Case C only. Table 1 summarizes the horizontal permeability of each region for the four cases. For Regions 1 and 2, horizontal permeability is ten times greater than vertical, and in the conduit and sealing zone permeability is isotropic. We assume a typical range of permeability reported in the literature for other volcanoes (10 −18   to 10 −10 m 2 , Saar and Manga 1999;Sruoga et al. 2004), because specific data are not available for Mt. Tokachidake. Although a wide range of variation in porosity is also likely within the actual edifice, we assume for simplicity a uniform porosity of 0.1. Other rock properties are assumed constant and listed in Table 2. We use the subsurface temperature distribution after conduit obstruction to calculate the magnetic total field at the ground surface, based on the temperature dependence of magnetization of rock samples from Mt. Tokachidake (original data from Soga 1997). We derive the following approximate relationship that converts temperature T [°C] (0-350 °C) to magnetization J [A m −1 ] for each model block: where J 0 is the magnetization at 0 °C, assumed to be 10 A m −1 here. The magnetic total field anomaly at any observation point (x, y, z) is calculated by integrating the contribution from each block using the prism model of Bhattacharyya (1964). Figures 7 and 8 show the distribution of temperature, gasand liquid-phase mass flow rates, and liquid saturation (Fig. 8) under quasi-steady-state conditions just before we impose conduit obstruction for each model case (A to D). In all cases, fluid ascending along the highly permeable conduit vents as a gas phase. Cases A and B, in which relatively permeable host rock is assumed, show a wider high-temperature region (100-300 °C), hot spring discharge at the piedmont, and a gas-dominated zone beneath the crater. Neither extensive convection nor hot spring discharge is observed in Cases C and D, in which we assume lower host rock permeability. Comparison between Cases A and B shows that the caprock zone limits the lateral extent of high temperatures. Figure 9 shows conditions 1 year after conduit obstruction (PCB reduced from 5 × 10 −10 to 5 × 10 −14 m 2 ). Cases A and C result in heat-pipe conditions: the gas phase ascends in the conduit and the liquid phase flows down around the conduit at a comparable mass flow rate. Figure 10 shows differences in temperature between the quasi-steady state prior to conduit obstruction and

Results of numerical simulations
1 year after obstruction. All cases result in heating below the obstruction level and cooling above the obstruction level. The difference between Cases A and B owes to the presence (A) or absence (B) of a low-permeability caprock.
We use these temperature changes to calculate magnetic total field changes at 1750 m above sea level (Fig. 11). All four cases show a decreasing magnetic field in the southern domain and an increasing field in the northern domain, consistent with our field observations in terms of polarity and spatial scale, although the calculated change exceeds actual records by up to a factor of two. This inconsistency in amplitude may owe to the simplistic approximation of the thermomagnetic property of rock (Eq. 1). Figure 12 shows the simulated temporal changes for a decade following conduit obstruction. We adjust the initial heat discharge rate through the block corresponding to the vent (25 × 250 m; the smallest cell size in the model) roughly to 100 MW. Although this adjustment results in some disagreement between estimated vent temperature in the actual records (Fig. 4) and the simulations (Fig. 12), we rely more on heat discharge than temperature, because the dominant H 2 O component of the plume carries most of the thermal energy through its large latent heat. We consider temporal changes in vent temperature to be more significant than its absolute values. In Fig. 12, Case A displays continuous demagnetization attenuating with time, vent temperature decreasing for the first 2 years and then stabilizing at a lower level, and heat flux decreasing over the first few months before recovering somewhat. In Case B, the behavior of vent temperature and heat discharge is similar in Case A, but whereas Case A displays gradual demagnetization that does not recover, Case B shows less overall demagnetization, with fluctuations over a time frame of several years. Escape of hot fluid from Region 2 is likely responsible for the more modest magnetic changes in Case B, which lacks a low-permeability caprock. Cases C and D exhibit early time reductions in heat flux, similar to Cases A and B, but heat flux returns to the initial level after 2 years. Early time decreases in vent temperature for Cases C and D exceed those in Cases A and B, but vent temperature rapidly recovers to the initial level, with a subsequent large-amplitude fluctuations in Case C. Such later-time recoveries in heat flux and vent temperature are inconsistent with field observations. The continuous demagnetization in Cases C and D is generally in agreement with field records, although demagnetization is faster in Case C due to the higher permeability in Region 2. Case A (relatively permeable host rock with a caprock around the conduit) is the best fit to the recent situation at Mt. Tokachidake.

Table 2 Rock properties used in all cases (A-D)
Property Value Another set of simulations varied the degree of conduit obstruction (Fig. 13). Larger reductions in permeability lead to larger demagnetization (Fig. 13, upper  panels). However, no further demagnetization is seen for permeability reduction beyond four orders of magnitude. Vent temperature and heat flux show no significant changes, or quick recovery in the first 2 years after obstruction, when the permeability decrease is less than three orders of magnitude. For even lower permeability, no such recovery occurs. The difference in the response can be interpreted as follows: conduit obstruction causes a temporary decrease in the gas discharge through the conduit, so that both temperature and pore pressure in the underlying conduit will increase. When the permeability reduction is sufficient small, the elevated pore pressure can sustain flow through the conduit. Matching the observations at Mt. Tokachidake since 2000 (continuous demagnetization, gradual depression of vent temperature and heat flux) requires a considerable degree of permeability reduction.

Fluctuation of heat supply from depth
In the previous section, we demonstrated that conduit obstruction could cause heat accumulation beneath the crater and reduction of heat discharge from the vent. In general, changes in the supply rate of hydrothermal or magmatic fluid are often suspected when any kind of elevated volcanic activity is observed. We have assumed a constant fluid supply from depth and have not yet investigated how a change in the supply rate would affect system behavior. Therefore, we perform additional simulations varying the fluid input from depth by plus or minus one order of magnitude at t = 0 (Fig. 14). A tenfold increase in the mass supply rate causes demagnetization with increases in vent temperature and heat flux. Abrupt tenfold reduction of the mass supply rate causes remagnetization accompanied decreases in vent temperature and heat flux. Neither of these cases is consistent with observations at Mt. Tokachidake.

Ground inflation at the crater
Continuous ground inflation at the crater has also been observed since 2006 (Volcanic Observations and Calculated temporal changes over 10 years following conduit obstruction: differential magnetic field between TKSM and TKNM, vent temperature, and heat flux Information Center, Sapporo District Meteorological Observatory, JMA 2016). Pore pressure distribution is an output of our simulations that is useful in interpreting the ground deformation (Fig. 15). For all the cases (A-D), very local and shallow pore pressure reduction of about 0.1 MPa occurs above the obstruction, with pressurization below the obstruction. The vertically elongated pressurized domain (up to 1 MPa pressure increase) is sharply bounded by the low-permeability caprock, especially in Case A. In Case B, the change is smaller and broad, because Case B lacks a low-permeability seal around the conduit. Case A shows some depressurization, on the order of 10 −2 MPa, on the downstream side, where hot spring discharge occurs, whereas Case B exhibits almost uniform pore pressure increase below the water table. While the pressurization in Cases A and B is ~0.1 MPa, pressurization reaches several MPa in most parts of the model space in the lower-permeability Cases C and D. This might be an overestimate, as the pore pressure changes in these calculations are comparable to the lithostatic load. Ground deformation due to shallow sources is in general strongly affected by topography. Although a STAR post-processor is available to calculate ground deformation due to subsurface pore pressure and temperature changes (e.g., Ishido et al. 2015), further development is required to treat detailed topography. Calculated temporal changes over 10 years following conduit obstruction for varying degrees of obstruction. Simulated conditions are otherwise the same as Case A: differential magnetic field between TKSM and TKNM, vent temperature, and heat flux. Legend indicates the permeability of the obstruction after permeability reduction. For instance, 5 × 10 −11 m 2 indicates that permeability decreases from the initial value of 5 × 10 −10 to 5 × 10 −11 m 2

Time [year]
In addition, given the simple assumption of a uniform porosity in this study and the difficulty of assigning the distribution of pertinent poroelastic parameters, we do not think it very meaningful to compare such simulated and observed deformation. However, it is clear that pore pressure increase at a shallow depth in and around the conduit could cause localized surface deformation.

Implications for the recent activity at Mt. Tokachidake
According to our numerical simulations, conduit obstruction decreases vent temperature and heat flux and results in heat accumulation below the level of obstruction. Conduit obstruction may be caused by either physical or chemical processes. Candidate physical processes include mechanical reworking of surface deposits (e.g.,

Decrease in supply
Increase in supply  Candidate chemical processes include deposition of native sulfur and/or hydrothermal minerals in the conduit. Hurst et al. (1991) studied cyclic behavior of the crater lake at Ruapehu (New Zealand) and suggested that the temperature-dependent viscosity of native sulfur played a key role in the cyclic activity. This mechanism may also be effective at Mt. Tokachidake. In June 1923, a molten sulfur pond (locally called "Yunuma") appeared south of the so-called central cone, and there was a sulfur mine at Mt. Tokachidake until the 1962 eruption. The discharge of sulfur dioxide from crater 62-2 was 210 t day −1 on July 7, 2007 (Mori et al. 2006) and 140 t day −1 in 2006 (Mori et al. 2013). At geothermal fields, buildup of solid silica often reduces the permeability of pipelines and wellbores. The solubility of amorphous silica decreases remarkably when it is cooled (Fournier and Rowe 1977), so cooling in the conduit may cause silica scaling that reduces permeability. This permeability reduction would cause further cooling above the obstruction level, generating positive feedback. However, silica scaling is inhibited in low pH environments (Gallup and Barcelon 2005), and the mid-slope hot spring at Mt. Tokachidake has pH between two and four (Takahashi et al. 2015). Many other types of hydrothermal alteration, such as clay formation, cause volume increases that tend to decrease permeability. It might be possible to clarify the dominant chemical processes under the specific conditions of Mt. Tokachidake by using other kinds of numerical simulation that incorporate fluid-rock interaction. Christenson et al. (2010) performed such simulations for a simple onedimensional case. Based on the records of CO 2 flux and temperature of the crater lake at Ruapehu, and numerical simulations, they discussed the 2007 eruption at Ruapehu. 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. In their calculations, deposition of sulfur reduced permeability from 10 −12 to 10 −17 m 2 . Although the same mechanisms are not likely to translate completely to Mt. Tokachidake (there is no stable lake in crater 62-2), meteoric water inflow to the conduit resulting from decreased hydrothermal upflow from depth could cause deposition of sulfur.

Validity of the two-dimensional approach
To this point in our study, we have dealt with 2-D models which exclude mass and heat flux in the Y-direction, possibly causing some overestimation of temperature and pore pressure. We performed additional 3-D calculations for a couple of representative cases (Fig. 16). Because 3-D simulation is quite time-consuming (10 times slower than the 2-D model), it is not practical for sensitivity analysis. Thermally insulating and hydraulically impermeable conditions are imposed at y = 0 and y = 4775 m. The hydrothermal input (350 °C, 170 kg s −1 ) at the bottom of the conduit (x = 3900-3925 m, y = 2375-2400 m, z = 0 m) is adjusted to reproduce a typical fumarolic heat discharge of about 100 MW before conduit obstruction. The other boundary conditions and numerical steps are the same as those in the previous 2-D modeling example. The combinations of permeability are those used for Case A (with caprock) and Case B (without caprock). The distribution of temperature at quasi-steady state prior to conduit obstruction is similar to that in the 2-D model in both Case A and Case B. In Case A, the lowpermeability caprock hinders lateral diffusion of hot fluid. Figures 17 and 18 show the differences in temperature between the quasi-steady state prior to conduit obstruction and that 1 year later on the Y-Z and X-Z planes, respectively. Some differences are recognized between the 2-D model and the 3-D model in Case B. In the 2-D model, Region 2 and its surroundings area diffusively heated (Fig. 10b). In the 3-D model, heating is concentrated within the conduit and at the level just beneath the permeability obstruction. Figure 19 shows the magnetic total field change calculated from the simulated temperature change. Although the change in temperature distribution is different between the 2-D and 3-D models, the magnetic total field changes show a similar bipolar pattern, again consistent with our field observations. Figure 20 shows modeled temporal changes for the decade following conduit obstruction. Essential features are similar to the 2-D model (Fig. 12). These limited comparisons suggest that the 2-D approach is reasonable for investigating the shallow system around the conduit at Mt. Tokachidake. Figures 21 and 22 show the 3-D distribution of pore pressure changes 1 year after conduit obstruction. In the deeper levels below the water table, pressurization is about one order of magnitude less than in the 2-D models. This is entirely due to lateral flow of fluid mass and heat in the Y-direction. We should be aware of this effect when we evaluate the ground deformation, since the differences are conspicuous. Detailed modeling of ground deformation may provide critical clues to understanding the shallow structure beneath the crater and the system's response to unrest.

Conclusions
We modeled a volcanic hydrothermal system by means of numerical simulation with three key observables as reference: the magnetic total field, vent temperature, and heat flux. Volcanic activity at Mt. Tokachidake since the early 2000s was considered as a case study. Field observations there have revealed continuous demagnetization, suggesting heat accumulation beneath the active crater, Calculated temporal changes for the 3-D model over 10 years following conduit obstruction: differential magnetic field between TKSM and TKNM, vent temperature, and heat flux gradual cooling of the crater fumaroles, and waning of the plume height from the crater (a proxy for heat loss). Our numerical simulations of the hydrothermal system reveal that conduit obstruction (i.e., permeability reduction in a conduit) could cause these changes. We also confirmed that the changes could not be consistently explained by fluctuation of heat supply from depth. A caprock structure that confines hot gases and hydrothermal water at shallow depth may play a key role in controlling the location of heating and pressurization, consistent with the concept of fluid confinement prior to phreatic (or hydrothermal) eruptions. Without such a caprock structure, hot fluids are prone to percolate diffusely, so that localized and long-lasting pressurization and demagnetization cannot be achieved. Conduit obstruction may be caused either by physical or by chemical processes, although the latter seems more likely in the case of Mt. Tokachidake. Further investigation is necessary to understand the details of the chemical processes, perhaps by using numerical simulations in which fluid-rock interaction is incorporated. We expect that a numerical study coupling physics and chemistry would contribute to better understanding of the preparation processes of phreatic/hydrothermal eruptions.