Modeling long-term volcanic deformation at Kusatsu-Shirane and Asama volcanoes, Japan, using the GNSS coordinate time series

Long-term deformation of Kusatsu-Shirane and Asama volcanoes in central Japan were investigated using Global Navigation Satellite System (GNSS) measurements. Large postseismic deformation caused by the 2011 Tohoku earthquake—which obscures the long-term volcanic deformation—was effectively removed by approximating the postseismic and other recent tectonic deformation in terms of quadrature of the geographical eastings/northings. Subsequently, deformation source parameters were estimated by the Markov Chain Monte Carlo (MCMC) method and linear inversion, employing an analytical model that calculates the deformation from an arbitrary oriented prolate/oblate spheroid. The deformation source of Kusatsu-Shirane volcano was found to be a sill-like oblate spheroid located a few kilometers northwest of the Yugama crater at a depth of approximately 4 km\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {km}$$\end{document}, while that of Asama was also estimated to be a sill-like oblate spheroid beneath the western flank of the edifice at a depth of approximately 12 km\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {km}$$\end{document}, along with the previously reported shallow east–west striking dike at a depth of approximately 1 km\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {km}$$\end{document}. It was revealed that (1) volume changes of the Kusatsu-Shirane deformation source and the shallow deformation source of Asama were correlated with the volcanic activities of the corresponding volcanoes, and (2) the Asama deep source has been steadily losing volume, which may indicate that the volcano will experience fewer eruptions in the near future.


