Quantitative relationship between plume emission and multiple deflations after the 2014 phreatic eruption at Ontake volcano, Japan

The phreatic eruption of Mount Ontake in 2014 caused local-scale subsidence and a mass discharge of water–vapor plumes from vents. A previous study of InSAR data analysis modeled the local subsidence as a deflation of a shallow hydrothermal reservoir (~ 500 m beneath the vents), and speculated that it was associated with plume emission continuing just after the eruption. In addition, combination of the InSAR and GNSS data implies that another, deeper deflation source (~ 3–6 km beneath the vents) contributes to the baseline contraction of the GNSS data. In this study, we estimated daily mass flux of the emitting plumes using photographed images, and compared the temporal behavior of the discharged mass with that of deflation of the two sources in order to clarify their association. The temporal profiles of the shallow deflation volume and the discharge mass both show evidence of decay, but with different characteristics; the deflation volume progress was approximated by a single exponential decay with a long relaxation time (379–641 days), whereas the discharge mass displayed a sum of a linear trend and an exponential decay with shorter relaxation time (47 days). This discrepancy, along with GNSS data, suggests the contribution of a deep deflation source with a short relaxation time (20–40 days). Estimation of mass balance between the emitting plume and fluids discharged from both shallow and deep sources revealed that more than 70% of the discharged mass came from the deep source. Based on the estimated mass balance, phase state of the shallow reservoir was estimated as a single-phase, liquid-rich reservoir. The fast decay of the deep deflation may reflect rapid depressurization due to violent fluid discharge at the onset of the eruption. In contrast, the slow decay of the shallow deflation suggests that it had a minor role in the eruption. However, such a wet reservoir has the potential to induce volcanic hazard such as snow-melting lahar for future eruptions, requiring monitoring the volcano, which will probably shift to pre-eruptive re-pressurized phase, until the future eruption.


Introduction
Phreatic eruptions occur frequently on volcanoes around the world. Despite their frequency, the triggering mechanisms and subsurface conditions leading to phreatic eruptions have not been well understood (Barberi et al. 1992;Oikawa et al. 2018). Germanovich and Lowell (1995) proposed a conceptual model suggesting that phreatic eruptions result from pressurization in volcanic, hydrothermal systems due to the injection of hot fluids from the magmatic region. Continuous geodetic observation may allow for the detection of pre-and syn-eruptive pressure changes within hydrothermal systems. However, it is often difficult to detect precursory geodetic signals because the source region of such deformations is usually shallow, and consequently the signal is spatially localized in a small area requiring monitoring of a very specific activity region. A frontier study of the detection of precursory ground inflation was achieved during the phreatic eruption of Mount Hakone in 2015 (Kobayashi et al. Narita et al. Earth, Planets and Space (2019) 71:145 2018), using space geodetic techniques such as Interferometric Synthetic Aperture Radar (InSAR). The high frequent observations started before the eruption, playing an important role in detecting ground deformations. However, even with such observations, it is still difficult to detect the precursors sufficient for understanding the mechanisms of phreatic eruptions.
Several studies have investigated ground deformations, usually subsidence, after phreatic eruptions, and most of these eruptions did not evolve into magmatic eruptions within the time frame of several years (Jousset et al. 2000;Lu et al. 2002;Nakaboh et al. 2003;Hamling et al. 2016;Doke et al. 2018;Miller et al. 2018;Narita and Murakami 2018). These subsidence are often accompanied by a continuous emission of water-vapor plumes, potentially linked with depressurization due to fluid discharge from shallow aquifers or hydrothermal reservoirs (Lu et al. 2002;Nakaboh et al. 2003;Miller et al. 2018). Although deformations are dominated by shallow, depressurized reservoirs, such fluids may not originate solely from shallow reservoirs. Several studies suggested that fluids were supplied from deeper parts. For example, a deflated, shallow reservoir was found after the 1995 phreatic eruption at Mount Kuju, south Japan, which was considered to be strongly related to the plume emission (Nakaboh et al. 2003). While the heat assessment revealed that the discharged heat from the shallow reservoir was insufficient to account for the total heat flux observed at the surface, it was inferred that the heat flux carried by fluid advection from a deeper, possibly magmatic, region was significant (Hashimoto et al. 2002). After the 2012 phreatic eruption at Mount Tongariro, New Zealand, InSAR and GNSS data revealed the existence of deflation sources at depths of 0.5 km and 2.5 km, respectively (Hamling et al. 2016). Furthermore, a strong and continuous gas emission at over 400 °C (Global Volcanism Program 2014) suggested a link between the discharge of both the shallow and the deep sources.
After the 2014 phreatic eruption at Mount Ontake, similar post-eruptive deflations at both shallow (~ 0.5 km) and deep (~ 3 km) depths were observed. The deformation was accompanied by active emitting vapor at the surface. Narita and Murakami (2018) speculated that the shallow deflation derived by InSAR might be caused by the mass discharge, but a quantitative relationship between the deflation sources and the emitting plume has remained unclear.
Mount Ontake, the second highest volcano (3067 m absl) in Japan, is one of the most active stratovolcanoes in the country. Hydrothermal activity has been recorded for over 300 years around Jigokudani Valley, where numerous post-and current-eruptive craters and fumaroles exist (Oikawa 2008). Historical eruptions occurred in 1979Historical eruptions occurred in , 1991Historical eruptions occurred in , 2007Historical eruptions occurred in , and 2014 in Jigokudani Valley, all of which were phreatic eruptions. Despite its small magnitude (VEI = 2), the 2014 eruption resulted in 58 deaths and five missing people. The most intense part of the eruption continued for 6 h on 27 September 2014, creating a series of linearly aligned crater chains. ALOS-2 InSAR observations revealed local-scale (~ 4 km) subsidence around new craters (Narita and Murakami 2018), with simultaneous plume emission occurring (Yamaoka et al. 2016) and continuing as of July 2017 (Fig. 1). The subsidence was attributed to a deflation source at a depth of 500 m below the vents, which might be caused by fluid discharge from a shallow, hydrothermal reservoir (Narita and Murakami 2018). However, it was observed that the shallow source could not account for all of the observed change in GNSS baseline length between Tanohara and Ochiai stations (Fig. 1b). Thus, it was inferred that a deeper, deflation source likely contributed to the baseline length change, and also possibly supplies water for the plume emission to the surface, similar to the eruptions at Mount Kuju and Mount Tongariro.
Previous studies at Mount Ontake suggested the presence of several pressure sources at depths greater than that of the InSAR-detected shallow source. There are several possible sources: a co-eruptive deflated source 6 km below the eruptive craters, which was derived from repeated leveling observations spanning the 2014 eruption (Murase et al. 2016); inflation sources in the 2007 eruption at 1 and 3 km depths below the craters (Takagi and Onizawa 2016); a VLP source in the 2007 eruption 2.4 km beneath the craters (Nakamichi et al. 2009), which were considered manifestations of activated hydrothermal systems. Each can be considered a candidate for supplying emitted vapors at the surface. A quantitative investigation on mass balance may provide further information about the ratio of contributions of these sources to the vapor emission.
The purpose of this study was to investigate the quantitative relationship between the shallow and deep deflation sources and the emitting plume after the 2014 Mount Ontake eruption. We focused on the detailed analysis of the relationship between mass discharge and deformation using new plume data that were made available to us, and compared it with previously discussed source models of shallow deflation (Narita and Murakami 2018). We also examined the GNSS data that may reflect a deep deflation.

