Precursory ground deformation of the 2018 phreatic eruption on Iwo-Yama volcano, revealed by four-dimensional joint analysis of airborne and spaceborne InSAR

We present detailed maps of local-scale 3D deformation preceding the 2018 phreatic eruption at Iwo-yama volcano (south of Kyushu Island, Japan), using a combination of airborne and spaceborne Interferometric Synthetic Aperture Radar (InSAR) data. The 3D and 2.5D deformation maps obtained at different periods allow us to successfully track their spatiotemporal evolution and to infer the transition of subsurface conditions responsible for the precursory deformation observed from 2014 to 2018. From 2014 to 2016, ground inflation depicted an axisymmetric pattern with the maximum displacement at the center of the deformed area. However, from 2016 to 2018, an inflation peak moved to the southern edge of the area deformed during 2014–2016 and became more localized, which was close to the newly generated vents in the 2018 eruption. Modeling of the inflations suggests that pressurization within a crack at a depth of 150 m beneath the Iwo-yama geothermal area caused the 2014–2016 deformation and had continued until the 2018 eruption. Modeling results highlight the persistence of the local ground inflation pattern just above the southern edge of the crack, which suggests the presence of a shallower inflation source contributing to the local inflation. Consequently, we interpret the sequence of these deformations as follows: from 2014, deeper-rooted fluid started to inject into a fluid-saturated crack at 150-m depth, which caused the 2014–2016 deformation. Then, after 2016, the crack inflation continued because of the continuous fluid injection and formed another pressurized part directly above the southern tip of the crack. Additionally, the results of the time-series analysis of the satellite InSAR data revealed that the local inflation started around April 2017 for which thermal activity including a mud emission became pronounced around the location of the local inflation. As a result of an episodic increase in supply rate of magmatic fluids from a deep magma reservoir from early 2018, a phreatic eruption finally occurred in the vicinity of the most deformed point, providing a clue for predicting future eruption sites, as was also observed in the Hakone 2015 eruption.