Introduction
Monitoring of volcanic deformation is a powerful tool to assess a potential volcanic hazard. Volcanic deformation may be a manifestation of pressure or volume changes in underlying magma chambers (e.g., Ozawa et al. 2004), may sometimes be related to intrusions in the form of diking (e.g., Hotta et al. 2016), or inflation/deflation of hydrothermal sources (e.g., Kobayashi et al. 2018). Many volcanoes exhibit deformation prior to eruption; monitoring the deformation allows us to anticipate the crisis beforehand and prepare for the ongoing hazard (e.g., Williams-Jones et al. 2002).
Kusatsu-Shirane volcano is an active volcano located in central Japan (Figs. 1,2). It is a composite volcano with notable groups of pyroclastic cones named Shirane-san, Ainomine, and Motoshirane-san starting from north to south (Fig. 2a). All historic eruptions up to 2018 were phreatic eruptions in the Shirane-san region with a volcanic explosivity index (VEI) (Newhall and Self 1982) of 1 or 2 (Japan Meteorological Agency 2013).
A summary of recent volcanic activities in the Shirane-san area, as per the Japan Meteorological Agency (Japan Meteorological Agency 2021a), is as follows: in 2014, the number of volcanic earthquakes around the Yugama crater and the area to the south of it increased from March to mid-August, and magnetic changes suggesting thermal demagnetization of rocks under the Yugama crater were observed from May to July. Inflation of a shallow source was observed from May to the end of November 2015 (Terada et al. 2021). In 2018, the number of volcanic earthquakes around the Yugama crater increased in April, and thermal demagnetization of rocks under the Yugama crater was observed from April to July. In addition, grayish water was observed in the Yugama crater from June to July. Inflation of a shallow source was observed from the end of April till the end of August. The number of volcanic earthquakes around the Yugama crater increased again in September. Since then, volcanic tremors, grayish water in the Yugama crater, and inflation of a shallow source have been observed several times.
In 2018, a phreatic eruption occurred at the Kagamiike-kita pyroclastic cone in Motoshirane-san, which  (Argus et al. 2011). Active volcanoes are shown as triangles. Circles represent GNSS sites used in the study. Gray circles are those more than 20 km away from Kusatsu-Shirane and Asama volcanoes that are used to estimate parameters representing tectonic deformation Munekane Earth, Planets and Space (2021)  is adjacent to the Kagami-ike crater (Fig. 2a), with one reported casualty (e.g., Himematsu et al. 2020). Asama volcano is one of the most active volcanoes in central Japan, located approximately 25 km south of the Kusatsu-Shirane volcano (Figs. 1, 2). It is a composite volcano that consists of three volcanoes, Kurofu, Hotokeiwa, and Maekake in order of age (Fig. 2b). Maekake is the youngest which started erupting approximately 10,000 years ago; its current activities are concentrated at the Kamayama crater (Japan Meteorological Agency 2013). It has experienced significant explosive eruptions (VEI > 4) in 1108, 1128, and 1783 (e.g., Japan Meteorological Agency 2013; Yasui and Koyaguchi 2004).
A summary of recent volcanic activities in the Asama volcano, as per the Japan Meteorological Agency (Japan Meteorological Agency 2021b), is as follows: From 2015 till 2018, the volcanic activity of the Asama volcano was elevated, resulting in an increase in the number of volcanic earthquakes, and enhanced fumarolic activities and SO 2 emissions. During the period, Asama volcano experienced two small eruptions in June 2015. In June 2020 , the volcanic activity was elevated again, resulting in the increase in the number of volcanic earthquakes and enhanced amounts of SO 2 .
In recent years, dense networks of geophysical observation systems including GNSS and tiltmeters have been deployed around these volcanoes by the Geospatial Information Authority of Japan (GSI), Japan Meteorological Agency (JMA), and the National Research Institute for Earth Science and Disaster Resilience (NIED) among others. Data from these networks, combined with the application of physical models, suggest that volcanic unrest in these volcanoes is often accompanied by inflation of shallow reservoirs or intrusion of a dike in the shallow part of the volcanic edifice. For example, the unrest of the Yugama crater in the Shirane-san cone in Kusatsu-Shirane volcano in 2014 and 2018 was accompanied by the inflation of a shallow reservoir, possibly in the form of a sill, beneath the Yugama crater at approximately 1.2-1.5 km above sea level. It resulted in deformation that was detected by tiltmeter and repeated GNSS surveys in the proximity of the Yugama crater (Japan Meteorological Agency 2018a;Terada et al. 2021;Tseng et al. 2020). Asama volcano has experienced small eruptions in and 2009(Japan Meteorological Agency 2013. During these periods of eruption, the western part of the volcano exhibited north-south extension, which was interpreted as a manifestation of the intrusion of an east-west striking dike whose top was a few kilometers below ground surface as revealed by continuous GNSS observation (Takeo et al. 2006;Aoki et al. 2013).
The existence of deeper sources, possibly a magma chamber that feeds heat, gas or magma to the shallower chambers, has been suggested by deformation detected by continuous GNSS (e.g., Japan Meteorological Agency 2018b) in the case of the Kusatsu-Shirane volcano. It has also been suggested by deformation detected through leveling survey in the case of the Asama volcano (Murase et al. 2007). However, for either volcano, the long-term behavior of such sources, especially relevance with the changes of the shallower chambers has not been clarified.
One of the difficulties in unraveling the behavior of deep sources, especially in the long term, lies in the treatment of tectonic deformation. Due to the coupling of the subducting Pacific Plate and the overriding North American Plate, this area was subjected to a north-east contraction of approximately 0.2 ppm/year prior to the Tohoku earthquake in 2011 (e.g., Sagiya et al. 2000). Following the earthquake, this area was subjected to a northeast-southwest extension exceeding 0.3 ppm/year in 2013 and decaying with time owing to the postseismic deformation after the Tohoku earthquake (e.g., Geospatial Information Authority of Japan 2014). As an example, the extension of 0.3 ppm/year is equivalent to 3 mm for a 10-km baseline. This value is unacceptably large when one deals with subtle volcanic deformation ( ∼ 1 cm ) detected using a widely spanning GNSS network whose range is typically a few tens of kilometers.
Hence, in this study, I first address the removal of the tectonic deformation from the observed GNSS coordinate time series. Then, with the processed coordinate time series, I determine the possible source locations of Kusatsu-Shirane and Asama volcanoes using the Markov Chain Monte Carlo (MCMC) method (e.g., Munekane et al. 2016), employing an analytical model that calculates the deformation from an arbitrary oriented prolate/ oblate spheroid (Cervelli 2013). Finally, I estimate the strengths of the estimated sources in terms of volume changes using linear inversion, and discuss the relationship between the estimated strengths and volcanic activities of both volcanoes.

Data and methods
I used the GNSS dual-frequency receivers deployed by GSI, JMA, and NIED within a rectangular area spanning 1 • × 1 • with both volcanoes located in its central area (Fig. 1). I processed the GNSS data using Precise Point Positioning (Zumberge et al. 1997) to obtain a daily GNSS coordinate time series. First, the GNSS ephemeris and clocks are estimated using the MADOCA software (Takasu 2013): the PPP analysis is then conducted using the RTKLIB (Takasu 2020) with the estimated ephemeris and clocks. The GNSS coordinate time series cover the period 2013-2020. They were resampled into 30-day bins to mitigate noises (97 epochs in total).

Removal of tectonic deformation
As stated in the introduction, the removal of tectonic deformation is critical in evaluating the long-term behavior of volcanic sources. A common strategy represents tectonic deformation through linear trends, and estimates these trends with the GNSS coordinate time series, when volcanic activities are not observed. However, this strategy may not work in this study as the tectonic deformation cannot be represented by linear trends as they are affected by postseismic deformation related to the Tohoku earthquake (e.g., Tobita 2016). Hence, in this study, I adopted representation of tectonic deformation by a low-order polynomial (e.g., Murakami 2005;Takayama and Yoshida 2007;Nakao et al. 2013).
The order of the polynomial was determined so that the resulting polynomial adequately represented tectonic deformation while avoiding overfitting. Considering the complexity of the tectonic deformation in the studied region, which is dominated by postseismic deformation from the Tohoku earthquake (Geospatial Information Authority of Japan 2014), I selected a second-degree polynomial to represent tectonic deformation. In this case, tectonic deformation was represented by the following equation: where T ij represents the tectonic deformation for the ith component at the jth epoch, and x , y represent a northward and eastward distance from a reference point, respectively. Note that the parameters corresponding to x 2 and y 2 are added to the Murakami's formula to make the equation quadratic for better representing tectonic deformation.
The GNSS coordinate time series of stations that are at least 20 km away from the center of each volcano (gray circles in Fig. 1) are used to estimate the coefficients of Eq. (1) at each epoch. To enhance the stability of the estimating process, the temporal gradient of each coefficient is constrained-that is, the equation is solved as follows: where d and G represent the displacement data and the geometrical factor defined in Eq. (1), and H i (i = E, N, U) are temporal gradient operators for parameters representing the East, North and Up components, respectively. m denotes the model parameters defined in Eq. (1). N, U) are hyper-parameters representing the East, North and Up components, respectively, which will be determined to minimize the Akaike Bayesian Information Criteria (ABIC) (Akaike 1980) described by the following equation: where where m represents estimated parameter values. Then . This procedure is repeated for variables β 2 E , β 2 N , β 2 U to conduct a grid search for minimizing ABIC , and the estimated parameter values for which ABIC is minimum are adopted.