Materials and methods
To analyze the relationship between the deflations and the plume emission, we compared the time series of the discharged water mass with that of deflation volume of the shallow source, which were derived from InSAR analysis results (Narita and Murakami 2018). We also compared them to available GNSS data, which show contraction over a relatively short baseline spanning the volcano suggesting a deflation at a deep part (Narita and Murakami 2018).
To see the basic characteristics of each time series, we fitted simple analytical functions to these data. We employed two types of functions: (1) an exponential function f 1 (t) = a 1 − e −t/τ , and (2) a combination of an exponential and a linear function f 2 (t) = a 1 − e − t τ + ct , where a , t, τ and c represent an amplitude of the exponential term, elapsed time (day) from 28 September 2014, relaxation time of the decay, and amplitude of the linear term, respectively.
For the deflation volume history of the shallow source, we used the source model proposed by Narita and Murakami (2018). As a best-fit source model, they estimated an almost spherical source with a deflation volume of 3.5 × 10 5 m 3 at a depth of 500 m between October 2014 and July 2015. We used dLOS time series of two orbits of the ALOS-2 satellite (Fig. 2a); path 126, right-looking ascending, and path 20, right-looking descending. Time series of the deflation volume were calculated by extrapolating the volume value between October 2014 and July 2015 assuming the same time series as that of dLOS of the two paths. We treated both the InSAR paths as independent representations of the time series of volume change of the shallow source. Note that the dLOS time series were plotted at the pixel showing significant dLOS (open circles in Fig. 2b, c).
In order to retrieve spatially long wavelength signals that deep pressure source may produce, we focused on the Tanohara-Ochiai (TO) GNSS baseline data, because it was the only baseline that showed a significant contraction after the 2014 eruption (Narita and Murakami 2018). The deflation of InSAR-derived shallow sources should also partly contribute to the baseline contraction. We calculated time series of the baseline contraction caused by the shallow deflation, using the contraction predicted by the finite element model of Narita and Murakami (2018). The contraction time series were same as those of the volume change in each path. We subtracted time series of the predicted baseline contraction, of which time functions were same as the dLOS time series, from the original TO baseline data, allowing us to recognize the remaining contraction that suggested a contribution from the deep part. We hereinafter refer this baseline data after the subtraction as the "corrected time series".
To obtain the discharge time series, we estimated mass flux of the discharged water contained in the plume by using plume rise methods (Briggs 1969;Kagiyama 1981), which allows for estimating the heat flux of the plume Q h , based on the following formulation: where u is the horizontal velocity of flowing plume and C 1 is coefficient depending on plume geometry. Both parameters can be estimated by tracking identical portion of the blowing plume over several photos. Note that we excluded the photos affected by upstream wind blowing along topographic slope, which could result in incorrect estimate of C 1 values. We estimated Q h from December 2014 to March 2017 on a daily basis using photo images taken from Takigoshi station (Fig. 3), which is operated by the Japan Meteorological Agency (JMA). The images were sampled at a rate of 10 s from 09:55 to 10:05 and from 14:55 to 15:05 on each day from 1st December 2014 to 31 March 2017. Although the images were taken every day during this time period, available images were limited because of (1) Q h = 2.6 × 10 4 u 3 C 3 1 , weather conditions and wind direction. As a result, only 5.4% of the total data (46 days/852 days) were available. Since heat flux of buoyant water-vapor plume mostly derives from latent heat discharged by vapor condensation that occurs immediately after emission from vents, the heat flux can be converted into mass flux by dividing by latent heat of vapor (Kagiyama 1981). We divided Q h by latent heat of vapor at 90 °C (2283 kJ/kg), which is boiling point at the vent elevation, and obtained discharge mass flux Q out . Since plume temperature after the eruption was relatively low (70-100 °C), as reported by Mori et al. (2016), the contribution of the sensible heat associated with vapor cooling was negligible. The photographic data sampled covered a period of 2 months after the eruption to the end of March 2017. From 1 day after the eruption (28 September 2014) to 22 November 2014, we referred to the data estimated by Terada (2014) and by Hashimoto and Tanaka (unpublished data). By combining these data sets, we are able to illustrate the complete time series of vapor emission for 2.5 years after the eruption. We integrated the linearly interpolated Q out and obtained the time series of cumulative discharge mass. Figure 4 shows the time series of the discharge mass Q out (t) and of the deflation volume of the shallow source ΔV s (t). Figure 5 shows the corrected time series in the TO baseline length ΔL d (t), which may only reflect the deep deflation signal. Each time series showed monotonic decay, but temporal characteristics of the shallow deflation were clearly different from that of the others. ΔV s (t) was approximated by a single, exponential decay curve f 1 (t), with a smaller AIC value than that of f 2 (t) in both of the paths (Table 1). In contrast, Q out (t) and ΔL d (t) required a combination of a linear and an exponential term f 2 (t), with a smaller AIC value than that of f 1 (t) ( Table 1). These time series showed fast decay for the first ~ 50 days after the eruption and then approached a linear decay (Figs. 4 and 5). Relaxation times τ for ΔV S (t) in both paths were 379-641 days, whereas Q out (t) and ΔL d (t) had significantly smaller τ values of 47 and 20-41 days, respectively (Table 1). Furthermore, temporal characteristics of Q out (t) and ΔL d (t) were significantly similar to each other ( Fig. 6a, b). These facts probably suggest a close relationship between the plume emission and the deep deflation.

