Imaging the Hydrothermal System of Kirishima Volcanic Complex With L‐Band InSAR Time Series

We present deformation measurements of the Kirishima volcanic complex from ALOS and ALOS‐2 Interferometric Synthetic Aperture Radar (InSAR) time series during 2006–2019. Shinmoe‐dake deflated ∼6 cm during the 2008–2010 phreatic eruptions and inflated ∼5 cm prior to the 2017 magmatic eruption. Iwo‐yama inflated ∼19 cm within the crater since January 2015 and ∼7 cm around the southern and western vents since 4 months before the 2018 eruption. These deformations can be modeled as ellipsoids at ∼700 m depth beneath Shinmoe‐dake and as a sphere on top of an ellipsoid at ∼130 and ∼340 m depths beneath Iwo‐yama. Combining geodetic, geoelectric, geochemical, and petrological analysis, we interpret the hydrothermal origin of the deflation at Shinmoe‐dake and inflation at Iwo‐yama; the hydrothermal–magmatic transition during the 2011 Shinmoe‐dake eruption; and water‐boiling and bottom‐up pressurization as driving mechanisms of the inflation at Iwo‐yama. The study highlights the imaging potential of InSAR time series on complex hydrothermal systems.

Here, we report a decade-long Interferometric Synthetic Aperture Radar (InSAR) deformation time series of the Kirishima volcanic complex, Japan, during 2006-2011 and 2014-2019, covering the 2008 and 2017 Shinmoe-dake eruptions and the 2018 Iwo-yama eruption. We estimate the location, geometry, and volume change of the pressure sources and combine them with previous observations from resistivity, geochemistry, and petrology to a new model for the shallow plumbing system of the central Kirishima volcanoes.

Geological Setting
The Kirishima volcanic complex (Japanese for foggy mountain) is located in southern Kyushu, where the Phillippine Sea Plate is subducted beneath the Amurian Plate (Wallace et al., 2009). The complex consists of more than 25 craters, cones, and lava domes produced by the southward migration of eruption centers in the last 330 ka (Figure 1a; Nakada et al., 2013). These volcanic centers form an elliptical 30 by 20 km northwest-trending zone with younger volcanism generally in the southeast (Chapman et al., 2009). The most active eruptive centers are Iwo-yama (Japanese for sulfur peak), Shinmoe-dake, and Ohachi (altitudes of 1,313, 1,421, and 1,408 m, respectively). Hydrothermal systems are widely distributed in shallow levels throughout Kirishima (e.g., Aizawa et al., 2014;Ohba et al., 1997). It has been suggested that the magmatic system beneath connects at depth with Aira caldera located 20 km to the south (Brothelande et al., 2018).
The 2011 Shinmoe-dake eruption is the first magmatic eruption in the complex since its 1716-1717 eruption. Other eruptions at Shinmoe-dake in 1822, 1959 and 1991 are phreatic (Imura & Kobayashi, 1991). Iwoyama, the youngest volcanic center, was formed in the sixteenth to seventeenth century and had a phreatic eruption in 1768 (Tajima et al., 2014). Ohachi had a series of eruptions between 1880 and 1923 (Global Volcanism Program (GVP), 2013).

