Coeruptive and posteruptive crustal deformation associated with the 2018 Kusatsu-Shirane phreatic eruption based on PALSAR-2 time-series analysis

Coeruptive deformation helps to interpret physical processes associated with volcanic eruptions. Because phreatic eruptions cause small, localized coeruptive deformation, we sometimes fail to identify plausible deformation signals. Satellite synthetic aperture radar (SAR) data allow us to identify extensive deformation elds with high spatial resolutions. Herein, we report coeruptive crustal deformation associated with the 2018 Kusatsu-Shirane phreatic eruption detected by time series analyses of L-band satellite SAR (ALOS-2/PALSAR-2) data. Cumulative deformation maps derived from SAR time series analyses show that subsidence and eastward displacement dominate the southwestern side of an eruptive crater with a spatial extent of approximately 2 km in diameter. Although we were unable to identify any signicant deformation signals before the 2018 eruption, posteruptive deformation on the southwestern side of the crater has been ongoing until the end of 2019. This prolonged deformation implies the progression of posteruptive physical processes within a conned hydrothermal system, such as volcanic uid discharge, similar to the processes observed during the 2014 Ontake eruption. Although accumulated snow and dense vegetation hinder the detection of deformation signals on Kusatsu-Shirane volcano using conventional InSAR data, L-band SAR with various temporal baselines allowed us to successfully extract both coeruptive and posteruptive deformation signals. The extracted cumulative deformation is well explained by a combination of normal faulting with a left-lateral slip component along a southwest-dipping fault plane and an isotropic deation. Based on the geological background in which the shallow hydrothermal system develops across Kusatsu-Shirane volcano, the inferred dislocation plane can be considered as a degassing pathway from the shallow hydrothermal system to the surface due to the phreatic eruption. We reconrmed that SAR data is a robust tool for detecting coeruptive and posteruptive deformations, which are helpful for understanding shallow physical processes associated with phreatic eruptions at active volcanoes.