Origin of the discrepancy of the relaxation times
Here, we discuss possible reasons that can cause such discrepancy among the relaxation times of ΔV s (t), Q out (t) and ΔL d (t). We start from a general relationship between volume change of a deformation source and the corresponding mass change. Assuming the source deflation is caused by mass discharge from the source, deflation volume V is related to corresponding mass outflow m as (Rivalta and Segall 2008;Segall 2010): where ρ, β f , and β c represent fluid density in the source, fluid compressibility, and compressibility of the deformed reservoir, which depends on reservoir shape and host rock rigidity, respectively. We further assumed that all of the mass discharged from eruptive vents was supplied from the source without any lateral leakage through the conduit wall. β c is unlikely to change within several years, if we assume no temporal change in ρ and β f in Eq.
(2), a time function of m depends on only that of V . In this case, if most of the plume mass came from the deflated source (Q out = Δm), then temporal changes and relaxation times of V and Q out should have been the same. However, our results do not support this conclusion. The clear discrepancy of the relaxation times can be

Fig. 3
An example photograph of plumes viewed from Takigoshi station taken on 15 August 2015. The photos were provided from JMA. By tracking of movements of the identical portion of the plume over consecutive images, we estimated wind velocity u and coefficient of the plume shape C 1 , in a local coordinate of x, horizontal distance from a reference point, and h, which is vertical distance from a reference point, and then plume heat flux Q h (see text for the detail of analysis procedure) Fig. 4 Time series of the discharge mass (Q out ) and the deflation volume (ΔV s ). Time plots correspond to the period after 28 September 2014. Plots of ΔV s are based on the estimation for the two independent orbits of ALOS-2 satellite (Narita and Murakami 2018). Dotted lines show best-fit functions of f 1 (t) resulted from curve fitting, whereas dashed lines show best-fit functions of f 2 (t) (see text for the detail). Note that for ΔV s , f 1 (t) and f 2 (t) fitting results are almost identical, which means no involvement of linear term. Relaxation time τ in exponential terms is indicated at best-fit time functions of each time series data attributed to a difference in water sources, which raises the possibility that another water source, rather than the shallow source, contributes to the discharge. In contrast, when the temporal change of ρ and β f are non-negligible, the time function of m may also depend on that of ρ and β f as well as ΔV. We checked an effect of temporal change in ρ and β f due to depressurization on the resulting time function of m (Additional file 1: Figure S1). We evaluated the possible changes in τ for m assuming reasonable thermodynamic conditions at a shallow depth (see Additional file 1: Text S1 for details). We calculated the values of τ for possible Path 20 combinations of parameter values within plausible range, and found that the smallest value of τ for Δm was 115 days (Additional file 1: Figure S1), which is larger than the observed value of τ for Q out (47 days). The discrepancy of the temporal behavior still remains even when considering the temporal changes in ρ and β f . Regardless of the cases with or without considering the temporal changes in ρ and β f , the discrepancy strongly suggests that most of the emitting plume derived from another water source, possibly one that is much deeper, which was inferred from the GNSS data, rather than the shallow deflation source showing significant subsidence signal in the InSAR data.