Introduction
High spatial resolution and no requirement of groundbased instruments for satellite-based Interferometric Synthetic Aperture Radar (InSAR) make it a standard tool for studies on volcanic deformations caused by deeper-rooted magmatic activities (e.g. Wicks et al. 2002;Miyagi et al. 2013;Pinel et al. 2014;Sigmundsson et al. 2015;Lundgren et al. 2017) and shallow hydrothermal systems (e.g., Hamling 2017;Doke et al. 2018;Kobayashi 2018;Narita and Murakami 2018). One of the few drawbacks of spaceborne InSAR is the lack of sensitivity to the North-South component of deformation due to inescapability from the near-polar orbit, which makes retrieval of 3D deformation challenging. Although the Multiple Aperture Interferometry (MAI) and pixel offset methods enable us to retrieve 3D deformation (e.g., Tobita et al. 2001;Bechor and Zebker 2006;Jo et al. 2015), the spatial resolution and precision of these methods are limited. These deficiencies prevent us from applying these techniques for retrieval of 3D components to small deformation both in scale (< 1 km) and magnitude (< 10 cm). Recent studies using interferograms acquired from more than three independent looking directions have demonstrated their effectiveness in the mapping of 3D deformation around volcanoes for larger magnitude signals (> 10 cm) and spatial scale (> several km) (e.g., Lundgren et al. 2013;Delbridge et al. 2016;Morishita et al. 2016;Narita and Murakami 2018).
At Kirishima Iwo-yama volcano, in southern Japan, local-scale deformation with a maximum displacement of ~ 15 cm has been observed over an area of ~ 500 m from 2014 to 2018 (Meteorological Research Institute 2018). Following this deformation, a phreatic eruption occurred in 2018. Because the magnitude and spatial extent of the deformation were relatively small, only satellite SAR data are insufficient to extract the 3D deformation. However, flexibility of flight paths and higher resolution of airborne SAR can mitigate this drawback and allow us to obtain 3D deformation with high precision even if the deformation signal is small (e.g., Lundgren et al. 2013;Delbridge et al. 2016). This study attempts to estimate high-quality 3D deformation at Iwo-yama using a couple of airborne and spaceborne InSAR datasets; repeated observations make four-dimensional analysis possible.
Iwo-yama volcano, a part of Kirishima volcano group, south of Kyushu Island (Fig. 1a), is located in the eastern part of Ebino plateau and near the northwest flank of Karakuni-dake (Fig. 1b). Iwo-yama was formed during the sixteenth-seventeenth century with an extrusive lava flow eruption (Tajima et al. 2014). The latest eruption was a phreatic eruption that occurred at the Iwo-yama east crater in 1768 (Fig. 1c), after which the volcano has been dormant for ~ 250 years (Tajima et al. 2014).
Ebino plateau is rich in thermal activity with numerous fumaroles and hot springs since the nineteenth century (Funasaki et al. 2017). From the 1930s to the 1960s, thermal activities were characterized by the heat discharge from numerous boiling mud pots, hot springs and fumaroles over Ebino plateau, including Iwo-yama (Funasaki et al. 2017). In the 1970s, fumarole activity in Iwo-yama was temporally intensified; the heat discharge rate of the fumarole plume was up to 61 MW, and the maximum height of the plume was up to 700 m (Kagiyama et al. 1979). However, the thermal activity started to decay in the 1990s, accompanied by local subsidence over Iwo-yama with a diameter of ~ 500 m detected by JERS-1 InSAR (Ozawa et al. 2003), and the main fumaroles disappeared in 2007 (Funasaki et al. 2017). Since December 2013, microseismicity has been occurring beneath Iwo-yama, followed by volcanic tremors and tilt changes since 2014. In 2015, precise leveling surveys and PALSAR-2 InSAR revealed ground inflations around Iwo-yama (Meteorological Research Institute 2018; Kyushu University 2019), which could be caused by the pressurization of a shallow hydrothermal system. The temporal evolution of the ground inflation and thermal activity was rather complicated. The activity was characterized by a repetition of rise and fall in fluid injection rates into the ground surface; cycles of uplifts and subsidence and of excitation and rest of thermal activity were observed. In addition to the progress of the ground deformation, thermal activity also showed temporal changes in its spatial distribution from December 2016 to April 2017. From 2014 to December 2016, the area of geothermal anomalies was located within the Iwo-yama edifice (Fig. 1c). From December 2016 to April 2017, the geothermal area extended to the south of the edifice (Fig. 1c), and new fumarole vents (including "H") and mud pots (including "F") were formed, and a mud emission occurred at jet fumarole vent A. Following these activations of the deformation and the thermal activity, a phreatic eruption finally occurred on April 19, 2019, forming new eruptive vents in the south of Iwo-yama (Fig. 1c). The erupted mass was only 1500 tons (Nagai et al. 2018), which tends to be smaller than other phreatic eruptions; 0.89-1.2 million tons in the Ontake 2014 eruption (Takarada et al. 2016), 20,000 tons in the Kuju 1995 eruption (Nakada et al. 1996), and 0.23 million tons in the August 2012 Te Maari eruption, which is converted from an erupted volume of 230,000 m 3 (Turner et al. 2014) by assuming an ash density of 1000 kg/m 3 . Another smaller eruption followed on April 26, 2019, forming new vents to the 400-m west of the Iwo-yama edifice (Fig. 1b).
The purpose of this study is to trace transient conditions in the shallow hydrothermal system of Iwo-yama, leading to the 2018 eruption, by examining a detailed 3D deformation field and tracking its spatiotemporal evolution. To achieve this, we utilized combinations of L-band airborne SAR (Pi-SAR-L2) and L-band spaceborne SAR (ALOS-2/PALSAR-2) data, both of which are operated by JAXA. In addition to the 3D deformation field, we additionally conducted a timeseries analysis using only PALSAR-2 data, which have a high temporal resolution of 1-2 months over Iwo-yama to compensate for the temporal resolution.

