**4.1. Evaluation of each candidate**

Waveform inversions were performed on the four data combinations resulting in 11 candidate source locations distributed from north to south (Fig. 10). The residual values of the candidates from data combinations (2) and (4) were smaller than those from data combinations (1) and (3), likely because they did not contain the NS component of N.KSHV.

The solution could not be uniquely identified, and multiple solutions were obtained from each data combination. In addition, many candidates did not explain the EW component of N.KSHV. These results may reflect the fact that only three stations were available and that the ground tilt motions dominantly affect ULP signal analyses. Here, we discuss the likelihood of each source location.

First, we evaluate the source candidates (b), (c), (f), and (k). These candidates are located away from the craters, fumaroles, hot springs (Fig. 10), and hypocenters of earthquakes (JMA 2018). They are also separate from a conductive region that has been interpreted as being a hydrothermal fluid reservoir (the C2 conductor identified by Matsunaga et al. (2020)). Therefore, candidates (b), (c), (f), and (k) are considered to be unlikely source locations.

Candidate (h) is located close to fumaroles and hot springs, where hydrothermal activity is more prevalent, and is not separate from the hydrothermal reservoir (Matsunaga et al. 2020). However, this source is slightly separated from the eruptive craters and appears to be inconsistent with the tiltmeter records near Yugama crater (Terada et al. 2021), suggesting that this candidate is unlikely to be the source location.

Candidate (a) is also close to fumaroles and hot springs, and only candidate (a) can explain the waveform of the NS component of N.KSHV (Fig. A1 in the Additional file). However, this candidate has three drawbacks. First, it is slightly separated from the eruptive craters. Second, it is located in a high-resistivity region, which suggests the existence of rock bodies from past volcanic activity (~0.6 Ma) (Hayakawa and Yui 1989; Matsunaga et al. 2020). This region is therefore unlikely to form a path of volcanic fluid. Third, the small residual region around candidate (a) is unnaturally narrow. Therefore, candidate (a) is unlikely to be the source location.

The remaining candidates are (d), (e), (g), (i), and (j). We estimated the sizes of the tensile cracks for these candidates. The moment tensor of a tensile crack is expressed as

$$M\left(t\right)=\left[\begin{array}{ccc}\lambda +2\mu & 0& 0\\ 0& \lambda & 0\\ 0& 0& \lambda \end{array}\right]\varDelta V\left(t\right)$$

,

where \(\lambda\) and \(\mu\) are the Lamé constants and \(\varDelta V\left(t\right)\) is the time history of the volume change of the crack (e.g., Aki and Richards 2002). The volume change for a rectangular tensile crack is expressed as

$$\varDelta V\left(t\right)=LW\varDelta d,$$

where *L*, *W*, and \(\varDelta d\) are the length, width, and aperture change of the crack, respectively. The source time functions shown in Fig. 9 are defined as *S*(t) =\(\lambda \varDelta V\left(t\right)\). We divided the maximum amplitude of the source time function (\({N}_{max}\) = max{*S*(*t*)}) by \(\lambda =\mu =\rho {V}_{p}^{2}/3\) (6.44 × 109 Pa) to investigate the maximum value of the volume change (\(\varDelta V\)). To estimate \(L\), \(W\), and \(\varDelta d\), we need to assume the ratios \(W/L\) and \(L/{\Delta }d\). Taguchi et al. (2018) investigated these ratios for shallow long-period Kusatsu-Shirane sources using a frequency analysis. The estimated crack geometries differed by event: \(W/L\) ranged from 0.45 to 0.83 with typical values around 0.70 and \(L/d\), where *d* is the crack aperture, ranged from 890 to 14,440 with typical values around 2000–7000. The \(L/d\) ratios of feeder dikes measured by field observations at the Miyakejima volcano are 800–5000 with typical values around 3000 (Geshi et al. 2020). Based on these observations, we assumed \(W/L\) = 0.7 and \(L/\varDelta d\) = 3000 to investigate the length, width, and aperture change of the tensile cracks at candidates (d), (e), (g), (i), and (j).

The estimated crack lengths of candidates (i) and (j) were 2200–2300 m (Table 2), which seem too large given source depths of 200–300 m with shallow dip angles. Therefore, candidates (i) and (j) are unlikely.

Candidate (e) is a sub-horizontal tensile crack opening in the N–S direction and located south of the eruptive craters (Fig. 10). In the case of candidate (e), the crack size is estimated to have a volume change (\(\varDelta V\)) of 1.04 × 106 m3, a width (\(W\)) of 1.15 × 103 m, a length (\(L\)) of 1.65 × 103 m, and an aperture change (\(\varDelta d\)) of 5.49 × 10−1 m from the maximum amplitude of the source time function, \({N}_{max}\) = 6.73 × 1015 N m (Table 2). This source is north dipping (Fig. A12), suggesting that it is not easy to connect to the eruptive vents. Therefore, candidate (e) is not viable.