Introduction
A phreatic eruption is a hazardous volcanic activity event that includes sporadic ejections of volcanic ash, steam and volcanic gases induced by transient pressure changes in the hydrothermal systems of active volcanoes. These episodic pressure changes are generally induced by injections of volcanic gases into a con ned shallow hydrothermal system or by the abrupt boiling of con ned overheated geothermal water due to sudden depressurization (Jolly et al. 2014;Narita and Murakami 2018). In general, phreatic eruptions cause less damage than magmatic eruptions, but can occasionally produce destructive hazards such as lahars that wash out infrastructure (e.g., Naranjo et al. 1986). An abrupt eruption can be a potential risk to tourists who are visiting volcanoes for mountaineering in the summer and skiing in the winter.
Ground deformation of active volcanoes is usually helpful for interpreting the physical processes of volcanic activity; however, most of ground deformation associated with phreatic eruptions is characterized as small and localized ones with a sudden onset. These characteristics hinder the detection of plausible deformation signals prior to phreatic eruptions due to sparse observation networks or limited measurement accuracy. Recently, satellite synthetic aperture radar (SAR) data have allowed us to identify displacement elds with high spatial resolution without ground-based instruments. Several studies have successfully detected coeruptive deformation associated with phreatic eruptions using satellite SAR data (e.g., Hamling 2017; Doke et al. 2018; Narita and Murakami 2018). If precursors of phreatic eruptions, such as overpressure in shallow hydrothermal systems, persist for a long time, satellite SAR data can identify the precursory deformation signals of a phreatic eruption (Kobayashi et al. 2018).
Kusatsu-Shirane volcano is an active volcanic complex located in central Japan (Figure 1 Institute 2018). However, some InSAR data were contaminated by decorrelation noises due to variations in the back-scatter characteristics on the ground, such as snow and volcanic tephra coverage. Here, we investigate both the coeruptive and posteruptive deformation signals associated with the 2018 Kusatsu-Shirane phreatic eruption by applying SAR time series analyses to L-band SAR data. We also propose a schematic model that explains the extracted coeruptive deformation. SAR data processing Page 4/20 In this study, we employed L-band SAR data acquired from the Phased Array type L-band SAR 2 (PALSAR-2) sensor onboard the Advanced Land Observation Satellite 2 (ALOS-2) to detect crustal deformation signals on Kusatsu-Shirane volcano (Additional le 1; Table S1). The L-band microwaves (wavelength: 23.6 cm) that PALSAR-2 uses are suitable for monitoring deformation signals on volcanoes covered by dense vegetation, such as those in Japan. In general, decorrelation problems can be caused by variations in scattering characteristics on the ground such as snow coverage or dense vegetation. Shorterwavelength microwaves, such as C-band (wavelength: 5.6 cm) and X-band (wavelength: 3.1 cm) microwaves, scatter on shallower parts of the snow/ice layer and on leaves or branches. Thus, we expect shorter-wavelength microwaves to suffer from decorrelation problems in this case. In contrast, longerwavelength microwaves, such as the L-band microwave, can penetrate deeper parts of a snow/ice layer and through dense vegetation. The region around Kusatsu-Shirane volcano is covered by dense vegetation in the summer and by snow in the winter (Additional le 1; Figure S1).Although the recurrent period of the ALOS-2 satellite (14 days) is greater than that of recent SAR constellations, such as Sentinel-1 (6 days) and Cosmo-SkyMed (1 day at a minimum), coherence using PALSAR-2 data tends to sustain, even with a temporal baseline for the PALSAR-2 InSAR greater than one year. Thus, we expect that the PALSAR-2 data are suitable for extracting deformation signals around Kusatsu-Shirane volcano through the year.
All InSAR data were generated using the GAMMA software (Wegmüller and Werner 1997). We corrected the topography-dependent fringes using a 10 m mesh digital elevation model released by GSI.
Tropospheric artifacts were corrected using zenith tropospheric delays provided by the Generic Atmospheric Correction Online Service for InSAR (GACOS; Yu et al. 2017;Yu et al. 2018). Longwavelength signals across the InSAR data were corrected by tting 2D polynomial functions. We discarded some PALSAR-2 data that were contaminated by strong ionospheric artifacts.
We employed a multi-temporal InSAR (MTI) analysis, one of the SAR time-series analyses, to infer spatiotemporal variations in the crustal deformation at Kusatsu-Shirane volcano (e.g., Schmidt and Bürgmann 2003). The MTI analysis infers mean displacement rates during each image acquisition interval using the InSAR data with various temporal baselines, assuming constant displacement rates during each image acquisition interval. Figure S2 in Additional le 1 shows a plot of the perpendicular baselines and SAR data combinations for the MTI analysis. We did not set the criteria for spatial and temporal baselines in estimating the displacement time-series, unlike the small baseline subset approach (Berardino et al. 2002). One reason for this is that L-band SAR data tend to avoid decorrelation problems even when using a pair of SAR images with a temporal baseline of more than one year. Another reason is that the ALOS-2 satellite has been operating within 500 m of the perpendicular baseline since the satellite was launched. The Laplacian operators for the smoothing temporal variations in line-of-sight (LOS) changes were optimized using the L-curve criterion (e.g., Hansen 1992; Additional le 1; Figure S3). We did not infer any temporal variations in the LOS changes for discarded pixels where the coherence of any individual InSAR data point was below 0.1. After estimating the time-series of the LOS changes using the MTI analysis, we extracted the cumulative deformation associated with the 2018 phreatic eruption until the end of 2019. Using pairs of cumulative coeruptive LOS changes in paths 19/125 and 19/126, we decomposed them into quasi-east-west (QEW) and quasi-up-down (QUD) components to better understand the spatial characteristics of the coeruptive LOS changes (Fujiwara et al. 2000).