Simple mass balance between the emitting plume, the shallow and the deep deflation source
It is possible that the corrected time series of the GNSS baseline, showing a high correlation with that of the plume emission, reflect mass discharge to the vents with ongoing depressurization. If we assume no temporal change in location and geometry of the pressure sources, then the length change in the TO baseline, ΔL d (t) should have been proportional to the volume change of the deep source, ΔV d (t). We again applied Eq.
(2) for the estimation of m of the deep source. However, it is almost impossible to assign plausible values for all the parameters in Eq.
(2), since both depth and geometry of the source are uncertain due to the insufficiency of the geodetic data. For simplicity, below, we assume that the time function of m for the deeper source, Q d (t), depends on only that of ΔV d . Then, Q d (t) becomes proportional to ΔL d (t).
Under the assumption that shallow and deep sources are depressurized due to mass discharge through the vents (Fig. 7), simple mass balance can be expressed as: where Q out (t) is discharged mass of the plume from vents, Q s (t) is discharged mass from the shallow reservoir and Q d (t) is discharged mass from the deep reservoir. Note that because of difficulties of evaluating meteoric water influx into the shallow reservoir, this model simply assumes little contribution of the influx into the shallow reservoir, such as snow-melt water, precipitation and preexisting groundwater, which may be sometimes crucial to the mass balance estimation at shallow part of volcanoes (e.g., Miller et al. 2018). Since all the time functions in Eq. (3) are known, we can estimate the quantitative relationships between those terms, that is, how much Q s (t) and Q d (t) contribute to Q out (t) . Note that Q s (t) is expressed as the time function of the shallow deflation volume �V s (t) and that Q d (t) is expressed as that of the baseline length change, �L d (t) . This required converting the unit of the two functions into mass, by multiplying conversion factors. Then, Q out (t) is rewritten as: where A is a conversion factor. The subscripts of 1 and 2 correspond to the shallow and the deep source, respectively (see Table 1 for the parameter values of �V s (t) and �L d (t) ). The first term corresponds to Q s (t) and the second term corresponds to Q d (t) . For simplicity, we Comparison between the discharge mass and the baseline contraction in Fig. 5a (a) and Fig. 5b (b), respectively. Gray-shaded areas indicate range of estimated error of the discharge mass assume that A 1 and A 2 , which include density and compressibility, are temporally constant. Regarding A 1 and A 2 as unknown model parameters, we inverted them by the least square fitting, and finally, obtained best-fit values of Q s (t) and Q d (t) by using observed values of Q out (t) , �V s (t) and �L d (t) . The fitting result shows that Q out (t) contains more that 70% Q d (t) during the entire period of time (Fig. 8a, b; Table 2). The mass ratio Q d (t)/Q out (t) a b Fig. 7 Schematic illustration of mass balance between the two deflation sources and the discharged plume. Q s , Q d and Q out correspond to extruded fluid from the shallow source, the deep sources and the discharged plume, respectively. ΔV s is deflated volume of the shallow source. ΔL d indicates baseline contraction due to the deep deflation. Deformation at the surface caused by deflations of deep and shallow sources brings displacements of GNSS stations and change in LOS distance toward the satellite Fig. 8 a Fitting results for discharge mass by combination of two components: (1) best-fit function for the corrected GNSS time series (as shown in Fig. 5) and (2) best-fit function for the deflation volume time series of the shallow source: ΔV s (t) (as shown in Fig. 4). Red and blue colors correspond to satellite path of 126 and 20, respectively. Solid, dashed and dotted lines indicate synthetic Q out (t) (= Q d (t) + Q s (t)), Q d (t) and Q s (t). Gray-shaded area indicates range of estimated errors of the discharge mass. b Temporal change in the mass ratio of Q d (t) to Q out (t). Red and blue colors correspond to satellite path of 126 and 20, respectively. Value of more than 1.0 is due to fitting error is small just after the eruption, then gradually increases later on path 126, which is due to fitting error, and is almost constant (> 90%) on path 20 (Fig. 8b). This strongly suggests that the deep deflation source is likely to be the main source for the discharge mass of the plume.