Data and interferometry
PALSAR-2 data were acquired from both ascending and descending paths; whereas, Pi-SAR-L2 data were acquired from three flight paths with headings from SE to NW (path EW), from NNE to SSW (path NS) and from SSW to NNE (path SN) (Table 1). We applied standard procedures of SAR interferometry to all the data except for additional processes for Pi-SAR-L2 data. We used the RINC software (Ozawa et al. 2016) for the entire process. We multilooked the images in the range and azimuth directions by 4 and 4 in PALSAR-2 and 5 and 25 in Pi-SAR-L2 data, respectively. This multilooking corresponds to a pixel resolution of approximately 10 m. We used the 10-m mesh Digital Elevation Model (DEM) of the Geospatial Information Authority of Japan (GSI) to remove topographic fringes. For unwrapping, we used the SNAPHU software (Chen and Zebker 2002). We reduced atmospheric noise by removing elevationcorrelated components calculated over the area except for the deformed area. Finally, geocoding was conducted so that the ground spacing was 10 m.
Pi-SAR-L2 data require several additional procedures to achieve higher coherence and to reduce the noise inherent in airborne InSAR. This noise originates from the geometrical heterogeneity of the flight trajectory of the airplane. The spaceborne SAR has a simple orbit because there are no perturbations due to atmospheric drag, so that it can be readily approximated by a simple quadric function of spatial location in the range and azimuth directions. Thus, co-registration of master and slave images for the spaceborne InSAR is straightforward. In contrast to the case of spaceborne InSAR, the speed and flight path of the airplane for airborne InSAR are constantly perturbed by random wind changes. In particular, the change in speed in the azimuth direction severely degrades the coherence in airborne InSAR. Since the observation flight trajectories selected over 11,000-m height is still heterogeneous due to the atmospheric behavior (Shimada et al. 2013), SAR image generation initiated even with the motion compensation method using the inertial navigation system leaves the fraction of doppler modulation in a SAR image and the low-frequency phase other than the geophysical information needs to be extracted using the accurate co-registration method. Unless we take proper care of this image distortion, it can cause significant coherence degradation. To cope with this problem, we precisely estimated the amount of the offset between master and slave images with a very high spatial sampling rate and then interpolated the offset distribution using cubic Bezier curves. The interferograms, thus, generated still have spatially long-wavelength fringes caused by the inaccurate flight position measured by the onboard navigation device. Because this study focuses on localized (~ 500 m) deformation around Iwo-yama (Meteorological Research Institute 2018), we removed the long-wavelength phase error using a simple spatial high-pass filter. This filtering process is justified by the fact that no long-wavelength deformation is found in the independent PALSAR-2 dataset. Thus, we obtained phase changes within a small spatial extent (< 500 m) by Pi-SAR-L2 from three different directions.

Procedure of time series analysis
Here, we briefly explain the procedure of the PALSAR-2 dataset time-series analysis. We employ the New Small Baseline Subset (NSBAS) method (e.g., López-Quiroz et al. 2009), implemented in GIAnT software (Agram et al. 2013). This allows us to estimate the time series of line-of-sight displacements (dLOS) by excluding the effect of DEM error, using the PALSAR-2 data acquired from the two orbits (Table 1). To evaluate the estimated error of the dLOS time series, we implemented the Jackknife resampling approach, which enables calculating the standard deviations of the estimated displacement between two adjacent acquisition times at each pixel, and then calculates the estimated error of the cumulative displacements by error propagation. To maintain high temporal resolution of the time series, we set a temporal Gaussian filter with a length of only 0.05 year, which creates weak temporal smoothing. The selection of interferograms used for the NSBAS analysis is based on the following criteria: (i) the perpendicular baseline is less than 300 m (Fig. 2), (ii) the mean value of spatial coherence over an interferogram is larger than 0.3, and (iii) the interferograms are not affected by heavy atmospheric noise. To handle criterion (iii), we found noisy interferograms as follows: first, the standard deviation (σ) of each interferogram was calculated over the non-deformed region. Then, the mean σ value was calculated among all interferograms. Finally, we regarded interferograms with more than twice the mean σ value as outliers and excluded them from the time-series analysis. In addition, we extracted temporally coherent pixels, that is, the pixels at which the temporal mean value of the coherence is larger than 0.3.
To select a reference pixel (zero-displacement pixel), we employ the following criteria suggested in Yunjun et al.  disturbance, (ii) temporally coherent, and (iii) close to the deformed region, to avoid spatially correlated components. Criterion (i) was satisfied by the interferogram selection. To satisfy (ii), we selected a highly coherent pixel with a temporal mean coherence value of > 0.4 (Fig. 3). To satisfy (iii), candidate pixels are located within a distance of 2 km from Iwo-yama. As a final step, we tested several candidates of reference pixels to iterate the NSBAS analysis. Among these candidates, we finally selected a pixel resulting in a less noisy time series as an appropriate reference pixel (Fig. 3).