Results
Coeruptive deformation associated with the 2018 phreatic eruption Figure 2 shows the cumulative LOS changes in paths 19, 125, and 126 for approximately two years since the 2018 phreatic eruption. In the path 19 MTI results, a positive LOS change was dominant on the southwestern side of the 2018 crater. The maximum amplitude of the positive LOS change was ~6 cm near the 2018 crater. The standard deviation of the inferred displacement velocity in the undeformed region in path 19 was ~1 cm (Additional le 1; Figure S4a). In the path 19 averaged InSAR data, we identi ed a displacement discontinuity with a WNW-ESE strike at the 2018 crater where the path 19 MTI data show missing data (Additional le 1; Figure S5). The four averaged individual interferograms in path 19 showed the same location for the displacement discontinuity (red arrows in Figure S5a in Additional le1). Thus, this displacement discontinuity is likely a plausible characteristic of the coeruptive deformation (Additional le 1; Table S2). The location of the displacement discontinuity is identical to that of the 2018 crater. Unlike the MTI data, the averaged data showed ~5 cm of negative LOS changes on the northeastern side of the 2018 crater where the path 19 MTI data showed missing data. The Next, we decomposed the data into coeruptive quasi-east-west (QEW) and quasi-up-down (QUD) components using the two pairs of cumulative LOS changes in path 19 from the ascending orbit and 125 and 126 from the descending track ( Figure 3). We also plotted the NE-SW cross-sections of both QEW and QUD components along pro le P-P', which crosses the 2018 crater in a NE-SW direction (Figures 3e  and 3f). In the QEW components, nearly the entire region of the displacement elds exhibited eastward movement (Figures 3a and 3c). The maximum amplitude of the eastward movement was ~4 cm near the 2018 crater. Eastward displacements were identi ed on not only the southwestern side but also the northeast side of the 2018 crater. The two QUD displacement maps show predominantly subsidence on the southwestern side of the 2018 crater and less than 2 cm of uplift on the northeastern side of the 2018 crater (Figures 3b and 3d). The cross-sections of the coeruptive QUD displacements exhibit an asymmetric pattern with subsidence of up to ~10 cm near the 2018 crater (Figure 3f). This asymmetric displacement pattern can be explained by a shear dislocation along a plane rather than a point source de ation, which would instead produce a concentric displacement eld in the vertical component. Following the 2018 phreatic eruption, the path 19 MTI results show 2 cm of uplift with a peak in mid-2018. However, we were unable to identify similar characteristics for the LOS changes in paths 125 and 126 during this period. While we found several differences between the temporal characteristics of the MTI results and those of the tiltmeter observations, the deformation signals in 2014-2019 at Yugama crater lake in all paths appear to have similar spatial extents, which is consistent with a non-vegetated region (Additional le 1; Figures S7-S9). We presume that the positive LOS changes at Yugama crater lake during 2016-2017 were derived from the residuals of the post InSAR process, or are the result of seasonal variations in soil moisture due to the sparse vegetation around the lake (Gabriel et al. 1989;Zwieback et al. 2015). We have not performed further analysis of the deformation at Yugama crater lake, but future analyses will improve our understanding of the shallow hydrothermal system across Kusatsu-Shirane volcano.

Discussion
Data tting of the cumulative LOS changes for the 2018 phreatic eruption The PALSAR-2 MTI results allowed us to extract both the deformation signals associated with the 2018 Kusatsu-Shirane phreatic eruption and the time series of LOS changes from 2014-2019 (Figures 3 and  4). We expect that the spatial characteristics of the coeruptive deformation can be explained by a dislocation along a plane, rather than by volume changes of a point source, because of the observed asymmetry of the coeruptive displacement eld. We used the analytical solutions of surface deformation caused by planar dislocation (Model A), volume changes of a point source (Model B), or both (Model C) to t the observed data (Mogi 1958;Okada 1985). We adopted a grid search algorithm to infer the following best-t parameters: width, length, dip angle, rake angle, the amount of slip for the planar dislocation and the location, depth, and volume change for the point source. Table S3 in Additional le 1 lists the ranges and intervals of the grid search parameters used in this analysis. We set the top location of the dislocation plane to that of the 2018 eruptive crater, xed the strike angle at 105 degrees, and set the shallowest depth at 50 m. We regarded the best-t solution as a combination of parameters that minimized the total root-mean-square (RMS) of the residuals between the observed and the computed deformations. We assigned a Poisson's ratio of 0.25. We resampled the LOS change maps in each path using concentric grids (e.g., Fukushima et (e.g., Fukushima et al. 2005). The 95 % con dence intervals of the best-t parameters were estimated by a bootstrap approach with 300 iterative nonparametric resamplings (Efron 1979).
Model C, a combination of a plane dislocation and an isotropic de ation, best ts the observed cumulative deformation associated with the 2018 eruption. It minimizes the total RMS residuals of the observed and computed deformations and exhibits a smaller AIC value than other models (Akaike 1974) (Table 1; Figure S10 in Additional le 1). The best-t geometry of a plane dislocation shows a width of 1500 m, a length of 500 m and a dip angle of 56º, i.e., the plane extends to the dip direction rather than the horizontal direction and dips toward the southwest. The bottom of the dislocation plane reaches at a depth of greater than 1200 m from the surface. The best-t model exhibits a rake angle of -64 degrees and a slip of 0.17 m on the southwest-dipping plane, which indicates normal faulting with a left-lateral component ( Figure S11a in Additional le 1). The total geodetic moment release due to a plane dislocation is 1.3×10 15 Nm, which corresponds to the moment magnitude of 4.0, assuming a shear modulus of 10 GPa (Kanamori 1977). In contrast, the location of the isotropic de ation source was inferred to be located near the eastern edge of the dislocation plane at a depth of 550 m (Figure S13b in Additional le 1). The best-t volume change (-30000 m 3 ) of the point source induces approximately 2 cm of LOS extension at a maximum and up to 0.5 cm of LOS extension within a 1.5 km radius. Thus, we can expect that the near-crater deformations were mainly driven by plane dislocation while wide-range displacements were caused by isotropic de ation. Other single source models (Model A and B) showed greater RMS residuals than those in Model C, while those models indicated similar features for the deformation source parameters as those in Model C (Table 1; Table S4, Figures S12 and S13 in Additional le 1).