Mechanism of the deflations
Deflation at a shallow part of a volcano is generally attributed to either pore pressure decrease due to fluid migration (Bonafede 1990;McTigue 1986;Miller et al. 2018) or to thermo-elastic response with cooling (Bonafede 1990;Furuya 2005; Wang and Aoki 2019). In the previous section, we have assumed only the depressurization as a plausible mechanism of deflations. Here, we discuss which mechanism is really dominant for the shallow and deep deflation source. The high correlation of the deep deflation with the plume emission ( Fig. 6) strongly supports that deep deflation can be induced by mass discharge to the surface through the eruptive vents. Considering gas emission of low temperature close to boiling point and low involvement of SO 2 gas, except immediately after the eruption ) led us to infer that this source may not be in a hot region like magma chambers, but in a part of the hydrothermal system. The thermo-elastic response for the shallow source can be insignificant. Temperature change ( T ) required for the thermo-elastic volume change is expressed as: where V t is the volume change with thermo-elastic deformation, α t is volumetric thermal expansion coefficient, V t is rock volume associated with the temperature change. Substituting a typical value of α t 1 × 10 −5 • C −1 , the observed value of V 7 × 10 5 m 3 and the constrained value of V t ≤ 10 8 m 3 , which is maximum volume when we assume spherical region of which center is at a depth of 500 m beneath the vents (see Additional file 1: Text S2), Eq. 5 yielded at least 700 °C of T , which is too high for the shallow region. In addition, if the thermo-elastic effect is dominant, geomagnetic field changes due to thermomagnetic effects should manifest around the deformed area (Rikitake and Yokoyama 1955;Hashimoto et al. 2002), but geomagnetic data near the Jigokudani Valley showed no significant change (JMA 2017). Thus, the thermo-elastic contraction does not seem to be a major contributor to the observed deflation.
Depressurization due to fluid mass migration can be a dominant mechanism of the shallow deflation as suggested in Narita and Murakami (2018). We assumed that all the fluid within the shallow reservoir discharged to the surface, but it now appears possible that the water extruded from the source flows laterally and diffuses outward within the edifice. Gravity observation might give us insights about the migration process. Note that meteoric water influx into the shallow reservoir is likely to contribute the shallow mass balance, as reported at other volcanoes (e.g., Miller et al. 2018; Calahorrano-Di Patre et al. 2019), although we could not evaluate the quantitative contribution of the influx. Up to July 2017, there was no report of gravity changes, which might be related to water migration in the edifice including meteoric water influx into the shallow reservoir. However, to completely exclude this possibility, repetition of detailed gravity measurements are necessary.