The 2008-2019 Activity
The unrest of Shinmoe-dake started with a substantial increase in seismicity and a rapid strain change 3 days before the first phreatic eruption on August 22, 2008, when a lake was present in the crater (Geshi et al., 2010;Yamazaki et al., 2020). Additional phreatic eruptions occurred from March to July 2010 (solid blue line/box in Figure 1i). The 2011 eruption started with a phreatomagmatic eruption on January 19 and three sub-Plinian eruptions on January 26-27, followed by stages of lava extrusion, Vulcanian and phreatomagmatic eruptions until September (Nakada et al., 2013). Shinmoe-dake had new magmatic eruptions on October 11-17, 2017, and March 1 to June 27, 2018 (solid orange line/box in Figure 1i; GVP, 2013).
In December 2009, more than a year after the first phreatic eruption, GPS and InSAR data showed inflation over the western flank of the volcanic complex that was attributed to an inflating pressure source beneath Ebino-dake at ∼10 km depth (dashed circle in Figure 1a; Figure 1h; Miyagi et al., 2013;Nakao et al., 2013). Rapid strain change is reported 3 days before this inflation (Yamazaki et al., 2020). The source deflated during the climactic phase of the 2011 Shinmoe-dake eruption and reinflated until November 2011 (Figure 1h; Nakao et al., 2013;Ueda et al., 2013).
About an hour and a half prior to the first sub-Plinian eruption, tiltmeter and broadband seismometer recorded localized inflation near the crater suggesting a shallow pressure source (Takeo et al., 2013). In the following 2 weeks, this source underwent a sequence of inflation-deflation cycles during the sub-Plinian, lava extrusion and Vulcanian stages. The deformation signals, synchronized with volcanic tremors or long-period events, were attributed to the pressurization of a shallow conduit beneath the crater (Nakamichi et al., 2013;Takeo et al., 2013). Localized deflation and inflation patterns were also observed from November 2011 to May 2013 and prior to the 2017 magmatic eruption (Miyagi et al., 2014;. The unrest of Iwo-yama started with an increase in seismicity in December 2013, followed by tremors in August 2014, thermal anomalies and weak fumarolic activity since December 2015, a steam blowout event on April 27, 2017, at the crater, and small phreatic eruptions on April 19 and 20, 2018, at two newly appeared vents on the southern and western side of the crater with plume heights of 100-200 m (Tajima et al., 2020; dashed blue lines in Figures 1g-1i). Localized pre-eruptive inflation suggests a very shallow pressure source at 150 m depth beneath the crater (Narita et al., 2020). Fumarolic activity and mud ejection continued from the southern and western vents as of December 2019 (Japan Meteorological Agency (JMA), 2019).

InSAR Data and Method
We use 2006-2011 ALOS (ascending track 424 and descending track 73) and 2014-2019 ALOS-2 (ascending track 131 and descending track 23) L-band stripmap SAR imagery to form small temporal and spatial baseline interferograms (see Figure S1 and Table S1 for detailed configurations). We remove the topographic phase using the Digital Elevation Model (DEM) released by Geospatial Information Authority of Japan (GSI, 0.4 arc second, ∼10 m) and unwrap the phase using the minimum cost flow method (Chen & Zebker, 2001). Ionospheric delays are not corrected.
We use the stripmap stack processor (Fattahi et al., 2017) of the ISCE software (Rosen et al., 2012) for interferogram processing and the Miami InSAR time series software in Python (MintPy) for time series analysis (Yunjun et al., 2019). We exclude low-coherence interferograms using coherence-based network modification with a custom area of interest around Shinmoe-dake (black empty square in Figure S2) for the average coherence calculation and thresholds of 0.7 for ALOS descending track 73 and 0.8 for the others. We correct for the stratified tropospheric delay using the ERA-5 atmospheric reanalysis model (Hersbach et al., 2020;Jolivet et al., 2011), for long spatial-wavelength phase components by removing linear phase ramps from all acquisitions and for topographic residuals due to the DEM error (Fattahi & Amelung, 2013; Figure S3). We use a temporal coherence threshold of 0.8 to eliminate unreliable pixels. Noisy acquisitions with residual phase root mean squares larger than the predefined cutoff (1 and 2 median absolute deviations for ALOS and ALOS-2, respectively) are excluded during the estimation of topographic residual and average velocity (empty circles in Figure 1g; Figure S4).
To obtain the optimal InSAR measurement for time periods of interest (Figures 1d-1f), we apply two extra steps in addition to the routine MintPy workflow ( Figures S5-S8). First, to maximize the number of reliable pixels, we exclude interferograms with acquisitions after the 2011 and 2017 eruptions, which are decorrelated by local processes inside the crater and/or by the newly deposited ash nearby (Figures S9 and S10). Second, to mitigate residual atmospheric turbulence, we estimate the linear line-of-sight (LOS) velocities for the time periods of interest and convert them to cumulative displacements instead of using the differential displacement between two acquisitions (see Figure S11 for a comparison).
Within the Shinmoe-dake crater, lava domes from the 2011 and 2017 eruptions dominate the topography change (DEM error) since the DEM generation before 2011. We use the locally referenced DEM errors from ALOS-2 ascending and descending InSAR time series between the 2011 and 2017 eruptions to measure the thickness and volume of the lava dome extruded from the 2011 eruption (see Figure S12 for details), providing an alternative approach to the SAR intensity simulation and the single-pass SAR interferometry (Ozawa & Kozono, 2013;Shimono et al., 2011).

Results
We obtain the quasi-vertical displacements from ascending and descending data (Wright et al., 2004) during 2006-2019 to examine the deformation at the Kirishima volcanic complex and find five distinct spatial patterns: three at Shinmoe-dake and two at Iwo-yama during three different time periods (Figures 1d-1f and S13-S15). Shinmoe-dake deflated ∼6 cm between the 2008 and 2010 phreatic eruptions (blue colors in The LOS displacement time series show temporal details (Figure 1g). At Shinmoe-dake, the eastern crater rim (point A) shows no noticeable deformation prior to the 2008 phreatic eruption (solid blue line) and ∼6 cm of linear LOS decrease (subsidence) between the 2008 and 2010 phreatic eruptions; the western crater rim (point B) shows ∼4 cm of LOS increase (uplift) prior to the 2017 magmatic eruption (solid orange line) and a net ∼4 cm of near-linear LOS decrease after the 2018 magmatic eruption. At Iwo-yama, the crater (point C) shows ∼19 cm of near-linear LOS increase during 2014-2019 while the southern and western vents (points S and W) show no significant displacement before December 2017 but ∼5 and ∼7 cm of LOS increases afterward. This horizontal expansion since December 2017 coincides with the increased seismic activity (Figure 1i). Since the end of 2018, Iwo-yama shows no noticeable displacement except for the vicinity of its western vent where linear LOS increase continues until August 2019.
Note that the relatively strong localized inflation on the western summit flank of Shinmoe-dake prior to the 2017 eruption (Figure 1e) is likely related to a potentially partially solidified fissure from which the previous 2008-2010 eruptions occurred (Geshi et al., 2010). We do not interpret the pre-eruptive, co-eruptive and post-eruptive deformation of the 2018 Shinmoe-dake eruption as shown in the GPS baseline change (Figure 1h) due to the overlap with the erupted ash/tephra from the 2017 eruption (Figures 1f, S7, S8, and S16). Similar ash/tephra deposits are also observed during the 2011 eruption from phase ( Figure S17) and coherence measurement (Jung et al., 2016).

Modeling
We use geophysical inverse models to constrain sources of deformation at Shinmoe-dake during 2008-2010 and 2015-2017 and at Iwo-yama during 2015-2017 and 2017-2019. We assume an isotropic elastic half-space and use the vectorized finite compound dislocation model (CDM), composed of three mutually orthogonal rectangular dislocations (with semiaxes , , X Y Z a a a ) with uniform opening and full rotational degrees of freedom (Nikkhoo et al., 2016; see also Beauducel et al., 2020). The CDM represents a generic ellipsoid eliminating the need to prespecify the source geometry such as sphere or sill. We substitute two semiaxes , Y Z a a of the CDM with their dimensionless ratios to X a as / Y X a a and / Z X a a , for easier constraints on the source shape if needed. For Shinmoe-dake during 2015-2017, we fixed the shape and orientation of the CDM due to the lack of near-field observations (we eliminated data points inside the crater affected by local processes). For Iwo-yama during 2017-2019, we use two CDMs and fix one of them with the geometry and location inferred from the previous 2015-2017 period. Although hydrothermal processes deform the ground in a thermo-poro-elastic fashion, simple elastic models are well suited to infer the source geometric features for deformation from distinct episodic unrest where poroelastic response dominates (Fournier & Chardot, 2012;Lu et al., 2002). If mechanically weak layers are present the sources will be deeper than estimated here (Manconi et al., 2007). We use a Poisson's ratio of 0.25.
We account for the elevation effect of the topography using the varying-source depth method (Williams & Wadge, 1998). The height value from GSI DEM is converted from ellipsoid to geoid (31.2 m difference), which is not negligible considering the shallow depth. To fulfill the positive depths requirement of the CDM solution, we apply a minimum free surface height constraint of 1,100 m for Shinmoe-dake and of 1,300 m for Iwo-yama. The impact of this constraint is negligible for all models except for Iwo-yama during 2017-2019, where the lower source of the two could be shallower (<56 m) than estimated here but within the reported 95% confidence intervals (±100 m; see Table S2 for a detailed evaluation).
We jointly invert the ascending and descending InSAR LOS displacement measurements using the Bayesian approach in GBIS (Bagnardi & Hooper, 2018). We subsample the data using a gradient-based adaptive quadtree method (Jónsson et al., 2002) in the near field and use uniform sampling in the far-field (where the signal-to-noise ratio is low; Table S3; Figures S18 and S19). We account for the data uncertainties using structural functions (Lohman & Simons, 2005; Table S3). We use uniform prior probability density functions (PDFs) bounded by geologically realistic values. The inversion algorithm samples posterior PDF of source model parameters with 1,000,000 iterations. The optimal (maximum a posteriori probability) parameter value and 95% confidence intervals are shown in Table S4. All free model parameters converged well except the semiaxis along the X-axis for Shinmoe-dake during 2015-2017 and the opening of the upper CDM for Iwo-yama during 2017-2019. There are trade-offs between some parameters (see joint PDF in Figures S20-S23) but do not affect the depth and derived volume change.
For Shinmoe-dake, the optimal solution for the 2008-2010 deflation is a slightly inclined prolate ellipsoid (Figures 2a-2e) beneath the northeastern crater section with the centroid at depth of ∼620 ± 50 m below the summit (black star in Figures 2c and 2d). For the 2015-2017 inflation, the optimal solution is a sphere beneath the crater center at depth of ∼730 ± 210 m below the summit (black star in Figures 2h and 2i). The estimated cavity volume changes for the two time periods are −124 ± 26 × 10 3 and 146 ± 95 × 10 3 m 3 , respectively.
For Iwo-yama, the optimal solution for the 2015-2017 inflation is a sphere with ∼70 m semiaxes at depth of ∼130 ± 10 m below the summit with an estimated cavity volume increase of 15 ± 2 × 10 3 m 3 . The optimal solution for 2017-2019 inflation is two CDMs on top of each other: one sphere with fixed parameters from the 2015-2017 period except for a free opening bounded by the 95% confidence intervals assuming a constant opening rate and one ellipsoid located at depth of ∼340 ± 100 m below the summit with elongated dimension along the east-west direction. The estimated total cavity volume increase is 76 ± 39 × 10 3 m 3 .

Discussion
The InSAR time series deformation provides new insights into the shallow volcanic system (within the volcanic edifice) of the central Kirishima volcanoes, in addition to a well-established deep source beneath Ebino-dake at ∼10 km depth below the surface.

Shinmoe-Dake: Hydrothermal Origin of the 2008-2010 Deflation
The deflation between the 2008 and 2010 phreatic eruptions is of hydrothermal origin, as most of the erupted material before March 2010 is hydrothermally altered (<1 vol% of juvenile material; Suzuki et al., 2013). This is consistent with geoelectrical-detected widespread shallow low-resistivity, water-saturated porous layers (e.g., Aizawa et al., 2014). Furthermore, most of the deflation occurred before December 2009 when the deep source began to inflate (point A in Figure 1g).
The prolate geometry from geodetic modeling could indicate the depressurization of a former magmatic conduit. The estimated cavity volume decrease (124 ± 26 × 10 3 m 3 ) is similar to the volume of erupted tephra (120 × 10 3 m 3 ; Geshi et al., 2010) from the August 2008 eruption. No noticeable co-eruptive deformation is observed from InSAR (point A in Figure 1g). A possible mechanism is that the August 2008 eruption unsealed the system and emptied the reservoir, which was immediately refilled by hydrothermal fluids; these fluids were released via steam emission (JMA, 2019) during 2008-2010, resulting in surface subsidence.

Shinmoe-Dake: Replacement of Hydrothermal System by Magmatic Body
The 2015-2017 inflation prior to a magmatic eruption was almost certainly magmatic, as indicated by the elevated SO 2 emission higher than the background level since 2011 (JMA, 2019). The inferred depth of 730 m is very similar to 2008-2010 (overlapping PDFs in Figure 3a). A similar source depth was also inferred for the Vulcanian stage of the 2011 magmatic eruption (Takeo et al., 2013) and for November 2011 to May 2013 magma extrusion period (Miyagi et al., 2014). These depths suggest that the deflation/inflation during the above four periods is from the same source. Together, the similar geodetic source depths, the appearance and dominance of SO 2 fluxes and juvenile materials from ash samples suggest that part of the porous rock mass that acted as a hydrothermal source in 2008-2010 was replaced by a magmatic body during the 2011 eruption and stayed as magmatic since then. This is similar to the 2014-2015 phreatic-magmatic transitions in Turrialba volcano inferred from gas emission data (de Moor et al., 2016) but at a much shallower level and transited only once.

Iwo-Yama: Triple-Source Hydrothermal System
The expanded inflation at Iwo-yama since December 2017 cannot be explained by the single source at 130 m depth (very well constrained with a 95% confidence interval of <10 m due to the shallow depth and coherent near-/far-field observation). Geodetic modeling suggests that an extra east-west oriented lay-down cigar-shaped source further down at depth of 340 m (Figures 2p-2t and 3b) with two sources on top of each other at very shallow depths prior to the eruption. Both sources are hydrothermal as evident from the strong fumarolic activity, steam emissions, and ejection of mud and hot water since 2014 (JMA, 2019).
The upper source at 130 m depth (1,180 m a.s.l.) has been inflating since at least January 2015 (point C in Figure 1g). There is a minor peak location change during 2016-2018 around the Iwo-yama crater (Narita et al., 2020). Geochemical analysis in 2016 shows that the water is meteoric (Tajima et al., 2020), suggesting the upward migration of magmatic fluid is not the primary direct driver. Geoelectric surveys have resolved a low-resistivity zone at 200-700 m depth (Tsukamoto et al., 2018), which is below the upper source thus cannot serve as caprocks for hydrothermal sealing to support persistent overpressure. The fumarolic temperature has continuously been around or above 96°C (Tajima et al., 2020), which is the water boiling point temperature at 1,180 m a.s.l. Thus, we conclude that the liquid-gas two-phase boiling of meteoric water is the primary driver of the persistent localized inflation around the Iwo-yama crater.
The lower source at 340 m depth hosts the southern and western vents and is responsible for the April 2018 eruption. Considering (1) the mix of meteoric and magmatic water and the sharp increase of SO 2 and H 2 concentrations in fumarolic gas from the geochemical analysis in 2018 before the eruption and (2) (Ohba et al., 2021;Tajima et al., 2020), we interpret the expanded inflation is driven by the bottom-up pressurization from the upward migration of magmatic fluid. The magmatic fluid is likely generated by the interaction between local meteoric water and magmatic vapor rising from the depth (Ohba et al., 2021).
Residuals of geodetic modeling show a localized uplift at the western vent after December 2017 ( Figure S24 and the difference between solid lines and empty circles in Figure 2t). The small spatial scale (∼150 m in diameter) and the exact same location from ascending and descending orbits suggest that the corresponding pressure source is even shallower than the one at 130 m depth. In total, Iwo-yama hosts a complex hydrothermal system with three tiny reservoirs at very shallow depths.

Conceptual Model of the Kirishima Plumbing System
The magmatic and hydrothermal plumbing system of the Kirishima volcanic complex is summarized in Figure 3c. A deep (∼10 km below the surface) magmatic source beneath Ebino-dake inflated between December 2009 and the 2011 eruption ( Figure S17) and inflated again between the 2017 and 2018 eruption ( Figure 1h). The shallow source beneath Shinmoe-dake (∼700 m below the summit) was hydrothermally evacuating gas and steam between the 2008 and 2010 phreatic eruptions, then turned magmatic during the 2011 eruption and has been accumulating magma since at least January 2015, feeding the 2017 and 2018 magmatic eruptions. The very shallow hydrothermal reservoir beneath Iwo-yama (130 m below the summit) has been boiling since at least January 2015. In December 2017, another slightly deeper hydrothermal reservoir (340 m below the summit) started accumulating magmatic fluid, feeding the April 2018 phreatic eruption. These shallow hydrothermal reservoirs beneath Iwo-yama are possibly connected to the subvertical conductor (2-10 km b.s.l.; Aizawa et al., 2014) via magma degassing and hydrothermal fluid migration (Ohba et al., 2021;Tsukamoto et al., 2018), thus, sharing the same origin of deep magmatic source with Shinmoe-dake.

Conclusions
We derived the surface deformation history of the Kirishima volcanic complex during 2006-2011 and 2014-2019 using InSAR time series analysis. Built upon the routine workflow, we demonstrated that excluding interferograms after destructive eruptions and using average velocity could increase the spatial coverage of near-field observation and further beat down the residual tropospheric turbulence after tropospheric correction, to improve the displacement estimate for time periods of interest. We also described a new approach for lava dome mapping using DEM errors from time series analysis.
Combining geodetic and petrologic data, we identified the hydrothermal reservoir beneath Shinmoe-dake subsiding after the August 2008 eruption, much earlier than the previously thought precursory magmatic inflation since December 2009. The compilation of multiple geodetic studies suggests that this hydrothermal reservoir is replaced by a magmatic body during the 2011 eruption. At Iwo-yama, the deformation time series reveals a precursory horizontal expansion for the 2018 eruption. We propose a hydrothermal system with three tiny reservoirs beneath Iwo-yama at very shallow levels to explain the observed deformation. Combining geodetic, geoelectric, and geochemical data, we conclude that the deformation is driven by liquid-gas two-phase boiling of meteoric water at the upper reservoir and bottom-up pressurization from upward migration of magmatic water at the lower reservoir. The study highlights the imaging potential of InSAR time series for complex hydrothermal systems and the importance of multidiscipline data for understanding volcanic processes.

Acknowledgments
The ALOS and ALOS-2 data are shared among the PALSAR Interferometry Consortium to Study our Evolving Land surface (PIXEL) and provided by JAXA under contract with the University of Tokyo. The ownership of SAR data belongs to JAXA and the Ministry of Economy Trade and Industry. We thank GSI for the DEM and GPS data and JMA for the earthquake data. We thank Jamie Farquharson, Heresh Fattahi, Mehdi Nikkhoo, Piyush Agram, and Bhuvan Varugu for discussions. This work was supported by NASA Headquarters under the Earth and Space Science Fellowship program (grant no. NNX15AN13H) and National Science Foundation's Geophysics program (grant no. EAR1345129).