Data and model interpretations
The most plausible interpretation of a plain dislocation dipping from the volcanic crater is a path formation of a sudden steam plume and/or a volcanic gas ejection from the shallow hydrothermal system to the surface. The geometry of the best-t plane dislocation can be considered as a degassing pathway (Additional le 1; Figure S14). The best-t width of the plane in Model C is 1500 m, which indicates that the bottom of the plane reaches an approximately depth of 1200 m. In contrast, the best-t point source was inferred to be emplaced at the middle part of the eastern edge of the dislocation plane ( Figure 12). Previous magnetotelluric surveys proposed a low-resistivity subsurface structure from Moto-Shirane volcano to Yugama crater lake at 1500-3000 m below the surface, with alternating thin laminae of low-and high-resistivity layers between the surface and large conductor (Nurhasan et  Considering that the observed asymmetric-shaped deformation can be interpreted as trapdoor faulting, we can presume the abrupt pressure changes occurred within the aquifer of the peripheral hydrothermal system before and after the phreatic eruption. Several historical WNW-ESE aligned craters detected by the GSI airborne SAR images are distributed perpendicular to the distribution axis of the pyroclastic cones on Kusatsu-Shirane volcano (GSI 2018a). Although phreatic eruptions are the dominant eruption type, a geological survey suggested that the sequence of pyroclastic cones formed due to magmatic eruptions ~3000 years ago (Hayakawa and Yui 1989). The orthogonal distribution of small craters and pyroclastic cones may suggest a rotation of the stress regime from a depth at which a magma body is emplaced, to a shallower depth where the hydrothermal system is developed. In addition, the distributions of small craters imply that the physical processes of historical phreatic eruptions were similar to those of the 2018 eruption. The left-lateral component of our best-t solution can be considered as a stress accommodation in the shallow brittle crust stimulated by breaking seals or weaknesses from normal stress due to episodic degassing. Considering that paths of ejected volcanic uid usually form along pre-existing fractures or low-energy pathways, the distribution axis of the aligned WNW-ESE craters implies a peripheral stress regime in the shallow part of the crust above the magma source.

Conclusions
We successfully extracted the coeruptive and posteruptive deformation signals associated with the 2018 Kusatsu-Shirane phreatic eruption based on ALOS-2/PALSAR-2 MTI data. The PALSAR-2 MTI results show that the coeruptive deformation elds are dominantly characterized by ~10 cm of subsidence and eastward movement on the southwestern side of the 2018 crater, which formed on the northern side of the Kagamiike-kita pyroclastic cone. We identi ed approximately 2-3 cm of posteruptive deformation that lasted for approximately two years following the 2018 eruption, while few plausible deformations were observed before the 2018 eruption. A combination of normal faulting with a left-lateral slip component on a southwest-dipping plane and an isotropic de ation yield a favorable solution that ts the spatial characteristics of the extracted cumulative deformation associated with the 2018 phreatic eruption. The inferred best-t dislocation plane can be interpreted as a pathway for volcanic uid transport from the deformation source to the surface. This interpretation is compatible with the geological background in which the shallow hydrothermal system develops across Kusatsu-Shirane volcano. To our knowledge, this is the rst report on crustal deformation associated with the 2018 Kusatsu-Shirane eruption using SAR data.
We recon rmed that satellite SAR data allow us to detect high-spatial-resolution crustal deformation on active volcanoes even if the eruption occurs at an unexpected site or in a region with a sparse groundbased network. Although the coeruptive deformation due to the 2018 Kusatsu-Shirane phreatic eruption was small in magnitude (~10 cm) and in spatial extent (~2 km in diameter), similar to deformation caused by other previous phreatic eruptions, the observation data contribute to an understanding of the shallow physical processes related to the 2018 phreatic eruption. Although we only present SAR image processing associated with the 2018 Kusatsu-Shirane eruption in this study, these data will be helpful for supporting other observational dataset.

Declarations
Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.

Competing interests
The authors declare that they have no competing interests.

Funding
This study is funded by "Integrated program for next generation volcano research and human resource development" led by the Ministry of Education, Culture, Sport, Science and Technology, Japan (MEXT).

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.