Implications for the thermodynamic state of the deflated reservoirs
A plausible mechanism of the deflation of both shallow and deep sources can be depressurization due to mass discharge from both sources. If β c , Δm and ΔV in Eq. 2 are known, then we can constrain the other unknown parameters, that is, thermodynamic parameters of water in the deflated reservoirs: fluid bulk density ρ and fluid bulk compressibility β f . However, this situation is different between the shallow and deep reservoirs. For the deep reservoirs, we only know m , and it is extremely difficult to constrain the other parameters in Eq. 2. In contrast, the values of β c and V of the shallow reservoir are known, and m is also constrained as below: Substituting Eq.
(2) into Eq. (6) yields where β = 1 + β f /β c . Using values of V ( 7 × 10 5 m 3 ) and Q out ( 1.3 × 10 10 kg ), the maximum value of ρβ becomes 6 × 10 3 kg/m 3 . The ρ and β f in the ρβ depend on thermodynamic parameters of water such as pressure, temperature and phase state. Among them the β f can greatly vary depending on the phase state, that is, a single vapor, single liquid, or two-phase liquid-vapor state (Grant and Sorey 1979). If we can estimate β f , we may be able to obtain some information about the phase state of the reservoir (Mastin et al. 2009;Kozono et al. 2013;Hreinsdóttir et al. 2014;Juncu et al. 2017Juncu et al. , 2019. To estimate the phase state of the shallow reservoir, we consider the following depressurization scenario. During the 2014 Mount Ontake eruption, new eruptive craters were formed over the area where no fumarolic activity had been identified before, and the vigorous plume emission started there at the onset of 2014 event. The shallow deflation also started just beneath the craters. We suggest that the shallow depressurization was triggered by the formation of new conduits, which suddenly formed (6) �m < 0.3Q out .
(7) ρβ < 0.3Q out �V , a pathway connecting the shallow reservoir to the surface. The shallow reservoir had not been hydraulically connected to the surface before the eruption, and it was a confined aquifer under lithostatic pressure. After the hydraulic connection was established, the depressurization of the reservoir began, and continued until its pressure P , reached the hydrostatic pressure (e.g., Ueda et al. 2018). Assuming that host rock density is 2000 kg/m 3 and water density is 1000 kg/m 3 , lithostatic pressure should be 10 MPa and hydrostatic pressure should be 5 MPa at a source depth of 500 m, with pressure drop P of 5 MPa. Under this scenario, we calculated possible values of ρβ with a P of 5 and 10 MPa, assuming various combinations of thermodynamic parameters related to ρ and β f ( Table 3). As a starting point, we used 1.4 × 10 −8 -10 −9 Pa −1 for β c , which was derived from deliberations on reasonable settings of hydrothermal shallow reservoir (see Additional file 1: Text S2 for the detail). These values should be regarded as time averaged from the end of September 2014 to the end of July 2017. For the calculation of thermodynamic parameters, we used the IAPWS Formulation 1995 (Wagner and Pruss 2002). In order to calculate β f at single phase (liquid and steam), we used equation of its definition: 1 ρ ∂ρ ∂P H at constant enthalpy H. For two-phase liquid-vapor state, we use the formulation suggested in Grant and Sorey (1979): where P is pressure, of which unit is bar, φ is porosity of host rock (0.1-0.5), S is volumetric liquid fraction (0.1-0.9) (see Table 3 for other parameters). This equation assumes that temperature drop of geothermal reservoir in porous rock with depressurization along boiling curve at constant enthalpy causes vaporization due to heat supply from surrounding porous rock to the liquid phase. We calculated ρβ values using various likely combinations of φ , S and β c values at hydrostatic and lithostatic pressure.
The calculated values of ρβ are for constant pressure, either 5 MPa (hydrostatic) or 10 MPa (lithostatic). Figure 9 shows the calculated values of ρβ for each pressure changing porosity, volumetric liquid fraction and compressibility of the deformation source. Note that a gap exists between the ρβ values of two-and single-phase, and this is because vaporization of saturated water for two-phase states with small pressure drops contribute to the larger compressibility more than that of single-phase gas (Grant and Sorey 1979). Although we did not carry out computation with tracing trajectory path of depression from 10 to 5 MPa, for two-phase states (S = 0.1-0.9), every ρβ value during depression should range in the region between hydrostatic and lithostatic pressure boundaries (gray-shaded areas (8) (1 − φ)ρ r C r + φSρ l C l φ 1.92 × 10 −6 P −1.66 , Table 3  Volumetric heat capacity of the wetted rock (1 − φ)ρ r C r + φSρ l C l Grant and Sorey (1979) ρ r Host rock density 2000 kg/m 3 -C r Isobaric specific heat of host rock 1 kJ/kg/°C -C l Isobaric specific heat of liquid water 5.0-6.1 kJ/kg/ °C a β c Compressibility of the pressure source 1.4 × 10 −8 -10 −9 Pa −1 c in Fig. 9). Thus, most of the estimated ρβ for the twophase state at any parameter combination exceeded the upper limit value of 6 × 10 3 kg/m 3 (a dashed black line in Fig. 9), except one case (indicated as a circle in the figure). This region corresponds to a case with φ = 0.5 , β c = 1.4 × 10 −8 Pa −1 and S = 0.1 , which is vapor-dominated two-phase state, but such a reservoir condition may be unlikely to be sustained after an eruption. As demonstrated by Ingebritsen and Sorey (1988), an impermeable structure, such as a cap rock of clayrich mineral, prevented the vapor from escaping vertically and laterally toward outside of the reservoir. This structure is essential to develop the vapor-dominated two-phase reservoir within the mountain edifice. However, the impact of the eruption might produce connected fractures in a permeable hydrological structure around the reservoir, likely resulting in low permeability conditions being destroyed. Thus, we consider the vapor-dominated two-phase state is rather unlikely. For single-phase, either liquid or steam state (S = 0.0 and 1.0), the ρβ values at any parameter combination are consistently several orders of magnitude smaller than the upper limit. This means both single-phase states are allowable from the constraint on ρβ , though the vapor state can be discarded. Note that this conclusion can be valid even if we use another model to calculate ρβ , a b Fig. 9 Calculated ρβ values for the shallow reservoir. Plots are made for each of two end members of β c : a β c = 1.4 × 10 −8 Pa −1 ( V = 10 7 m 3 ) and at b β c = 1.4 × 10 −9 Pa −1 ( V = 10 8 m 3 ) (see also Additional file 1: Text S2). ρβ are computed as a function of S, (0.0 ≤ S ≤ 1.0) at each pressure. Details of these parameters setting are described in Additional file 1: Text S1. Solid and open symbols correspond to two-phase and single-phase state, respectively. Solid circles and triangles correspond to host rock porosity φ = 0.1 and 0.5 in the two-phase state, respectively. Red and blue symbols correspond to conditions of lithostatic and hydrostatic pressure at the 500-m depth, respectively. A black dashed line shows an upper limit value of ρβ , which is derived from observations. Gray-shaded area indicates possible range of allowable ρβ values, within which even calculated values with pressure drop should range which assumes the deflated poroelastic reservoir (Additional file 1: Text S3).
To maintain the single-phase steam state at a pressure of 5-10 MPa, high temperatures (> 264-311 °C) are required. However, there are several suggestions that the temperature at shallow part of Mount Ontake is relatively low. Temperatures at a shallow depth of hundreds meters beneath the vents inferred from a mineralogical study were relatively low (< 250 °C) (Minami et al. 2016). Plume temperature after the eruption was 70-100 °C JMA 2017), which implies little involvement with gas reservoirs at high temperature (> 300 °C). Asai et al. (2006) revealed the presence of a groundwater table more than 2000 m above sea level, which is close to the elevation of the shallow deflation source. Geo-electromagnetic observations estimated the low resistivity layer near the depth of the shallow source (Abd Allah and Mogi 2016), which suggests the presence of a liquid-rich zone at the depth of the shallow source. Considering these observations, a dry and superheated vapor reservoir is unlikely to exist at the shallow depth of 500 m below the surface. We suggest that a single-phase, relatively cool liquid state is the only possible thermodynamic state of the shallow reservoir.
In geothermal drilling, the state of liquid-dominated geothermal reservoirs can change to partial two-phase states in small regions near the borehole due to decompression (Grant and Bixley 2011). However, at Mount Ontake, we inferred that the state of the reservoir after the 2014 eruption was a totally liquid-dominated reservoir. If we take insights from drilling into consideration, we can infer that within the shallow deflated reservoir, a small part of the shallow reservoir near the vents exists in a two-phase state. However, the two-phase region is small, with little impact on the bulk compressibility of the reservoir. In addition, most of the volume of the reservoir may be filled with unsaturated single-phase liquid.
Liquid-rich reservoirs at a shallow depth of volcanic edifices have the potential to induce various hazards, such as hydrothermal eruptions, phreatic eruptions due to sudden decompression with landslide (Ohba 2011;Procter et al. 2014) and snow-melting lahar (e.g., Uesawa 2008). In the past, landslides have occurred at Mount Ontake, including one associated with the 1984 Nagano earthquake. Strong seismic vibration can cause large-scale landslides and subsequent decompressioninduced hydrothermal eruptions. If a large amount of fluid emission occurs in winter, snow-melting lahar may follow. Furthermore, a liquid-rich reservoir favors occurrence of more violent phreatic or hydrothermal explosions while a vapor-rich one can induce less violent explosions (e.g., Mayer et al. 2015). Thus, understanding the reservoir state is crucial to evaluating the risk of such hazards. While the estimation of the reservoir state based on a combination of ground deformation and discharge mass data has been made in several geothermal areas (Hreinsdóttir et al. 2014;Juncu et al. 2017Juncu et al. , 2019, there are few case studies regarding active volcanoes. It is necessary to increase detailed studies for other active volcanoes with high risks of phreatic eruption and landslide.

Depth of the deep deflation source
Here, we discuss the depth of the deep deflation source that is a main mass source of the emitting plume. Since the deep source likely corresponds to part of the hydrothermal system, which has been suggested by geophysical (Nakamichi et al. 2009;Kato et al. 2015;Takagi and Onizawa 2016), geochemical ) and mineralogical (Minami et al. 2016) studies, the deep deflation may reflect mass emission from the hydrothermal system. Deformation signals from the deep deflation source were observed only in the TO baseline. The insufficiency of the wide-range geodetic data makes estimating the detailed source model (e.g., location, geometry and volume change of the source), extremely difficult. There are several candidates of this deeper source. Takagi and Onizawa (2016) identified two inflation sources at 1 km and 3 km depths beneath Jigokudani Valley during the 2007 eruption, and Narita and Murakami (2018) concluded that the source at 3 km depth was plausible as the deep source. Murase et al. (2016) modeled co-eruptive subsidence at the east flank using leveling data as horizontal crack emplaced 3 km below sea level, (i.e., 6 km below the eruptive craters) (Additional file 1: Figure S3). To observe whether this model could account for the GNSS data, we tried applying this model to the TO baseline data. We searched the best-fit value of the crack closure while fixing the crack geometry and the location. While a possible solution may have been found, it should be regarded as only one possible candidate and far from conclusive (Additional file 1: Figure S4). Hence, the results from this study and Narita and Murakami (2018) imply that the possible pressure source might be emplaced from 3 to 6 km below the craters. This range of depth is almost consistent with the hypocenter zone of VT seismicity after the 2014 eruption, and there is no seismicity at greater depth. This seismicity suggests low temperature below brittle-ductile transition (~ 400 °C) at 3-6 km depth. Furthermore, this zone overlies the dike intruded in the 2007 eruption (Fig. 10), which may allow for the development and pressurization of a volcanic hydrothermal system over the dike. Narita and Murakami (2018) inferred that the shallow source might play a minor role in supplying water for the eruption, based on lack of co-eruptive deflation corresponding to fluid discharge with violent ballistics emission (Yamada et al. 2015). Our results may support this inference because post-eruptive temporal behavior of the shallow source clearly differs from that of the deep source and the plume emission, (e.g., time functions of decay and relaxation times). This difference may be attributed to how the sources contributed to the release of violent fluids at the time of the 2014 eruption (Fig. 10). The fluids were likely discharged from the deeper source due to release of overpressure accumulated before the eruption. The built-up process of the overpressure prior to the eruption might have corresponded to the signal detected by stacking various combinations of GNSS baseline data (Miyaoka and Takagi 2016). It may also be related to stress field changes (Terakawa et al. 2016). As a result of the fluid release at the time of the eruption, the deep source may rapidly depressurize for an early period of time after the eruption (approximately 2 months). Subsequently, steady deflation of the deep source was followed and continued until July 2017. This steady deflation may correspond to an early stage of exponential decay with an extended relaxation time of several decades, which was observed at other volcanoes, (e.g., 20 years ground deflation after the 1995 eruption at Mount Kuju) (Hashimoto et al. 2002).

