## 5.1 Phase 2 source

We previously considered that the location and shape of the pressure source remained almost constant and continued to deflate during Phase 2. Therefore, we first estimated the parameters of this pressure source, including its horizontal position and shape, using the crustal deformation data from Phase 2 prior to estimating the Phase 1 source parameters. We adopted an integrated Markov Chain Monte Carlo (MCMC) and inverse analysis approach as the calculation method. More specifically, the optimal solution determined by the MCMC method was used as the initial value for the nonlinear inversion to improve the efficiency of the calculation. The same MCMC method as in Munekane et al. (2016) was employed in the calculation, whereby the posterior probability distributions of the pressure source parameters were randomly sampled. We applied the varying depth model (Williams and Wadge 1998, 2000) in this analysis and assumed a half-infinite homogeneous elastic solid.

Three sub-Plinian eruptions and lava outflow were confirmed during the 2011 Shinmoe-dake eruption, and a single Mogi model pressure source (Mogi 1958) was estimated to the northwest of the crater at 7–10 km bsl based on the GNSS and tiltmeter crustal deformation observations (Nakao et al. 2013; Ueda et al. 2013). The Mogi model source was estimated at almost the same location for the 2018 eruption (JMA 2018), and this pressure source is inferred to have been a magma reservoir that repeatedly accumulates and supplies magma to Shinmoe-dake crater. Therefore, we used a finite spherical pressure source model (McTigue 1987) because this model assumes that a pressure source isotropically inflates and deflates, as in the Mogi model.

The optimal solution in our calculation is a deflation source at 7.1 km bsl to the northwest of Shinmoe-dake crater during Phase 2 (source A, Fig. 6a). Both the horizontal position and depth are determined in almost the same location as the deflation source that was estimated using GNSS data (JMA 2018), but with a smaller volume change (Table 1). The reason for the difference in volume change is that the analysis period for the GNSS-based volume is 1–10 March, whereas the analysis period for the strain- and tilt-based volume is only 6–8 March.

Table 1

Comparison of the estimated model parameters in this study and JMA (2018).

| Latitude (\(^\circ\)N) | Longitude (\(^\circ\)E) | Depth (km) | dV (× 106 m3) |

Spherical Source (source A) | 31.9331 | 130.8139 | 7.1 | −4.8 |

Mogi model Source (JMA 2018) | 31.9331 | 130.8167 | 7.2 | −7.1 |

**5.2 Time-dependent inversion for the Phase 1 and Phase 2 source parameters using the ensemble Kalman filter**

We tested two models with different numbers of set spherical pressure sources for the Phase 1 source based on the Phase 2 source (source A) estimation results due to the observed |E2/E1| changes during Phase 1.

We first considered the case of two spherical pressure sources, whereby magma migration between the two sources is assumed during Phase 1. However, we found that this model was not appropriate for this study, as explained below. We assumed that source A was one of the two pressure sources, such that its position (horizontal position and depth) was fixed; the volume change of source A and all of the parameters of the other pressure source were then estimated via the MCMC method. These parameters were estimated via the deflation due to two pressure sources, sources A and A'. Source A' was modeled ~ 7 km to the northwest and 3 km shallower than source A; the modeled parameters for these two sources are listed in Table 2. Magma movement between these two pressure sources was unlikely during Phase 1 due to the modeled locations of these two pressure sources relative to Shinmoe-dake crater. Therefore, the two-source model is not suitable for modeling Phase 1 of the 2018 eruption.

Table 2

Estimated model parameters during Phase 1 for the case of two spherical sources.

| Latitude (\(^\circ\)N) | Longitude (\(^\circ\)E) | Depth (km) | dV (× 106 m3) |

source A (sphere) | 31.9331 (fixed) | 130.8139 (fixed) | 7.1 (fixed) | −0.27 |

source A' (sphere) | 31.9713 | 130.7465 | 3.8 | −0.04 |

We then considered the case of a single spherical pressure source with an ascending deflation front during Phase 1. Note that we considered a single stationary spherical pressure source that deflated during Phase 2 (“The Phase 2 source” subsection). Therefore, we modeled a spherical pressure source with a continuously ascending deflation front during Phase1. The results of the Phase 2 pressure source estimation show that the E1 and E2 azimuths at ISA are generally radial and transverse to source A, respectively. If we use the Mogi model for simplicity and assume that E1 and E2 are radial and transverse to the Mogi model, respectively, then the |E2/E1| ratio can be expressed as:

$$\left|\frac{E2}{E1}\right| \approx \left|\frac{{r}^{2}+{h}^{2}}{-2{r}^{2}+{h}^{2}}\right|$$

4

,

where \(r\) is the horizontal distance from ISA to the source and \(h\) is the depth of the source. If we assume that \(r\) is constant over time, then the observed convergence of |E2/E1| to ~ 0.7 from 3–4 during the period from the onset of Phase 1 to the early stage of Phase 2 (Fig. 3b) corresponds to a decrease in depth \(h\). Therefore, if we assume a single pressure source during Phase 1, then the continuous ascension of the deflation front to source A can be used explain the strain observations.

We specifically want to know how the pressure source changed over time from Phase 1 to Phase 2. Here, we introduced a time-dependent inversion analysis (Segall and Matthews 1997) using the ensemble Kalman filter (EnKF) to solve this problem.