Determination of source positions/geometries by MCMC
The GNSS coordinate time series within 20 km of Kusatsu-Shirane and Asama volcanoes (white circles in Fig. 1) are used to determine source positions/geometries. First, I remove tectonic deformation estimated using Eq. (1) from these coordinate time series. Next, to reduce the computational burden of the MCMC deformation source analysis, I further downsampled the deformation data by keeping the data at epoch #17 (May 2014) and then every fifth epoch after the first, up to epoch #92 (June 2020) (16 epochs in total). Then the MCMC deformation modeling analysis is conducted on selected epochs to determine source positions/geometries and strengths in terms of volume changes. Source positions/ geometries are assumed to be common for all periods. Details of source models employed in the MCMC deformation modeling analysis are descried in the following subsection; the MCMC method employed in this study is described in detail in a previous study (Munekane et al. 2016).
In the MCMC deformation source analysis, I applied the varying-depth model proposed by Williams and Wadge (1998) to include the effects of topography with a homogeneous half-space model, in which a different source depth is assumed at each point for which a solution is desired.

The Kusatsu-Shirane source
To represent the deformation source around Kusatsu-Shirane volcano, I used an analytical model for an arbitrarily oriented spheroid with a point-source approximation where the ratio of the horizontal axis to the vertical axis (axis ratio) is estimated, instead of considering the absolute semi-axis dimensions (Cervelli 2013;Xue et al. 2020). This source is intended to represent a deep source that is supposed to be located in the north or north-western part of the Yugama crater (Japan Meteorological Agency 2018b). Notably, there exists a shallow sill-like source in the proximity of the Yugama crater at approximately 1.2-1.5 km above the sea level (Japan Meteorological Agency 2018a; Terada et al. 2021). I assumed that the deformation caused by this source is negligible at the GNSS stations considered in this study, and have excluded this source in the MCMC estimation. The validity of this assumption is discussed later.

Asama sources
I used a shallow dike to represent the expansion source that was frequently observed at the western flank during times of volcanic unrest (e.g., Takeo et al. 2006;Aoki et al. 2013). I adopted the source parameters from Takeo et al. (2006) to represent the source, and estimated only strengths in the analysis. Furthermore, the analytical model for an arbitrarily oriented spheroid with a pointsource approximation (Cervelli 2013;Xue et al. 2020) was adopted to represent a possible deep source (e.g., Murase et al. 2007).