Estimation of 3D and 2.5D deformation field
We computed 3D deformation on a pixel-by-pixel basis based on a conventional weighted-least-squared inversion (Wright et al. 2004). As the inversion weight, we used phase variances calculated over the non-deformed area in each interferogram (Table 1) (Table 1). We estimated the 3D deformation using interferograms between 2014 and 2016 and between 2016 and 2017, obtained from a total of five independent orbits (Table 1). Because no Pi-SAR-L2 data were available from 2017 to 2018, we alternatively estimated the 2.5D deformation of quasi-vertical and quasi-eastward components (Fujiwara et al. 2000) using the ascending and descending PALSAR-2 data. While InSAR images acquired during the same period should ideally be used for a precise calculation of the 3D deformation, the master date of the descending PALSAR-2 path (January 2015) was approximately 3 months after that of the ascending path (end of September 2014), which can give inaccurate results. As discussed later, however, this offset had a negligible effect because no significant deformation was observed during the 3 months between September 2014 and January 2015.

Interferograms
All the interferograms from 2014-2016, 2016-2017, and 2017-2018 show the dLOSs toward the satellite and airplane ( Fig. 4) in an area with a diameter of a few hundred meters, suggesting ground inflation. The airborne interferograms from the SN and EW paths between 2014 and 2016 (Fig. 4a, c) have incoherent lines, likely due to fluctuations in the flight path caused by sudden changes in wind speed. These incoherent areas do not affect the 3D decomposition of the deformation around Iwo-yama because they are separated from the deformed region by ~ 500 m.

3D and 2.5D deformation field
We estimated 3D deformation during 2014-2016 and 2016-2017 and 2.5D deformation during 2017-2018, using the interferograms presented in Fig. 4. The 3D deformation field between 2014 and 2016 indicates axisymmetric inflation with its axis at the center of the Iwo-yama fumarole area (Fig. 5a-c). The maximum uplift is located at the center of the contour lines of the uplifted region (Fig. 5a). This deformation field is likely to be modeled by an axisymmetric source with simple geometry, such as a spheroidal source. A general feature of the deformation field from 2016 to 2017 is similar to that from 2014 to 2016, with minor differences ( Fig. 5d-f). For example, the location of the uplift peak is clearly shifted southward by ~ 70 m: moved from the peak of the 2014-2016 deformation field (hereafter, peak A; the red circle in Fig. 6) to that of the 2016-2017 deformation field (hereafter, peak "B"; the blue circle in Fig. 6). The displacements of the EW and NS components were also concentrated around peak B (Fig. 5d-f), which are very close to the 2018 eruption vents. The extent of the deformed area of this local inflation around peak B is only ~ 100 m, indicating a shallowing of the pressurization compared with the depth of the 2014-2016 inflation source. These deformation features with different scale ground inflation, observed during 2014-2016 and 2016-2017 suggest that a single inflation source is responsible for the 2014-2016 deformation and subsequently continued until at least 2017 to produce peak A (Fig. 6a, b). In addition, from 2016 to 2017, a shallower source became active and produced the local inflation around peak B, which is 70 m to the south of peak A (Fig. 6b, c). Another important feature of the deformation around peak B is more asymmetric than that around peak A, that is, the westward  On April 15-18, 2017, a new fumarole vent H with vigorous fumarole discharge was formed, which is expressed as "jet fumarole" in Tajima et al. (2020). On April 26, 2017, the mud emission created jet fumarolic vent A with a diameter of 1.5 m (Tajima et al. 2020). These imply that increasing fluid injection to shallower levels promotes fracturing of the host medium at near-surface depths. By numerical modeling of progressive inflation of a sill-like source accompanying the host-rock fracturing, Holohan et al. (2017) show that as the fracturing progresses, finite and plastic, rather than elastic, deformation becomes dominant, and the deformation pattern tends to exhibit significantly asymmetric one. Focusing on the property of the fracture-rich host rock and the axis-asymmetric horizontal displacements at Iwo-yama geothermal area, it is highly possible that the localized deformation is probably no longer elastic, but rather inelastic. We will further mention the possibility of the inelastic deformation during this period, later in the modeling section for the period of 2016-2018.
From September 2017 to April 2018, just before the eruption, the quasi-vertical displacement pattern was similar to that from 2016 to 2017, suggesting ongoing ground inflation from 2014 around peak A and no migration of the peak location during this period (Figs. 5g,h;6c). Note that this deformation field of quasi-vertical and quasi-eastward components contains a relatively large noise component, which may be attributed to contamination of the atmospheric disturbance (Fig. 4l). However, we can still recognize the general characteristics of the deformation and the above interpretation remains valid.

