Waveform inversion of the ultra-long-period seismic event associated with ground tilt motion during an eruption of Mount Kusatsu–Shirane, Japan, on January 23, 2018

We conducted waveform inversions of an ultra-long-period (~ 240-s) event associated with the phreatic eruption of Mount Kusatsu–Shirane on January 23, 2018. We used broadband seismic and tilt records from three stations surrounding the eruption site. The horizontal components of the broadband seismic records were severely contaminated by tilt motions. We applied a waveform inversion algorithm to account for both the translational and tilt motions. To reduce the number of free parameters, we assumed a tensile crack source and conducted grid searches for the centroid location and orientation of the crack. The results showed a rapid inflation of 105 m3 of the crack, followed by a slow deflation starting 8–11 s prior to the onset of the eruption. The source location and crack orientation were not uniquely determined. The most likely source is a north–south-opening sub-vertical crack near the eruptive craters. This ultra-long-period event may represent volcanic fluid migration from depth to the surface through a vertical crack during the eruption.


Introduction
Very-long-period (VLP; 2-100-s period) signals are sometimes observed associated with volcanic activity. Their source mechanisms have been investigated at various volcanoes (e.g., Chouet et al. 2003Chouet et al. , 2005Ohminato et al. 2006;Nakamichi et al. 2009;Haney et al. 2013;Maeda et al. 2015;Jolly et al. 2018). Especially long-period signals, with a period exceeding 100 s, were observed during the phreatic eruption of Mount Aso on September 22, 1994 (Kawakatsu et al. 2000). In addition, Sanderson et al. (2010) observed signals coincident with eruptions at Santiaguito volcano. However, observations  (NIED). The triangles and the blue square represent the summits and a recently active crater lake, respectively. The red square indicates the region depicted in Fig. 2. The inset indicates the location of Mount Kusatsu-Shirane at a regional scale of signals in this frequency band at volcanoes are relatively rare compared with VLP signals. These signals were referred to as ultra-long-period (ULP) signals by Chouet (1996). ULP signals are thought to represent pressure changes in the fluid in the source or in a rising explosive magma batch (e.g., Kawakatsu et al. 2000;D' Auria et al. 2006). The analysis of VLP or ULP waveforms can lead to estimations of the source location, geometry, and time history of the force exerted during the volcanic activity.

Graphical abstract
Mount Kusatsu-Shirane is an active quaternary volcano located in the northwestern part of the Gunma Prefecture, Japan (Fig. 1). This volcano consists of three pyroclastic cones: Shirane, Ainomine, and Motoshirane. Since 1882, the eruptive activity has primarily occurred at the Shirane pyroclastic cone, which includes the Yugama crater lake . A hydrothermal system has developed beneath the volcano at Mount Kusatsu-Shirane, and fumaroles exist around Yugama crater. So far, the resistivity structure beneath the volcano has mostly been estimated around the Shirane pyroclastic cone. Nurhasan et al. (2006) conducted an audio-frequency magnetotelluric survey at Mount Shirane and detected a 300-1000-m-thick conductive layer that underlies a resistive layer located on the eastern slope. This conductor was interpreted as being a smectite-rich layer within Pliocene rocks; this layer acts as a low-permeability cap and separates the fluids supplied to the hot springs and fumaroles in the summit area from those on the eastern flank of the volcano (Matsunaga et al. 2020;). In addition, Matsunaga et al. (2020) conducted broadband (320-0.0005 Hz) magnetotelluric surveys in 2015 and 2016 around Mount Motoshirane. They detected a highly conductive layer at depths of 1-3 km, extending from beneath Yugama crater to Mount Motoshirane. They interpreted this layer as a hydrothermal fluid reservoir that supplies fluid to the Yugama crater lake and the hot springs on the eastern and western flanks of Mount Kusatsu-Shirane (Matsunaga et al. 2020).
Seismicity at the Mount Kusatsu-Shirane volcano, especially beneath Yugama crater, became active in 2014 and then ceased in the middle of 2015 (Japan Meteorological Agency (JMA)2018). Then, on January 23, 2018, a phreatic eruption suddenly occurred at Mount Motoshirane, 2 km south of Yugama crater. Mount Motoshirane experienced repeated magmatic eruptions until approximately 1500 years ago (Ishizaki et al. 2020). The eruption onset time is estimated to be 10:02:09 Japan Standard Time (JST) from the records of the air-shock . The phreatic eruption formed eruptive craters lining up in the west-northwest-east-southeast (WNW-ESE) direction on the northern part of Kagamiike-kita crater, also called the main crater (Kametani et al. 2021) (Fig. 2). In addition, there is one small eruptive crater on the west side of Kagamiike-kita crater, called the west crater, and at least four small eruptive craters on the southwest side of Kagamiike crater (The University of Tokyo 2018), called the south crater (Kametani et al. 2021) (Fig. 2). These craters were detected by airborne synthetic-aperture radar (SAR) images (Geospatial Information Authority of Japan (GSI) 2018a) and aerial observations (The University of Tokyo 2018).
At the time of the eruption, short-period borehole seismometers, surface broadband seismometers, and borehole tiltmeters installed around Mount Motoshirane detected signals of volcanic earthquakes, tremors, and tilt motion. The volcanic tremors and tilt motion started 3 and 2 min prior to the eruption, respectively. The records of the broadband seismometer consisted of a ULP signal that exceeded the natural period of the seismometers and resembled those of the tiltmeters. The volume changes associated with the eruption were estimated from eight tiltmeters installed by the Tokyo Institute of Technology and the National Research Institute for Earth Science and Disaster Resilience (NIED) using the model of Okada (1992). The result was a 0.35-m opening (10:00-10:02 JST) and a 0.25-m closing (10:02-10:20 JST) of a 1700-m-long and 850-m-wide crack at 675-1525 m above sea level (a.s.l.) just below the new crater row . The volume changes consisted of a 5.1 × 10 5 -m 3 expansion at 10:00-10:02 JST and a 3.6 × 10 5 -m 3 contraction at 10:02-10:20 JST. Himematsu The green lines indicate the outlines of the eruptive vents estimated by the Geospatial Information Authority of Japan (GSI 2018b) et al. (2020) detected coeruptive and posteruptive crustal deformation associated with the phreatic eruption based on an L-band SAR data time-series analysis. They estimated that the sources that caused the deformation were a plane dislocation and a point source deflation. Yamada et al. (2021) estimated the sources of the volcanic tremor immediately prior to the eruption, using the amplitude source location method (Battaglia and Aki 2003;Kumagai et al. 2010), to be located 1 km north of Mount Motoshirane at a depth of 0.5-1 km. Yamada et al. (2019) applied the waveform inversion method of Nakano and Kumagai (2005) to the broadband seismic records. They estimated an opening crack at a depth of 700 m beneath Kagamiike-kita crater, from the 09:59-10:02 JST data with a frequency band of 0.02-0.1 Hz, and an inclined cylindrical seismic mechanism 1000 m from Kagamiikekita crater, from data of the 80 s following 10:02 JST with a frequency band of 0.033-0.1 Hz. However, the source mechanism of the ULP signal has yet to be investigated.
In this study, we estimate the source mechanism of the ULP event observed during the eruption of Mount Kusatsu-Shirane using waveform inversions of the records of the broadband seismometers and tiltmeters.

Data
The seismic event associated with the eruption of Mount Motoshirane was recorded by broadband seismometers (with a natural period of 240 s) and borehole tiltmeters installed at three stations, Hoshimata (N.KSHV), Nikenya (N.KSNV), and Yazawahara (N.KSYV), of the Fundamental Volcano Observation Network (V-net) (NIED 2019). As shown in Fig. 1, the V-net stations surround Mount Motoshirane and are located approximately 4-6 km from the eruptive craters. They are distributed at altitudes of 1000-1500 m a.s.l. Broadband seismometers are installed at the surface, while tiltmeters are installed at a depth of 200 m. The data were sampled at 100 Hz (broadband seismometers) or 20 Hz (tiltmeters), telemetered to NIED, and published on their website. In this paper, the station codes of the borehole tiltmeters are denoted as N.KSHVB, N.KSNVB, and N.KSYVB to distinguish them from the surface broadband seismic stations.
The eruptive seismic event started at 09:59 JST (Fig. 3). A low frequency (< 0.1 Hz) translational motion started at 10:00 JST and reached a maximum amplitude prior to the eruption (Fig. 4). The east-west (EW) components at N.KSHV and N.KSNV showed a prominent eastward movement, while the EW component at N.KSYV showed a prominent westward movement. The north-south (NS) components at the three stations are rich in long-period signals. A comparison of the raw (Fig. 3) and low-passfiltered (0.1 Hz) waveforms ( Fig. 4) indicates that the waveforms were significantly affected by long-period components. The dominant period of the waveforms during the event exceeds 200 s, indicating that it was a ULP event.
We used three-component broadband seismograms and two-component borehole tiltmeter records at the V-net stations and carried out waveform inversions of the ULP event. From individual raw seismic records covering the time period of 09:30-10:30 JST, we applied a low-pass filter of 10 s and selected a time window of 09:58:20-10:10:00 JST (grey shaded region in Fig. 4). After selecting the time window, we applied a cosine taper to the last 200 s of the window. Even though we did not apply a high-pass filter, we used waveforms that were not corrected for the instrumental response (refer to the following subsection "Waveform inversion"), meaning that the passband was effectively limited to the natural period of the seismometer (240 s).
The horizontal components of broadband seismometers may depict tilt motion in addition to translational motion (Rodgers 1968). In the periods beyond the natural period of the broadband seismometer, responses to tilt motions become larger than those to translational Fig. 4 Low-pass-filtered records. a Low-pass-filtered (10-s) velocity seismograms recorded by broadband seismometers, and b low-pass-filtered (10-s) tilt records converted into velocity waveforms. The red line and the gray band have the same meanings as those in Fig. 3 motions (Aoyama 2008). Because this occurs for the analyzed event, the signals appear to be affected by changes in the tilt. The waveforms of teleseismic event that occurred in Alaska (M 7.9) did not reveal any significant amplitude differences among three-component broadband seismograms, so we did not normalize them as having no local site effect.

Waveform inversion
A waveform inversion for many source parameters requires many observation stations. However, available seismic and geodetic data are limited. Therefore, we follow the approach of Nakano and Kumagai (2005), in which the source mechanism is assumed to reduce the number of free parameters to account for the limited number of stations. Here, a tensile crack source is assumed, because it is the most commonly estimated source mechanism for VLP events at volcanoes (e.g., Chouet and Matoza 2013). The normal vector of the crack was represented by two angles, θ and ϕ , which are measured from the vertical and counterclockwise from the east, respectively. The free parameters to be solved are the centroid source location, θ , ϕ , and the source time function (M(t)).
To estimate the source time function, M(t), we applied the waveform inversion method of Maeda et al. (2011), which accounts for tilt motions. This method uses the representational theorem, which can be written in the frequency domain for a point source as where ω is the angular frequency, p is the direction of the force, q is the direction of the arm of the moment, n is the station and component of the seismic displacement, U is the observed seismogram, M is the source time function, G trans is the Green function of the translational component, G tilt is the Green function of the tilt component, I trans is the translational response of the seismometer, I tilt is the tilt response of the seismometer, and x, y and z, are positive directions to the east, north, and up, respectively. M pq (ω) was obtained via the least-squares method using this equation, and the source time function ( M pq (t) ) was obtained via the inverse Fourier transformation of M pq (ω) . To obtain a single optimum source location and crack orientation independent of the frequency, we calculated the root mean square (RMS) residual value in the time domain summed over the analysis time window, even though the waveform inversion at each candidate (1) p, q = x, y, z source location was conducted in the frequency domain. The formula for this RMS residual is where u obs n (k�t) and u syn n (k�t) are observed and synthetic velocities at the same time sample k of the same trace n , respectively. We carried out grid searches for the location and crack orientation to investigate their optimal values to minimizeE total . The intervals and ranges of the grid searches are summarized in Table 1.
We calculated Green's functions using a finite-difference method, which considers a free-surface boundary condition, including a three-dimensional (3D) topography (Ohminato and Chouet 1997) and perfectly matched layers for the boundary (Festa and Nielsen 2003). We assumed a compressional wave velocity ( V p ) of 2900 m s −1 , a shear wave velocity ( V s = V p / √ 3 ) of 1674 m s −1 , and a rock density ( ρ ) of 2300 kg m −3 (Onizawa et al. 2005). The period range used in this analysis (> 10 s) corresponds to wavelengths of more than 16.74 km. Given these long wavelengths, the effect of small-scale structural heterogeneities is negligible. We used a mesh size of 20 m × 20 m × 20 m, consisting of a 3D mesh with 461 × 461 × 241 nodes, 441 × 506 × 241 nodes, and 505 × 461 × 241 nodes for the calculations of Green's functions using the reciprocity theorem at N.KSHV, N.KSNV, and N.KSYV, respectively. We used a time step of 0.002 s. The 3D topography used in the calculation was based on a 0.4-s grid digital elevation model produced by GSI.
The maximum amplitude of the NS component at N.KSHV appeared at a different time from those of the other five horizontal waveforms at the three stations ( Fig. 4). The cause of this time lag is unknown. However, given the long wavelength of the signal, it is unlikely that only one component would behave differently from the others, implying that the time difference might be due to an equipment malfunction or local effects at the installation site. In addition, it is unclear whether using seismometers alone or both seismometers and tiltmeters is better for the analysis. Taking the above situation into account, we performed the waveform inversion using the following four combinations of data: (1) the three-component broadband seismometers installed at the three stations (9 components in total), (2) the same as (1) but  The tiltmeter records were converted to velocity responses to align their dimensions with those of the seismograms and to treat them equally in the inversion analysis. This conversion was realized by first multiplying the tilt record by the gravitational acceleration ( g = 9.8 m s −2 ) to obtain the apparent acceleration due to the tilt, then integrating the resulting acceleration seismogram to obtain the apparent velocity, and finally multiplying this by the response of the broadband seismometer.
The above process followed that of Maeda et al. (2017). The waveforms we used in this analysis have a dominant period longer than 200 s, which suggests that the response characteristics of the tiltmeter have little influence.

Results
Here, we show the spatial distribution and crack orientation dependence of the residuals between the observed and synthetic waveforms in Figs. 5, 6, 7, and 8, as explained in the following subsections in detail. The origin of the 3D coordinate system was set at sea level at 138.5394° E, 36.6283° N. This horizontal coordinate position is near the center of the main crater row. The results show that the source location and the crack orientation cannot be uniquely identified. Instead, waveform inversions using each of the four data combinations result in 2-3 candidate sources in the grid search (11 candidates in total; Table 2). Below, we explain the results from the four individual data combinations.

Broadband seismograms (9 components)
Three candidate source locations with similar residual values, labeled (a), (b), and (c) in Fig. 5, were obtained from the waveform inversion using all of the three-component broadband seismograms at the three stations (combination (1) in the subsection "Waveform inversion"). The locations of all the candidates are separate from the crater center. Candidate (a) is located near the western end of the eruptive craters. The region of small residuals around this candidate is extremely narrow. This solution indicates a shallow-dipping crack that opens in the northeast (NE)-southwest (SW) direction. Candidates (b) and (c) indicate southern locations of the eruptive craters and show cracks opening in the NE-SW and north-northeast-south-southwest directions, respectively. Candidate (a) is the only solution that explains the time difference in the maximum amplitudes between the NS component of N.KSHV and the other components (Additional file 1: Figs. S1-S3).

Broadband seismograms (8 components)
Three candidate source locations with various distances from the craters, labeled (d)

Broadband seismograms and tiltmeter records (15 components)
Two candidate source locations at various distances from the craters, labeled (g) and (h) in Fig. 7, were obtained by the waveform inversion using all of the three components of the broadband seismometers and the two horizontal components of the tiltmeters at the three stations (combination (3) in the subsection "Waveform inversion"). Candidate (g) is located southwest of the eruptive craters, and the small residual solutions around this candidate are limited to

Broadband seismograms and tiltmeter records (14 components)
Three candidate source locations with various distances from the craters, labeled (i), (j), and (k) in Fig. 8, were obtained by the waveform inversion using eight traces of the broadband seismometers (excluding the NS component at N.KSHV) and all six traces of the tiltmeters (combination (4) in the subsection "Waveform inversion"). Candidate (i) is located approximately west of the eruptive craters, and the small residual solutions around this candidate are limited to shallow-dipping NE-SW-opening cracks on the western side of the crater. Candidate (j) is located southeast of the eruptive craters, and the small residual solutions around this candidate are distributed broadly in space but the orientation is restricted to a shallow-dipping NE-SW-opening crack. Candidate (k) is located away from the eruptive craters to the south. The

Source time functions
The source time functions for these candidate sources show similar temporal variations with different maximum amplitudes (Fig. 9). They show rapid inflation followed by slow deflation of the crack. The inflation starts approximately 2 min prior to the start of the eruption and switches to deflation 5.7-15.3 s prior to the start of the eruption. The maximum moment amplitudes differ between the candidates, ranging from 1.56 × 10 15 N m to 1.89 × 10 16 N m (Table 2).

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), as well as hypocenters of earthquakes (JMA 2018). They are also separated 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 far from the hydrothermal reservoir (Matsunaga et al. 2020). However, this source has some distance from the eruptive craters and appears to be inconsistent with the tiltmeter records near Yugama crater , 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 (Additional file 1: Fig. S1). 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 where and µ are the Lamé constants and �V (t) 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 where L, W, and 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) = �V (t) . We divided the maximum amplitude of the source time function ( N max = max{S(t)}) by = µ = ρV 2 p /3 (6.44 × 10 9 Pa) to investigate the maximum value of the volume change ( V ). To estimate L , W , and d , we need to assume the ratios W /L and L/�d . Taguchi et al. (2018) investigated these ratios for shallow long-period Kusatsu-Shirane sources using a frequency analysis. The  (Geshi et al. 2020). Based on these observations, we assumed W /L = 0.7 and L/�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 ( V ) of 1.04 × 10 6 m 3 , a width ( W ) of 1.15 × 10 3 m, a length ( L ) of 1.65 × 10 3 m, and an aperture change ( d ) of 5.49 × 10 −1 m from the maximum amplitude of the source time function, N max = 6.73 × 10 15 N m (Table 2). This source is north dipping (Additional file 1: Fig. S12), 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 V = 3.38 × 10 5 m 3 , W = 7.92 × 10 2 m), L = 1.13 × 10 3 m, and d = 3.77 × 10 −1 m based on the maximum amplitude of the source time function, N max = 2.18 × 10 15 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 (Additional file 1: Figs. S4 and S7). Subsequently, there is no difficulty in considering candidates (d) and (g) as viable sources.
The range of small residual values around each candidate varies from the other. Several characteristics review the broadness of the small residual value around each candidate. Around candidates (b), (c), (e), (f ),   Table 4 Summary of the second waveform inversion results.
The notation is the same as in Table 2 Candidate ( and (k), which are located in the center of the observation network, there is a wide spatial distribution of small residual values. On the other hand, the distribution of small residual values for crack orientation at each candidate location is narrow. A similar trend can be seen for candidates (i) and (j) between N.KSHV and N.KSYV. Around candidates (a) and (h), which are located outside the observation network, the distribution of small residual values for crack orientation at each candidate location is also narrow. Significantly, the spatial distribution of small residual values for space is extremely narrow around candidate (a). These candidates were treated as unlikely candidates in the above discussion, and their characteristics are thought to be evident in the way the residual values are distributed.
In other words, a small distribution of residual values in space or crack orientation is considered to represent the unreliability of the solution. In candidates (d) Therefore, an N-S-opening sub-vertical tensile crack near the eruptive craters (Additional file 1: Fig. S13) 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).

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 Fig. 12 Location and geometry of the tensile crack for candidate (d1). The red arrows indicate the opening direction. The stars represent the centroid of candidate (d1) Fig. 13 Waveform fitting between the observed (red) and synthetic (blue) waveforms for candidate (d1) 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 (Additional file 1: Figs. S4 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 (Additional file 1: Figs. S7 and 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.

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 . 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. Although our solution may show the same source by Terada et al. (2021), there are several differences between the two solutions. First, the result of Terada  (d1) and (g1). The notation is the same as that used in Fig. 9 et al. (2021) constrained a solution uniquely. Second, Terada et al. (2021) used a tilt data set, whereas we used seismic and tilt data sets. Third, Terada et al. (2021) artificially separated the two-time windows, whereas the time window we used did not separate the phases. The time of maximum amplitude of the source time function obtained in this study and the time of the boundary divided by Terada et al. (2021) are different. Forth, the number of stations used in our analysis was less than that used in Terada et al. (2021). As a result, the location and mechanism of our Spatial distribution and crack orientation dependence of the residuals between the observed and synthetic waveforms obtained by a second grid search about candidate (g). The locations of candidates (g) and (g1) are indicated by straight lines on the plane map. See Fig. 5 for details concerning the plotting format solution may be less constrained than those of Terada et al. (2021), while transition from expansion to contraction may be better constrained in our study.

Conclusions
We conducted waveform inversions of the ULP event associated with the phreatic eruption of Mount Kusatsu-Shirane on January 23, 2018. Eleven candidate source locations were obtained from the waveform inversion, indicating a limitation of ULP source analyses using few seismic stations. After examining the various likelihoods of the 11 obtained candidates, 9 of these candidates were classified as unlikely based on comparisons with observations and 2 of the candidates were similar to each other. Therefore, a suitable source candidate appears to be an N-S-opening sub-vertical crack near the eruptive craters. This ULP event may represent volcanic fluid migration at a velocity of 49-100 m s −1 from depth to the surface through this crack during the eruption.  Fig. S1. Waveform fitting between the observed (red) and synthetic (blue) waveforms for candidate (a). Fig. S2. Waveform fitting between the observed (red) and synthetic (blue) waveforms for candidate (b). Fig. S3. Waveform fitting between the observed (red) and synthetic (blue) waveforms for candidate (c). Fig. S4. Waveform fitting between the observed (red) and synthetic (blue) waveforms for candidate (d). Fig. S5. Waveform fitting between the observed (red) and synthetic (blue) waveforms for candidate (e). Fig. S6.. Waveform fitting between the observed (red) and synthetic (blue) waveforms for candidate (f ). Fig.  S7. Waveform fitting between the observed (red) and synthetic (blue) waveforms for candidate (g). The notation is the same as that in Fig. 16.  Fig. S8. Waveform fitting between the observed (red) and synthetic (blue) waveforms for candidate (h). The notation is the same as that in Fig. 16.  Fig. S9. Waveform fitting between the observed (red) and synthetic (blue) waveforms for candidate (i). The notation is the same as that in Fig. 16. Fig.  S10. Waveform fitting between the observed (red) and synthetic (blue) waveforms for candidate (j). The notation is the same as that in Fig. 16. Fig.  S11.. Waveform fitting between the observed (red) and synthetic (blue) waveforms for candidate (k). The notation is the same as that in Fig. 16.  Fig. S12.. Location and geometry of the tensile crack for candidate (e). The notation is the same as that in Fig. 12. Fig. S13. Location and geometry of the tensile crack for candidate (d). The notation is the same as that in Fig. 12.