Linear inversion for source strengths
After determining parameters regarding positions/geometries source parameters for the Kusatsu-Shirane and Asama sources, the strengths parameters, in terms of volume change, were evaluated using linear inversion. All the data were used for this process. The strengths were subject to temporal gradient constraint to stabilize the estimated parameters; that is, I solved the following equations: where d and G denote the deformation data and Green's functions, respectively. H i (i = 0, 1, 2) represents the smoothing matrix which gives the temporal gradient of strengths for the ith source, and β i (i = 0, 1, 2) correspond to a hyperparameter. The indices 0, 1, 2 represent the Kusatsu-Shirane, Asama shallow (dike), and Asama deep sources, respectively. m represents a vector containing m j i , which stands for the strength of the ith source at epoch j.
Note that in the calculation of the Green's functions G , I applied the varying-depth model proposed by Williams and Wadge (1998) to include the effects of topography with a homogeneous half-space model.
Equation (6) is solved similarly to Eq. (2). Namely, Eq. (6) is first solved with the least-square method for fixed β 2 i : where m is the solution for Eq. (6). The ABIC for the corresponding solution is given by replacing m with m in the following equation: where This procedure is repeated for variables β 2 0 , β 2 1 , β 2 2 to conduct a grid search for minimizing ABIC and the estimated parameter values for which ABIC is minimum are adopted. Figure 3 shows the selected GNSS coordinate time series at non-volcanic and volcanic sites (whole GNSS coordinate time series are shown in Additional file 1: Figure  S1a, b. It can be seen that the coordinate time series agree well with estimated tectonic deformation for nonvolcanic sites. At volcanic sites, some deviations between GNSS coordinate time series and tectonic deformation can be observed. For example, some deviations are observed during 2013-2015 and after 2018 in the northsouth components at the site kshv, which are of volcanic origin as explained later. To evaluate the goodness of the fit, I estimated the reduced χ 2 , which is defined as where n and n deg are the number of the observed deformation data and the degree of freedom, respectively. x obs

Removal of tectonic deformation
(10) represents the ith observed deformation, and x cal i represents the ith model-calculated deformation, respectively. σ i stands for error for x obs i . The reduced χ 2 is calculated to be 2.76 for non-volcanic sites.
The optimum values of β i (i = E, N, U) are found to be (β 2 E , β 2 N , β 2 U ) = (10., 30., 30.) . The cross-sections of ABIC for each axis at the optimal point are given in Additional file 1: Figure S1c. Table 1 shows the estimated parameters for the source near the Kusatsu-Shirane volcano. It is located a few kilometers west of the Yugama crater at a depth of approximately 4 km (Fig. 4a). It is to be noted that the source is a highly oblate spheroid (axis ratio is over 200) with a dip angle close to 90 • , which may be approximated as the penny-shaped crack (Fialko et al. 2001). Figure 4a shows the location of the source, and  Table 2 shows the estimated parameters for the sources under the Asama volcano. In Table 2, source parameters representing a shallow expansion source of the volcano (Takeo et al. 2006), fixed during the inversion, are also shown. A deep source is found under the western flank of the Asama volcano at approximately 12 km depth, which is also a highly oblate spheroid (axis ratio is over 100) with dip angle close to 90 • , that may be approximated to a penny-shaped crack (Fialko et al. 2001). Figure 5a shows the location of the sources and Fig. 5b shows the observed and calculated deformation between May 2014 and June 2020.

Source determination by MCMC
The reduced χ 2 (Eq. 10) is 2.30. The corner plot representing the posterior PDF of the estimated parameters is given in Additional file 2: Figure S2. Note that only parameters related to the source positions/geometries are plotted for brevity. Figure 6 shows the estimated strengths in terms of volume changes for sources under Kusatsu-Shirane and Asama volcanoes.   The deep source of the Asama volcano shows a general deflation with a plateau from the beginning of 2015 to the end of 2017. It lost approximately 20 × 10 6 m 3 in volume during the analysis period. Figure 7 shows the observed and calculated coordinate time series for selected sites (whole GNSS coordinate time series are given in Additional file 3: Figure  S3a, b. Reasonable fits between observed and calculated coordinate time series are observed. The reduced χ 2 (Eq. 10) is calculated to be 2.61.

Discussion and conclusions
The Kusatsu-Shirane source I suggest this correlation implies that the Kusatsu-Shirane source acts as a heat-source for the shallow source that triggers volcanic activity. In fact, previous works (Matsunaga et al. 2020a, b) presented the threedimensional electrical resistivity structure under the Kusatsu-Shirane volcano. They found a widely spreading conductor (C2a) at a depth of approximately 1-1.5 km observed over an area stretching from the Yugama crater to Motoshirane-san, which has been interpreted as a large-scale hydrothermal fluid reservoir. Further, Matsunaga et al. (2020b) revealed the deep extension of the C2a conductor at least up to the depth of (7.5 km ) (C2b), which they interpreted as a possible magma chamber that drives hydrothermal fluid migration in the C2a conductor. Considering this structural information, the deep source revealed in this study may be interpreted as a magma chamber that corresponds to the C2b conductor in Matsunaga et al. (2020b). This source may provide additional heat to the shallow hydrothermal fluid reservoir during its inflation period that corresponds to the C2a conductor (Matsunaga et al. 2020a, b) to trigger the  Earth, Planets and Space (2021) 73:192 volcanic activities at the Yugama crater and possibly at the Motoshirane-san.

Asama sources
In Fig. 6b, it can be observed that the volume change of the Asama shallow source shows an oscillatory pattern and inflation occurred from the beginning of 2015 through the beginning of 2018 with a short period of deflation in late 2016. Inflation was observed again in mid-2020. Correspondingly, I observed the enhanced volcanic activities in the Asama volcano, as summarized in the introduction. This correspondence supports the view that the shallow source, a nearly east-west-trending dike, acts as a magma deposit, triggering various volcanic anomalies such as increase in volcanic earthquakes (Aoki et al. 2013). Aoki et al. (2013) reported the S-wave velocity structure in and around the Asama volcano and found the widely spread low-velocity zone under the western flank of the volcano at a depth of 5-10 km , which was interpreted as a magma chamber that feeds magma to the shallow dike. The deep source identified in my study was located at the bottom of the low-velocity zone. Hence, I may consider the source as a part of the magma chamber from which the magma migrates upward without causing considerable deformation.
Furthermore, Murase et al. (2007) found a deep source located under the western flank of the Asama volcano at a depth of 6-8 km from the leveling data between 1902-2005. The horizontal location of this source is within the margin of error of the deep source in this study. However, there is a notable difference in the geometries of the two  sources (sphere vs. spheroid) and the depth (6-8 km vs. 12 km ). It requires further investigation to determine if these sources may be considered the same. It is noted that in Murase et al. (2007), the authors did not take into account the effect of the topography, which could be a possible cause of the discrepancy in the estimated depth.
In case these two refer to the same source, the volume changes of the source as revealed in Fig. 6b may have an important role in volcanic monitoring since Murase et al. (2007) claim that the temporal changes in the pressure or volume of the source exhibit a strong positive correlation with the eruption frequency, implying that Asama volcano might experience reduced eruptive activities in the near future because the volume of the deep source rapidly decreased during the study period (2013-2020), as experienced in the period 1943-1967 when fewer volcanic eruptions were observed. For disaster mitigation, it may be important to monitor the volume changes of the deep source using GNSS to estimate inflation, which may be an indicator of the rekindling of volcanic activities.
Then what physical mechanism is causing the deep source to reduce in volume with time? One of the candidates is degassing. Kazahaya et al. (2015) investigated the magma budget of the shallow magma plumbing system of Asama volcano, and pointed out that to explain the large volume of degassed magma ( 4 × 10 8 to 9 × 10 8 m 3 from 1975 to 2011), a significant amount of magma was likely supplied to the shallow dike-like conduit from the mid-crustal magma reservoir while the degassed magma descended back to the mid-crustal magma reservoir. The deep source in my study likely corresponds to the supposed mid-crustal magma reservoir in Kazahaya et al. (2015) and it loses its volume through this magma exchange with the shallow source.

Possibility of close relationships between magma systems
Given the proximity of the two volcanoes, it would be interesting to discuss the possibility that both magma systems are fed by the same magma source and interact. For example, Brothelande et al. (2018) used the GPS time series and deformation modeling to show that Aira caldera and Kirishima, two adjacent volcanic centers in south Kyushu, Japan, interacted during Kirishima unrest in 2011. Namely, whereas the magma source under the Aira caldera had been inflating steadily, it deflated during the eruption of Kirishima in 2011. They interpreted the deflation as the result of magma withdrawal from the Aira source through a common deep reservoir. It is interesting to see if this type of inter-connectivity can be seen in other mutually adjacent volcanic systems. In Fig. 6a, b, I showed the volume changes of the Kusatsu-Shirane source and the Asama deep source, respectively. If these sources are interconnected through a common deep reservoir, it is expected that both changes would positively correlate. However, from Fig. 6a, b, one can see that the Kusatsu-Shirane source is generally inflating while the Asama deep source is generally deflating. Hence, it may be conjectured that both sources are not currently interconnected though it is not sufficient to exclude the existence of a common deep reservoir. In a previous study (Brothelande et al. 2018), the authors note that the volcanic systems are not connected all the time because magma pathways open and close. Hence, at this point, it may be equally probable that the magma pathway from Kusatsu and Asama is temporarily closed or it does not exist at all.

Influence of 2018 Motoshirane-san eruption and inflation of a shallow source in the vicinity of the Yugama crater
Here, I evaluate the magnitude of the deformation from deformation sources omitted in this study, to see their impacts on our model estimations.
The first is the deformation caused by the 2018 eruption of the Motoshirane-san. Himematsu et al. (2020) reported that co-and post-eruptive deformation which is localized around the Kagami-ike crater. Here, I evaluate co-and post-eruptive deformation expected at our GNSS observation sites based on the deformation model presented in Himematsu et al. (2020), whose parameters are summarized in Additional file 4: Table S1. As a result, I found that the maximum co-and post-eruptive horizontal deformation was as small as 2 mm at j423 and the maximum vertical deformation is 0.4 mm at j423. Since these magnitudes are well within the error bars of the deformation observed at respective sites (Additional file 3: Figure S3a), I may conclude that the deformation caused by the 2018 Motoshiranesan eruption is safely ignored in our model estimation.
The second is the deformation caused by inflation of a shallow source in the vicinity of the Yugama crater in the Shirane-san area (Fig. 2a). When volcanic activities at the Yugama crater increased in 2014 and 2018, deformation caused by inflation of the source were detected by tiltmeters and repeated GNSS surveys in the proximity of the Yugama crater (e.g., Japan Meteorological Agency 2020).
Regarding the 2014 unrest, Terada et al. (2021) presented the physical model based on the changes in tiltmeters, whose parameters are summarized in Additional file 5: Table S2. I evaluated the expected deformation on our GNSS sites based on the proposed model, and found that the maximum horizontal and vertical displacement on out GNSS sites caused by the shallow source during the 2014 unrest are as small as 1.1 mm at j423 and 0.4 mm at j424, respectively. Concerning the 2018 unrest, though no model has been proposed, it is reported that the magnitude was much smaller (more than a factor of two) than the 2014 unrest in terms of observed tilt changes (e.g., Japan Meteorological Agency 2020). Hence, I may conclude that the deformation at out GNSS sites caused by the shallow source during the analysis period will not exceed 2 mm in the horizontal components and 1 mm in the vertical component, which is well within the error bars of the deformation at the relevant sites and hence may be safely ignored.

Effect of viscoelastic rheology
In this study, I assumed that the deformation sources were spheroidal cavities embedded in a homogeneous elastic half-space. However, previous studies have shown that viscoelastic rheology may affect the amount and time history of observed deformation (e.g., Newman et al. 2001;Yamasaki et al. 2018). Newman et al. (2001) employed a viscoelastic shell model with a spherical pressure source surrounded by a concentric, spherical, Maxwell viscoelastic shell in an elastic full/half-space (e.g, Dragoni and Magnanensi 1989). In contrast, Yamasaki et al. (2018) investigated an oblate-spheroidal sill embedded in a two-layered medium that consisted of an elastic layer underlain by a viscoelastic layer. Both models have common features as follows: (1) uplift caused by instantaneous inflation will be followed by subsidence even if the source strength remains constant, which may lead to overestimation of volume decreases in the pure elastic model, and (2) during gradual inflation, the magnitude of inflation-induced uplift will be reduced by the relaxation, which may lead to underestimation of volume increases in the pure elastic model. Hence, a volume decrease following a volume increase, as observed during 2016-2017 in the Kusatsu source (Fig. 6a), may be overestimated. Additionally, the total volume increase, as observed during 2014 in the Kusatsu source (Fig. 6a), may be underestimated, because of the viscoelastic relaxation. Hence, detailed quantitative interpretation of the volume changes of a source needs to include the effect of viscoelastic rheology.