Role of the two deflation sources in the 2014 phreatic eruption
While the present status of the hydrothermal system at Mount Ontake is an open system, it will likely shift to a closed system due to several factors, including the progression of mineralogical sealing. This will initiate the re-pressurization of the edifice (Christenson et al. 2010;Hamling et al. 2016;Tanaka et al. 2017). High-frequency InSAR observations and continuous near-field monitoring will allow for detecting ground inflation preceding phreatic eruptions such as the 2015 Hakone eruption (Kobayashi et al. 2018), which will provide important insights about the pressure conditions within active volcanoes. This information will be also useful for hazard potential evaluation.

Conclusions
Quantitative investigations of the relationship between the shallow and the deep deflations and the plume emission after the 2014 eruption at Mount Ontake revealed that the main source of water-vapor plume is not a shallow source detected by InSAR, but a deeper deflation source detected by GNSS data. Furthermore, this deep source approximately 3-6 km below the craters might play a key role in supplying fluid for phreatic eruptions. In contrast, the shallow deflated source only appears to have minimal involvement in the eruption process. This shallow reservoir was assessed to be a liquid-dominated reservoir under thermodynamic consideration based on the estimated mass balance. Although such a wet reservoir at a shallow depth was not likely to be a main water source of the eruption, it still has the potential to contribute to future eruptions or other hazardous phenomena, such as snow-melting lahar. To assess the pressure condition within the volcano and to evaluate potential hazards, continuous remote and near-field observations are necessary in the future.
Additional file 1: Text S1. The Effect of fluid compressibility on time function of the plume discharge. Figure S1. Time function of ∆m(t) when considering temporal change in ρβ (dashed lines) for combination of different values of source compressibility (β c ), porosity (φ) and volumetric liquid fraction (S). These are plotted for both of independent satellite paths.
Solid line corresponds to ∆m(t) without temporal changes in ρ(t) and β f (t) which means that fluid is incompressible. ∆m(t) on the vertical axis is normalized by its maximum value. τ min shows the smallest value of the relaxation time among the plotted curves in each figure. Text S2. Realistic value of Compressibility of Deformation source ( β c ). Text S3. Effect of poroelasticity on estimate of thermodynamic state at the shallow reservoir. Figure S2. Calculated ρβ p values for the shallow reservoir using the formulation of Juncu et al. (2019). Blue and red colors indicate hydrostatic and lithostatic pressure condition at 500-m depth, respectively. Circles and triangles correspond to reservoir porosity of 0.1 and 0.5, respectively. Solid and outlined symbols show two-or single-phase state, respectively.  Figure S4. Time series of length change in the four baselines from 28 September 2014 and length change predicted by the crack model of Murase et al. (2016). Red, yellow and blue bars correspond to crack closing of 0.5, 1.0, 1.5 m, respectively. Miyaoka and Takagi (2016) highlight that Otiai-Outaki baseline includes non-negligible amount of tectonic deformation, which appears as compressional change along the direction of northwest and southeast. We evaluated and subtracted the tectonic trend included in Otiai-Outaki baseline from the raw data of Otiai-Outaki baseline using GEONET data around Mount Ontake. We chose the GEONET data that include no volcano-related and seismicity-related deformation, only tectonic deformation, from 2012 to 2018. Table S1. The source model of Murase et al. (2016).