We solved two equations, the equation of state and the observation equation, during our time-dependent inversion analysis. The equation of state describes the change in model parameters between the current and previous epochs, and the observation equation links the model parameters to the observed values. We dealt with the nonlinear problem of simultaneously estimating the source location and pressure change by treating the source parameters as an ensemble. Here, the equation of state is written as:

$${{X}}_{{k}|{k}-1}^{{l}}=F\left({{X}}_{{k}-1|{k}-1}^{{l}}\right)+{{\delta }}_{{k}}^{{l}}$$

5

,

where \({X}_{k|k-1}^{l}\) is the state-space vector at time \(k\), \({X}_{k-1|k-1}^{l}\) is the state-space vector at time \(k\)–1, and \({\delta }_{k}^{l}\) is the system noise. The observation equation is written as:

$${{y}}_{{k}}=h\left({{X}}_{{k}}^{{t}{r}{u}{e}}\right)+{\in }_{{k}}$$

6

,

where \({y}_{k}\) is the cumulative crustal deformation vector and \({\in }_{k}\) is the observation error. The Kalman filter is defined from Equations (5) and (6) as:

$${{X}}_{{k}|{k}}^{{l}}={{X}}_{{k}|{k}-1}^{{l}}+\stackrel{-}{ {{K}}_{{K}} }\left({{y}}_{{k}}+{{r}}_{{k}}^{{l}}-h\left({{X}}_{{k}|{k}-1}\right)\right)$$

7

,

where \(\stackrel{-}{ {K}_{K} }\)is the Kalman gain and \({r}_{k}^{l}\) is the pseudo-observation error. Note that the \(l\)-th state space vector \({X}_{k}^{l}\) is expressed as Eq. (8), since the cumulative crustal deformation is used as data:

$${{X}}_{{k}}^{{l}}= \left[{{m}}_{{k}}^{{l}} {{c}}_{{k}-1}^{{l}} \right]$$

8

,

where \({m}_{k}^{l}\) is the model parameter at the current epoch and \({c}_{k-1}^{l}\)is the cumulative crustal deformation up to the previous epoch.

Munekane and Oikawa (2017) have demonstrated the effectiveness of using EnKF in the time-dependent inversion analysis of volcanic crustal deformation problems. Here, the depth and pressure of the spherical source were both varied over time. We first defined a set of pressure source parameters and simulated the source parameters using the same arrangement of actual observation stations that obtained the real data prior to performing the calculations using the real data (Additional file 1: Fig. S2). After calculating the strain and tilt changes from these set pressure source parameters via a forward analysis, we added white noise levels of 1 × 10− 8 and 2 × 10− 8 to the simulated strain and tilt data to obtain the pseudo-observed values. The trend of the pressure source depth over time was consistent with the simulation data, and the estimated volume change was within the margin of error. These results confirmed the effectiveness of the time-dependent inversion using EnKF.

The model setup for the real data was as follows. A 70-h analysis period was defined from the start of Phase 1 to the end of Phase 2 (14:00 on 5 March to 12:00 on 8 March), with each epoch being 1 h in length, resulting in 70 epochs. The total number of data observations was 770, as the tilt data consisted of two-component observations from four stations and the strain data consisted of three-component data from one station. Here, we assumed that the crustal deformation was caused by the movement of a finite spherical pressure source, which was used as the source model. The horizontal position (latitude and longitude) of the spherical pressure source was fixed based on the estimated position during Phase 2, and the depth and pressure change (i.e., volume change) were then treated as the unknown parameters that were estimated at each epoch, resulting in 140 unknown parameters for the entire analysis period. We fixed the horizontal position in the model owing to the limited number of observation data and the fact that the depth provides more important information on the pressure source movement than the horizontal position in this study.

It is necessary to determine the optimal balance between the smoothness and residuals of the unknown parameters to ensure the convergence and reliability of the inversion results. Here, we introduced \({\alpha }_{dep}\) and \({\alpha }_{cap}\) as smoothing parameters for the depth and pressure change, respectively, since these two parameters are free in the inversion. The \({\alpha }_{dep}\) and \({\alpha }_{cap}\) were incorporated into Eq. (6) to obtain:

$$\left[\begin{array}{c}{{y}}_{{k}}\\ 0\\ 0\end{array}\right]=\left[\begin{array}{c}h\\ {\alpha }_{dep}{{H}}_{{d}{e}{p}}\\ {\alpha }_{cap}{{H}}_{{c}{a}{p}}\end{array}\right]\left({{X}}_{{k}}^{{t}{r}{u}{e}}\right)+{\in }_{{k}}$$

9

where \({H}_{dep}\) and \({H}_{cap}\) represent the smoothing matrices that provide the temporal parameter gradients for the depth and pressure change, respectively. We calculated the maximum likelihood based on Segall and Matthews (1997) and obtained the optimal hyperparameter values when the maximum likelihood was minimized:

$$\left(Maximum likelihood\right)= \sum _{k=1}^{S}log\left|{{V}}_{{k}}\right|+Slog\left[\sum _{k=1}^{S}{{\nu }}_{{k}}^{{T}}{{V}}_{{k}}^{-1}{{\nu }}_{{k}}\right]$$

10

,

where \(S\) is the total number of epochs, \({V}_{k}\) is the covariance matrix, and \({\nu }_{k}\) is prediction residual, which is the difference between the observed data at time \({t}_{k}\) and that predicted by the state-space vector conditioned using data up to time \({t}_{k-1}\).