Source modeling
The deformation sequence described above can be summarized as follows: from 2014, pressurization started probably at a shallow depth, and from 2016, another source began to be pressurized at an even shallower depth near the southern edge of the previously deformed area during 2014-2016, and then the inflation of both pressure sources continued until the 2018 eruption.
We use forward models of a sill-like (Okada 1985), spherical (Mogi 1958), and dipping prolate source (Yang et al. 1988) in an elastic, homogeneous and isotropic half-space medium. To invert for the source parameters, we employed Geodetic Bayesian Inversion Software (Bagnardi and Hooper 2018), which enables us to perform an efficient search of the optimal solution within certain ranges of model parameters by applying the Markov Chain Monte Carlo method.
For source modeling, we first focus on the 2014-2016 deformation, which is characterized by its axisymmetric pattern. From 2014 to 2016, best-fit sources are emplaced at a shallow depth of 100-150 m with an inflation volume ( V ) of 4 − 8 × 10 3 m 3 ( Table 2). The estimated sources lie above a low-resistivity region, which was interpreted as a clay-rich impermeable layer (Tsukamoto et al. 2018). The best-fit source model is a near-vertically dipping prolate ellipsoid (Fig. 7j-l). However, its pressure change ( P ) value is too large to be realistic, as we show below. In the following discussions, we assume that the rigidity ( µ ) of the host rock is as low as 0.1-1.0 GPa, a realistic value for hydrothermally altered volcanoes and geothermal fields at shallow depths (e.g., Lynne et al. 2013;Juncu et al. 2019). The estimated value of �P/µ (2.23) requires P to be more than 230 MPa, which greatly exceeds the tensile strength of hydrothermally altered rock (1-10 MPa, Kumar et al. 2011). Similarly, the possibility of a spherical source can be rejected. P of the Mogi source is expressed as �P = �V µ/ r 3 π , where r is the radius. Assuming r = 20 m, which is the upper limit derived from the Mogi source assumption ( r/d < 0.2 where d is the source depth), P becomes 17-170 MPa, which is also physically unrealistic. Thus, only a sill-like source can be a possible candidate. The P value of the sill source can be 3.4-34 MPa, based on the relationship �P = µ�u/(L(1 − ν) (Gudmundsson 1990), where u is opening (1.7 m), L is the mean size of the sill (68 m), and ν is Poisson's ratio (0.25). When µ is 0.1 GPa, ΔP is likely to be 3.4 MPa. Thus, inflation of a sill-like crack can be the most physically plausible source when the host rock consists of compliant and fracture-rich material.
Note that the sill model fits the observations more poorly than other models (Table 2); the modeled horizontal displacement is significantly smaller than the observations (Fig. 7d-f ). The heterogeneity of the surrounding medium can be a plausible explanation for this inconsistency. Considering a simple two-layer medium in which the rigidity of the upper layer is smaller than that of the lower one, the ratio of horizontal to vertical displacement on the ground surface can be larger than that in a homogeneous medium, even if the same crack geometry is applied (Fialko et al. 2001). Thus, inversion of this distorted ground displacement pattern may prefer a spherical or prolate source to the true solution with a sill. The near-surface medium of Iwo-yama is composed of soft layers, such as a fracture-rich and hydrothermally altered lava dome and deposits of pumice, pyroclastic flows, and the sector collapse of Karakuni-dake (Geshi and Kobayashi 2016), suggesting the high possibility of selective amplification in horizontal displacement due to heterogeneity of the medium. Therefore, we regard the sill-like crack as a plausible model.
In recent studies that target volcanic deformations observed at shallow volcanic hydrothermal systems, sill-like crack sources have frequently been applied as a first-order approximation of a fluid-saturated layer with pore-fluid pressure change (e.g., Hamling 2017; Kobayashi 2018; Kobayashi et al. 2018;Miller et al. 2018); this study does the same. However, a more realistic physical process of these cases is thermo-poroelastic deformation, which has been better modeled by a coupling analysis of hydrothermal fluid flow and ground deformation (e.g., Currenti and Napoli 2017). Although this approach is also essential in the case of Iwo-yama, no mass and heat discharge of fumaroles have been estimated, which is necessary to assign heat and mass flow boundary conditions and to validate the fluid-flow model. More realistic

Table 2 Parameters of best-fit source models (August 2014-August 2016)
Square brackets show 2.5 and 97.5 percentiles of posterior probability density functions of the parameters 1 Residual is defined as: where Obs is the observed displacement, Cal is calculated displacement, i is the pixel index, and N is the number of pixels used for the calculation 2 Ratio of minor axis b to major axis a 3 Dip angle δ is defined as δ = 0° if the ellipsoid is horizontal and δ = 90° if it is vertical 4 Strike is defined as clockwise from the north 5 ΔP/μ corresponds to ratio of pressure change (ΔP) to host rock rigidity (μ) 6 Position of the southwest corner of the horizontal crack  Narita et al. Earth, Planets and Space (2020) 72:145 modeling of the deformation time series along with heat discharge remains as a future study.
We then model the 2016-2018 deformation under an assumption that two sources continued to be active until 2018; one is a sill-like source to explain the entire   (Table 3), which is also above or upper limit of the conductive sealing layer, similar to the depth of the 2014-2016 model. Sill size and location are almost same as that of the 2014-2016 model (Fig. 7m-u; Table 3), which indicates that the same source has consistently produced the deformation around peak A. The best-fit Mogi source is emplaced at a depth of ~ 30 m, which is greatly shallower than the sill source, and is above the Although the localized inflation observed during 2016-2018 is well reproduced by the elastic models ( Fig. 7maa), it is still possible that this deformation is inelastic because of its large strain. Since cumulative uplift during 2014-2018 is up to 0.2 m at peak B and the deformed length of the host rock is approximately 30 m, a maximum vertical axial strain is 0.007, which indicates that the local uplift is dominated by the regime of finite strain rather than infinitesimal strain. Experimental studies have shown that hydrothermally altered and fracture-rich rock with an axial strain of more than 10 −3 can deform in an inelastic behavior (e.g., Heap et al. 2015), which can support the possibility of inelastic strain around peak B.

Timing of migration of the inflation peak
Here, we examine the exact timing of migration of the inflation peak, observed during 2016-2017, based on the time series of the PALSAR-2 data. The comparison of locations of the uplift peaks during the three periods (2014-2016, 2016-2017, and 2017-2018) indicates that the peak location migrated southward by ~ 70 m from peak A during 2014-2016 to peak B during 2016-2018 (Fig. 6). Peak B is only 30 m to the west of the south vents of the 2018 eruption (Fig. 6b, c). To see the detailed time series during 2016-2017, we present the time series of dLOSs of both paths from 2014 to 2019, as sampled at peak A (the red point in Figs. 6a and 8a, b) and peak B (the blue point in Figs. 6b,c and 8a, b). From September 2014 to March 2017, the dLOSs of both paths at peak A were consistently larger than those at peak B, suggesting that the uplift at peak A was still dominant. Then, in April 2017, path 131 exhibited larger dLOSs at peak B and path 23 exhibited comparable dLOSs at both peaks (Fig. 8), suggesting that the deformation at peak B became prominent around April 2017. Note that displacements from September 2014 to January 2015 are small enough that the time difference in paths 131 and 23 used for the 3D estimation does not affect the 3D deformation field. This transition of the deformation pattern almost coincided with the increase in thermal activity around peak B (the dotted line in Fig. 8c, d), such as an appearance of mud pot F on March 19-21, 2017, the formation of jet fumarole vent H on April 15-18, 2017 (Fig. 6), and the mud emission at jet fumarole vent A on April 26, 2017. Temperature measurements at 1-m depth within the geothermal area show that the geothermal anomaly area (GAA) with more than 50 °C extended to the south of peak A, close to fumarole vent H and mud pot F, from March 2017 (Tajima et al. 2020). Deformation acceleration accompanied by intensified fumarole activity was also observed in Midagahara volcano during 2007-2010 where the acceleration of ground inflation was attributed to the increase in pore fluid pressure due to activated fumaroles overlapping on the most deformed area (Kobayashi 2018). Moreover, the apparent equilibrium temperature (AET) of the fumarole gases rapidly increased by approximately 200-300 °C from January to May 2017 (Tokai University et al. 2018). By considering these thermal activations, we infer that the prominence of the local inflation after April 2017 was triggered by the increased injection rate of hydrothermal fluid from greater depths, such as another pressure source at 700-m depth just beneath a low resistivity, likely impermeable clay-rich layer (Tsukamoto et al. 2018). This inference can also be supported by deflation of the 700-m source observed from February 2017 to October 2017 (Kyushu University 2019).

Interpretation of temporal evolution of the ground inflation
Here, we address factors controlling the temporal features of the dLOSs for the entire observation period (2014-2019) by comparison with other observations. The dLOSs of both ascending and descending tracks (paths 131 and 23) show three episodic accelerations of the ground inflation at peak B (Fig. 9); the first acceleration event is late 2015-mid 2016 with dLOS velocity of 0.8 cm/month, the second one is late 2016-mid 2017 with the velocity of 1.2 cm/month, and the last one is late 2017-mid 2018 with the velocity of 1.4 cm/month (calculated in the path 131).
The first event seems to coincide with the geothermal manifestation from December 2015 to March 2016 ( Fig. 9), including the first expansion of the GAA, the appearance of fumaroles, and an increase in spring water temperature (Tajima et al. 2020). These geothermal activations are interpreted to have initiated when hydrothermal fluids started to rise from the inflation source at 700-m depth, which is considered to be a hydrothermal reservoir just beneath the impermeable layer (Tajima et al. 2020). Based on these inferences, we suggest that the first acceleration resulted from fluid injection from the deeper hydrothermal reservoir into the crack source at 150-m depth, emplaced in a shallower part of the hydrothermal system. Because the increase in the spring water temperature during this period was only 6 °C (Tajima et al. 2020), the contribution of the thermal expansion of the crack-filled fluid to this deformation may be neglected, and the pore-pressure increase within the crack source, responding to the injection rate increase, may be the dominant deformation mechanism. The second acceleration, as mentioned in the previous section, can be caused by episodic fluid injection at a higher rate than the first one, from the deeper source. Tajima et al. (2020) inferred that the temperature of 130 °C of jet fumarole H can be attributed to leakage from a shallower part up to 350-m depth, rather than the deeper hydrothermal reservoir at a temperature of more than 200 °C (Tsukamoto et al. 2018). Considering this inference, the mechanism producing the second acceleration is that hydrothermal fluids rose from the deeper reservoir, injected into a shallower depth around the crack-source depth, and extruded fluids within the crack into the surface, which could cause the formation of jet fumarole H, the mud emission, and the rapid expansion of GAA (pink-shaded zone in Fig. 9a).
The third acceleration event, just before the eruption, seems to be correlated with the deeper-rooted magmatic activity. First, geochemical observations of the fumarole gases have suggested the contribution of magmatic fluids. The He/CH 4 ratio and AET of the fumarole gases rose rapidly from 6 months and 1 month before the eruption, respectively (the black and striped zones in Fig. 9a). These observations are interpreted as responding to the increased flow rate of magmatic, hot fluids into the hydrothermal system just before the eruption (Tokai University et al. 2018), coincident with the last dLOS acceleration. In addition, the timing of these activities coincides with an inflation of a magma chamber at 7-km depth to 4-km west of Iwo-yama (Fig. 1a) from July 2017 to January 2019 (Fig. 9b). Although the inflation time series was disturbed by co-eruptive deflation of the March 2018 Shinmoe-dake eruption (Fig. 9b), inflation continued after the eruption and stopped in January 2019 (Japan Meteorological Agency 2019), which is correlated with the dLOS sequence for the same period. Since this magma chamber is considered to connect to Iwo-yama as well as Shinmoe-dake (Aizawa et al. 2014), magmatic fluids are also supplied to the Iwo-yama due to magma supply into the inflated magma chamber, which could cause changes in gas composition and temperature and the last acceleration of the pressurization of the entire shallow hydrothermal system. Synchronized inflation between a shallow hydrothermal system and a magma chamber, leading to eruptions, has also been observed before the Hakone 2015 eruption .

Sequence of transition of subsurface conditions leading to eruption
Based on the above inferences, we summarize and interpret the spatiotemporal evolution of subsurface pressure conditions as follows (Fig. 10). After 2015, hydrothermal fluid was injected into a sill-like crack, embedded just above or in the uppermost part of the sealing layer with a partly fractured conduit zone (Tsukamoto et al. 2018). As the pressurization within the crack progressed, the stress concentration at the southern tip of the crack prompted the gradual formation of vertical fractures reaching the shallower part, and geothermal activity flourished above the southern edge of the crack from March to April 2017. It is also possible that the formation of the shallower pressurization could  Tajima et al. (2020). The red and blue circles indicate dLOS time series of path 131 and 23 at peak B (same as in Fig. 8). b Plot of cumulative volume change of the magma chamber at 7-km depth after November 2013, derived from GNSS baseline observation (JMA 2019). The red arrow indicates the timing of the Shinmoe-dake eruption with co-eruptive deflation of the magma chamber be attributed to pre-existing pathways, as suggested by precedent fumarole activity in the 1970s in the identical location as the 2018 eruptive vents (Kagiyama et al. 1979). At the end of April 2017, hot fluids rapidly rose from the deeper hydrothermal reservoir through the pathway, producing the local pressurization directly above the southern tip of the crack and episodic thermal manifestations, including the formation of jet fumarolic vent H and the mud emission at jet fumarolic vent A. Because of the rapid increase in the injection rate of magmatic fluids derived from the magma chamber at 7-km depth starting from Early 2018, the eruption on April 19, 2018, finally occurred at vents formed close to the most deformed site (peak B).
This sequence (migration of the uplift peak and eruption occurrence at the vicinity of peak) is remarkably similar to the case of the 2015 phreatic eruption at Hakone volcano, except for the temporal feature of the pre-eruptive ground inflation. In the Hakone 2015 eruption, the inflation started 10 months before the eruption onset with a constant rate of ~ 0.5 cm/month, followed by an exponential acceleration to ~ 10 cm/ month 2 months before the eruption , which might provide hints for predicting the timing of the eruption because of its simple time function. In contrast, at Iwo-yama, there were three episodic accelerations, and the last one eventually led to the eruption (Fig. 9a). Note that the speed of the last inflation was no faster than that of the others (Fig. 9a), which implies difficulty in assessing the possibility of the eruption occurrence and much more of its timing. However, the good correlation between the location of the eruption onset and the most deformed area suggests that frequent InSAR monitoring would help to predict future eruption sites. To ensure the generality of this phenomenon and its effectiveness, it is essential to make steady efforts to increase case studies of the pre-eruptive stage of phreatic eruptions.

Conclusions
This study presented small-scale 3D deformation fields preceding the 2018 Iwo-yama phreatic eruption. The 3D deformation estimation revealed clear temporal changes in the deformation fields; the most deformed site clearly migrated to the eruption vents. Modeling results suggest that pressurization of a sill-like crack at 150-m depth started in 2014, which produced an axisymmetric inflation pattern from 2014 to 2016 and that after a mud emission in April 2017, local uplift started just above the southern tip of the crack, possibly indicating the formation of another, shallower pressurized zone. As the pressurization progressed and the injection rate of magmatic fluids rapidly increased, a phreatic eruption eventually occurred in April 2018 close (~ 30 m) to the locally uplifted site. This sequence is significantly similar to the pre-eruptive stage of the 2015 Hakone eruption; however, the temporal features of the inflation processes clearly differed. In the future study, it is necessary to address the origin of this difference. Our results show that