Candidates (d) and (g) are N–S-opening sub-vertical tensile cracks located close to the eruptive vents in Kagamiike-kita crater (Fig. 10). The locations and orientations of these two cracks are similar. In the case of candidate (d), the crack is thought to have dimensions of \(\varDelta V\) = 3.38 × 105 m3, \(W\) = 7.92 × 102 m), \(L\) = 1.13 × 103 m, and \(\varDelta d\) = 3.77 × 10−1 m based on the maximum amplitude of the source time function, \({N}_{max}\) = 2.18 × 1015 N m (Table 2). The crack size of candidate (g) is similar to that of candidate (d) (Table 2). The observed waveforms were well explained by the synthetic waveforms obtained from these candidates (Figs. A4 and A7 in the Additional file). Subsequently, there is no difficulty in considering candidates (d) and (g) as viable sources.

Therefore, a N–S-opening sub-vertical tensile crack near the eruptive craters (Fig. A13 in the Additional file) is the most likely source. The orientation of this crack is consistent with those estimated by Terada et al. (2021) and Yamada et al. (2019). The estimated volume change is also similar to that estimated by Terada et al. (2021) and two orders of magnitudes greater than that estimated by Yamada et al. (2019), suggesting that most of the volume change occurred more slowly than the upper limit period used in the latter study (50 s). This is consistent with the time scale of the source time function (Fig. 9).

**4.2. Details of the possible candidates**

Based on the above discussion, we conducted a second search with finer grid intervals for candidates (d) and (g) (Table 3). The spatial and crack orientation ranges for the second search were determined to cover the grid nodes for which the residuals were less than 1.05 times the smallest residual in the first search (Table 3). The second grid search obtained new source locations around candidates (d) and (g). The results of the second search are summarized in Table 4.

Around the location of candidate (d), there are two candidates labeled (d1) and (d2) (Fig. 11). Candidate (d2) has a smaller residual value than candidates (d1) and (d). However, candidate (d2) has similar drawbacks to candidate (a). Therefore, candidate (d2) is unlikely. Candidate (d1) is located close to candidate (d), and the orientation of its tensile crack is very similar to that of candidate (d). This candidate is steeply dipping to the north with the estimated strike of the crack being consistent with the orientation of the crater row (Fig. 12). The waveform fits and source time function obtained for candidate (d1) are displayed in Figs. 13 and 14, respectively. Comparisons of the waveform fits (Figs. A4 in the Additional file and 13), source time functions (Figs. 9 and 14), and crack sizes (Tables 2 and 4) of candidates (d) and (d1) show that there is no significant difference between the two candidates.

Near the location of candidate (g), there is a candidate labeled (g1) (Fig. 15). Candidate (g1) is located close to candidate (g), and the orientation of its tensile crack is very similar. This candidate is steeply dipping to the north. The observed waveforms were well-fitted by the synthetic waveforms obtained from candidate (g1) (Fig. 16). There are few differences between candidates (g) and (g1) with respect to their waveform fits (Fig. A7 in the Additional file and Fig. 16), source time functions (Figs. 9 and 14), and estimated crack sizes (Tables 2 and 4).

The above implies that the source mechanism of the phreatic eruption is a steeply dipping northward and N–S-opening crack. All of the estimated parameters, including the time difference between the eruption onset and the maximum moment, the maximum moment amplitude, the volume change range, and the crack size, were similar for the first and second searches.

**4.3. Relationship between the source mechanism and the phreatic eruption**

The rapid inflation and slow deflation shown by the source time function (Fig. 14) may be due to fluid migration from depth to the eruptive craters through the estimated crack during the eruption (Fig. 17). The estimated crack is located above the C2 conductive zone that extends beneath Mount Shirane and Mount Motoshirane and is interpreted as being a hydrothermal fluid reservoir (Matsunaga et al. 2020). Terada et al. (2021) suggested that a pre-existing sub-vertical crack along the fluid pathway opened during the eruption. The centroid of this crack is located at 1100 m a.s.l., which is between candidates (d1) (1200 m a.s.l.; Table 4) and (g1) (1000 m a.s.l.; Table 4). The above two studies suggested that hydrothermal fluid or heat is supplied from the depth of the hydrothermal reservoir. The fluid migration from the crack to the main crater requires a certain amount of time; this may correspond to the time difference of 8.9–10.7 s between the eruption onset and the maximum moment for the candidate sources (d1) and (g1) (Table 4). The time of the maximum moment amplitude indicates the time of the start of the volcanic fluid leakage from the crack. To estimate the fluid migration velocity, we assume that volcanic fluids first flow out from the edge of the crack closest to the crater. In the case of candidate (d1), the distance from the nearest crack edge to the main crater is approximately 440 m and the time difference is 8.9 s, giving a velocity of 49 m s−1. In the case of candidate (g1), the distance is 1070 m and the time difference is 10.7 s, giving a velocity of 100 m s−1. Therefore, the volcanic fluid may have migrated at a velocity of 49–100 m s−1 from the edge of the source crack nearest the main crater. These values are comparable to the initial velocity of the ballistic ejecta, which was estimated to be more than 52–82 m s−1 based on the spatial distribution of the ejecta (Ogawa et al. 2019). The phreatic eruption formed several eruptive craters including the main crater, west crater, and south crater. In addition, the subsidence after the eruption detected by Himematsu et al. (2020) suggests that fluid under the subsurface of Mount Motoshirane may have moved considerably. Therefore, the calculated velocities can be interpreted as fluid migration.