S-wave modelling of the Showa-Shinzan lava dome in Usu Volcano, Northern Japan, from seismic observations

∗ akiko-t@eri.u-tokyo.ac.jp 1 Earthquake Research Institute, The University of Tokyo, Tokyo, Japan 2 Institute of Seismology and Volcanology, Hokkaido University, Hokkaido, Japan 3 Now at Tokyo Metropolitan Government Bureau of Construction, Tokyo, Japan 4 Graduate School of Science and Technology, Hirosaki University, Aomori, Japan 5 Now at JX Nippon Oil & Gas Exploration Corp, Tokyo, Japan 6 Now at Kobe University, Hyogo, Japan 7 Department of Surveying and Geo-Informatics, Southwest Jiaotong University, Chengdu, China 8 Department of Earth and Space Sciences, Southern University of Science and Technology,

the strong inhomogeneity within the small area and additional technical difficulties. We therefore focus on one-dimensional models to mainly discuss the relationship between S-wave velocity and density in this area (Section 6.3), discuss the detailed difficulty of obtaining a three-dimensional model (Section 6.4), and leave the three-dimensional modelling for a separate study.

USU VOLCANO AND THE SHOWA-SHINZAN LAVA DOME
Usu Volcano is one of the arc volcanoes in the Japanese islands associated with the subduction of the Pacific plate (Fig. 1). The volcano sits near the rim of the Toya caldera ( Fig. 1) that was formed circa 112,000-115,000 years ago (e.g., Miyabuchi et al. 2014). Usu Volcano formed about 30,000 years ago after a quiescence of Toya caldera for more than 60,000 years. After its formation, the volcano was dormant until the first eruption in historical time in 1663 with a Volcanic Explosive Index (VEI; Newhall and Self 1982) of 5. Recent eruptions after 1663 took place in 1769 (VEI=4), 1822 (VEI=4), 1853 (VEI=4), 1910 (VEI=2), 1943-1945(VEI=2), 1977-1982, and 2000 (VEI=2). Volcanism of Usu Volcano is characterized by the emergence of silicic magmatism with 68-74 wt of SiO 2 (e.g., Tomiya and Takahashi 2005) not only at the summit but also on the flank. Indeed, while lava domes associated with 1977-1982 eruption emerged at the summit, those associated with 1910, 1943-1945, and 2000 eruptions emerged from northern, eastern, and western flanks, respectively. The geological map of Usu volcano is available in Soya et al. (2007).
The Showa-Shinzan lava dome emerged during the 1943-1945 eruption. It is the only dome with a mound among lava domes that recently emerged in Usu Volcano, consisting of a dome outcrop without vegetation around the summit with a maximum elevation of 398 m, north (northwestern) and east roofs with vegetation at elevations of ∼250 m ( Fig. 1c and Fig. 2). Based on seismic and geodetic observations associated with the eruption by Minakami et al. (1951) and Kizawa (1957Kizawa ( , 1958, Yokoyama and Seino (2000) and Yokoyama (2002Yokoyama ( , 2004 suggested that 1943-1945 unrest is divided into four stages: intense seismicity, magma migration, explosions, and dome forming. The eruption was preceded by a swarm of earthquakes from 28 December 1943 at Yanagihara (YH; Fig. 1b), with the largest earthquake of M 5.0 on 9 January 1944 at a depth of 11 km Showa-Shinzan S-wave Modelling 5 (Kizawa 1958). Seismicity was highest in the first period of the seismic swarm, waning over time. Yokoyama and Seino (2000) and Yokoyama (2002Yokoyama ( , 2004 suggested that intruded magma up to ∼1 km below sea level drove the seismic swarm. They also suggested that an eruption did not take place at the time because the depth at the top of the intruded magma is far below that of the aquifer. Indeed, the earthquakes there in July 1944 are located at depths between 1.5 and 5.0 km (Minakami et al. 1951). The area around YH was uplifted by 23 meters in the first four months since the onset of the seismic swarm and by 60 meters by May 1945 (Minakami et al. 1951).
The magma then migrated northward by about 2 km to Fukaba (FB; Fig. 1b) before the first eruption on 23 June 1944 to generate a mound there. The eruption was phreatomagmatic, leading Yokoyama and Seino (2000) and Yokoyama (2002Yokoyama ( , 2004 to suggest that the eruption resulted from an interaction of ascended magma with groundwater. Eruptions continued until early December 1944, when a solidified lava spine appeared at the top of the mound. A series of eruptions gave way to dome forming by the end of 1944. The dome forming uplifted the surface with a rate up to 0.5 m/day between May and October 1945 (Yokoyama and Seino 2000). By the end of 1945, the lava dome and the mound grew to the height of 280 and 170 m from the base with the lava dome's diameter approximately 300 m. The total volume of the emerged lava dome is about 2.2 × 10 7 m 3 above the mound, and the total volume of deformation is 2.3 × 10 8 m 3 .
After the eruption, the dome has been subsiding at least since the 1960s (e.g., Yokoyama and Seino 2000;Wang and Aoki 2019). Recent space geodetic observations reveal that Showa-Shinzan has been subsiding more or less constantly with a rate of about 16 mm/yr in the last ∼25 years Wang and Aoki 2019). Wang and Aoki (2019) suggested that this subsidence is due to the thermal contraction of the emerged lava dome. Indeed, the heat brought at the time of the dome emergence has not been wholly released; the temperature of the fumarole was still 190 • C as of early 1997, decreasing from 1000 • C in 1945 and ∼500 • C in 1960 (Minakami et al. 1951;Yokoyama and Seino 2000).
The temperature of 1000 C • at the fumarole right after the emergence of the dome suggests that the intruded lava must have been at least partially molten because it is the temperature above the solidus of dacite (e.g., Holtz et al. 2001). On the other hand, the river gravel on the intruded lava dome suggests that its surface must have been solidified or viscous enough to carry the river gravel upon intrusion.
The structure of the Showa-Shinzan lava dome has been investigated for decades. An early comprehensive survey by Nemoto et al. (1957) found a local positive Bouguer gravity anomaly of up to 1 mgal in Showa-Shinzan. They also inferred that the lava dome has a P-wave velocity as high as about 4 km/s at around sea level or ∼400 m beneath the summit, in contrast to 1.7-1.9 km/s around it at the same depth level. These observations suggest that the intruded material during the 1943-1945 eruption is more rigid than the surroundings. A few decades later, Miyamachi et al. (1987) utilized explosions by a firework in 1984 to find that the P-wave velocity inside the lava dome decreased to 1.8-2.2 km/s, although the velocity is still higher than surrounding bodies.
They interpreted this velocity decrease as due to thermal metamorphism. Nishida and Miyajima (1984) suggested the existence of a column-like intruded magma above the Curie temperature (560 • C) based on magnetic surveys. Geochemical measurements of fumaroles led Symonds et al. (1996) and Hernández et al. (2006) to suggest that the mound, as well as the lava, consists of magma intruded during the 1943-1945 eruption. However, Goto and Johmori (2014) suggested from the resistivity structure that the sizable magma with a diameter of ∼400 m exists only beneath the lava dome as a quasi-spherical body, and the mound does not contain magma. This suggestion is consistent with Nishida and Miyajima (1984), who did not find any high-temperature body above the Curie temperature beneath the mound.
Recent developments in muon radiography allow us to delineate detailed density variation beneath the lava dome. Tanaka et al. (2007) and Tanaka and Yokoyama (2008) delineated a twodimensional density variation beneath the lava dome from a single muon detector to find a highdensity body with a diameter of about 400 m beneath the lava dome. This high-density body represents the intruded magma, consistent with the resistivity structure (Goto and Johmori 2014).
More recently, Nishiyama et al. (2014) combined muon radiography with gravity measurements at the surface to make the density distribution three-dimensional. Nishiyama et al. (2017) further added a second muon detector to improve the precision of the image. They imaged a subspherical high-density body with density up to ∼ 2.6 × 10 3 kg/m 3 which is inferred as the lava dome.
Since the S-wave velocity models in the lava dome have not been estimated in previous studies, we conducted a passive seismic survey to obtain S-wave velocity models and compare them with density models obtained in previous studies.

CHARACTERISTICS
We deployed 21 1-Hz seismometers, U01-U21 (Lennartz 3D Lite Mark III; Fig. 1c and Table 1), around the Showa-Shinzan lava dome on 8-10 May 2018 and retrieved them on 6-8 June 2018. Seismic waveforms were recorded by Hakusan LS8800 data loggers with a dynamic range of 24 bits and a sampling frequency of 200 Hz. The seismometer and data logger's energy consumption is so low that a seismic site could survive for about a month with eight Panasonic EVOLTA D batteries under the conditions of daily minimum and maximum atmospheric temperature of 5-10 • C and 15-20 • C, respectively. Each seismometer was buried, when possible, to avoid cultural or unwanted noises as much as possible. We also utilized the seismic site V.USSW ( Fig. 1c and Table   1) with another 1-Hz seismometer operated by Japan Meteorological Agency for the analysis. Table 2 show the nominal frequency response of 1-Hz seismometers that we deployed in the campaign (Lennartz 3D Lite Mark III). There is no frequency response information for each instrument provided by the production company, unlike the cases in broadband seismometers. Since the typical interstation distance is short of ∼400 m comparable to the wavelength at 2 Hz, 20-degree relative phase-response error can cause ∼ 20/360 = ∼ 6% error in phase-velocity measurements. If the uncertainty of natural frequency itself is large, for example, the error in phase response can be especially large around the natural frequency where phase dependence on frequency is large. For these reasons, we calibrated instrumental responses of 1-Hz seismomters at Earthquake Research Institute, the University of Tokyo, between 26 March and 10 April 2018. A broadband seismometer (Güralp CMG-3T) with a natural frequency of 120 s is used as a reference.

While Figures 3a, b and
We used a one-day record on 27 March 2018 because the day was seismically quiet. First, we divided the whole one-day record into 10-minute segments. Second, a segment is accepted for further analysis only if the maximum amplitude of the seismogram of each 1-Hz seismometer is between 0.5 and 2 times that of the broadband seismometer. This procedure is to avoid glitches in the seismogram of the 1-Hz seismometers. Third, the Hanning taper is applied to the accepted segments. With an assumption that the broadband seismometer recorded the true ground motion, the spectrum of the i th 1-Hz seismometer at the j th segment X j i (f ) can be represented bȳ where Y j (ω) represents the spectrum of the referenced broadband seismometer at the j th segment.
R i is the response function of the i th seismometer given by where h is the dumping constant, ω 0 is the natural frequency of about 1 Hz, ω 1 is the other cut-off frequency, and τ is time delay. We perturb these parameters from the initial values; the initial value of ω 0 , ω 1 and h are given by the catalog (h = √ 2/2, ω 0 = 2π (rad/s), and ω 1 = 1.083 (rad/s)). We also assumed that the other initial values as τ = 0.015 (s), and A = 1.0.
We inferred five parameters, ω 0 , ω 1 , h, A, and τ by minimizing the squared spectrum difference S between that from the 1-Hz seismometers and predicted from the broadband sensor over the frequency band of ω min = 0.2π (0.1 Hz) to ω max = 3π (1.5 Hz): Figures 3c-h denote the amplitude and phase responses of the seismometers to the nominal response in the catalogue, indicating that the phase difference among seismometers reaches about 5 • between 0.1 and 1 Hz. Although the corresponding phase-velocity error is not so large of 1-2%, correction of the frequency response for each sensor enables us to measure the phase with a precision smaller than 1 • , and allows us to utilize the phase information in the ambient noise cross-correlations.
For understanding the characteristics of ambient noise in our dataset, we plotted the spectral contents of each seismic record. Figure 4 shows an example of spectral density functions of seismic sites at the summit (U06), at the roof (U18), and at the base (U07), which are least to most susceptible to cultural activity. They indicate that all sites recorded signals down to ∼0.1 Hz, suggesting Showa-Shinzan S-wave Modelling 9 that extracting the wavefield higher than 0.1 Hz is possible. The noise level of all sites between 0.1 and 10 Hz is between the low and high noise models (Peterson 1993) for all components (Fig. 4). Figure 5 shows spectrograms of the same sites as Figure 4, clearly indicating that all sites exhibit diurnal variations in noise amplitude, especially below 10 Hz. These variations are likely due to cultural activitybecause the Showa-Shinzan area is a sightseeing area, and human activity will not be lower at the weekend compared to weekdays. The noise levels at these sites are comparable during the daytime, but those at U18 and U07, which are located closer to cultural activity, are higher than U06 at night. This observation might indicate that cultural activity is not negligible, even at night. The strong noise, on May 18-21, for example, have a correlation between below 1 Hz and above 10 Hz. This might correspond to natural noise such as oceanic disturbances below 1 Hz and wind disturbances above 10 Hz.

PHASE-VELOCITY MEASUREMENT
This section describes phase-velocity measurements, including the methods. To measure phase velocities, we use one of the conventional seismic interferometry methods but with assuming layered structures. As a result, we obtain average phase velocities of the fundamental-mode Rayleigh wave at three areas: summit, roof (east roof in Figure 2) and base.
To extract seismic wavefield between two seismic sites, we take spatial auto correlation (SPAC) of seismic records (e.g., Aki 1957;Takeo et al. 2013). The SPAC method is very similar to seismic interferometry in terms of calculating cross spectra, which is equivalent to calculating crosscorrelation functions to extract empirical Green's function between stations (Shapiro and Campillo 2004) by seismic interferometry. Unlike phase-velocity measurement for each pair of stations in ambient noise tomography, however, the SPAC method will measure phase-velocity velocities for each group of seismic stations.
For the calculation of cross spectra, the vertical components of the whole dataset were first down-sampled to 20 Hz and divided into segments of 102.4 seconds with overlaps of 51.2 seconds.
We then rejected segments whose mean-square amplitude at 1-5 Hz is more than 10 times that of the former segment to avoid contamination of deterministic signals such as earthquakes. After these procedures, we took ensemble averages of the cross spectrum of vertical records of i th and j th site by spectral whitening method (Bensen et al. 2007) as where F i (ω) denotes the seismic record of the i th site in the frequency domain, ω represents angular frequency, * denotes the complex conjugate, and denotes the ensemble average of 102.4s segments.
Figures 6a, 7a, and 8a represent extracted wavefield in time domain aligned with interstation distance at the summit, roof, and base areas, respectively, showing that the wave propagation within at least a few hundred meters can be extracted at the summit and roof areas, and more than one kilometre at the base. For clarification of the signals, waveforms for each region are filtered at the range where uncertainties of phase-velocity measurements (described later) are lower than 0.1 km/s. We should note that the quality of each waveform is not enough to measure phase velocity for each pair of stations and conduct ambient noise tomography. The short interstation distance equivalent to the wavelength also makes the measurement for each pair difficult. We, therefore, focus on measuring average phase velocities for each areas by the SPAC method (Aki 1957), which is capable of measuring phase velocities in the frequency domain even if the quality of each cross spectrum is low and when the interstation distance is equivalent to the wavelength.
Assuming a laterally homogeneous medium, the cross spectrum (equation 4) between seismograms obtained at two sites depends on interstation distance between the i th and j th site d ij and the phase velocity of the Rayleigh wave c is given by the SPAC method (Aki 1957) as where J 0 is the Bessel function of the first kind at the 0th order. A(ω) is a constant representing the source intensity estimated by a least-squares fitting for each frequency. The residual between the real component of the observed cross spectrum, Re S obs ij (ω, c) , and the modelled cross spectrum, , as a function of phase velocity and angular frequency can be defined by Showa-Shinzan S-wave Modelling 11 where w ij (ω) corrects the amplitude variation due to the normalization in equation (4) as (Takeo et al. 2013) Figures 6c, 7c, and 8c show residuals as a function of frequency and phase velocity. In principle, phase-velocity measurement is possible by connecting those with the smallest residuals at each frequency. However, this approach would make the modelled dispersion physically unrealistic due to 2π ambiguity and low quality of data. Therefore, we assumed a layered structure as follows to obtain physically realistic dispersion curves.
The S-wave velocity structure is modelled with seven layers within the shallowest 1 km, and constant at depths below 1 km ( Table 3). The layer thicknesses are almost equal intervals in log scale by considering the shape of surface-wave sensitivity kernels. The S-wave velocities greater than 300 km and the attenuation structure are fixed to the Preliminary Earth Reference Model (PREM; Dziewonski and Anderson 1981). S-wave velocity scales the P-wave velocity and density, according to Brocher (2005). The phase velocity corresponding to the modelled velocity structure c(ω) = c(ω; β) is calculated by DISPER80 (Saito 1988), where β denotes the model parameter vector, i.e., S-wave velocities at 8 layers. The S-wave velocity at each layer is constrained to be larger than 80% of that of the layer directly above to avoid numerical instability at high frequencies.
The seismic structure of Showa-Shinzan is so heterogeneous that we cannot employ equation (5) to obtain a one-dimensional reference velocity of the whole area of interest. We thus divided the whole region into three areas, summit, roof and base to derive the Rayleigh wave dispersion for each area (Figs. 6b, 7b, and 8b). Figures 6c, 7c, and 8c denote the distribution of R(ω, c) and the measured phase velocities for each area. Figure 9a depicts the phase velocities of all areas. The corresponding tentative S-wave velocity models are shown in Figure A1. Although the frequency ranges of the accurate phase-velocity measurements are different between areas, the corresponding depth range of high sensitivity is similar: 20-500 m below the surface at the summit area (Fig. 6d), 10-500 m below the surface at the roof area (Fig. 7d), and 50-500 m below the surface at the base area (Fig. 8d). The stronger heterogeneity at depths shallower than ∼10-50 m for each area might have caused incoherence phases between different pairs of stations, and large uncertainties of phase-velocity measurements at high frequencies ( Fig. 6c, Fig. 7c and Fig. 8c).

ONE-DIMENSIONAL S-WAVE VELOCITY MODELS
During the phase-velocity measurement, we assumed layered structures (Fig. A1). Since the structures are tentative as described in Appendix A, we obtain final models by only using phase-velocity measurements obtained in the previous section and by further introducing a vertical smoothing parameter, ε, as described in Appendix B. The inversion allowed us to delineate the S-wave velocity structure down to ∼500 m from the surface for each area, as already seen in sensitivity kernels (Figs. 6d, 7d, and 8d). As expected from the phase velocity dispersion of the Rayleigh wave, the summit area exhibits the highest Swave velocity with ∼2 km/s at sea level, whereas the base area exhibits the lowest S-wave velocity with ∼1 km/s at sea level. The S-wave velocity at the roof area is intermediate down to ∼100 m below sea level and comparable to that at the base area at greater depths. These results do not depend on the choice of ε, as shown in Figures A2d-f. At depths greater than ∼500 m from the surface, we recognize the non-negligible effect of vertical smoothing. At this depth range, S-wave velocity becomes lower with ε = 1 (Fig. A2f) than that with ε = 0.1 (Fig. 9b and Fig. A2e) due to insufficient phase velocity measurements at lower frequencies and stronger vertical smoothing.
The boundary between layers is sharp due to our assumption of the layered structure. Most of the boundaries have no physical meaning because the uncertainty of estimated S-wave velocity values overlaps between adjacent layers. The boundary at a depth of 100 m from the surface for the summit region (red in Fig. 9b), however, seems significant compared to the uncertainties even with a lower smoothing parameter (Fig. A2d). The meaning of this boundary will be discussed in Section 6.3.

Phase-Velocity Measurement Method
For the ambient noise analyses at periods longer than ∼10 s, various automatic phase-velocity measuring method has been developed, such as spatial auto-correlation method (Aki 1957) or zero crossing (Ekstrom et al. 2009), and used for regional ambient noise tomography studies. There have been, however, still technical difficulties for the phase-velocity measurement in the highly heterogeneous regions, especially at frequencies higher than ∼0.1 Hz. At high frequencies, the interstation distance becomes much larger than the wavelengths, which causes high attenuation and low SN ratio of the extracted surface-wave signals, and causes severe problems of the phase unwrapping ambiguity (e.g., Liu et al. 2016). The strong dispersion at high frequencies further makes the phase ambiguity serious and hamper the group-velocity measurement. For example, in oceanic regions, previous studies carefully prepared spline functions for each region, whose intervals depend on the frequency, to measure multi-mode phase velocities (Takeo et al. 2013(Takeo et al. , 2014. This study, therefore, developed an original method to measure average phase velocities of surface waves from cross-correlation functions by assuming layered structures. This approach is similar to those for the analyses of teleseismic surface waves at periods longer than ∼30 s to analyze multi-mode waveforms (Yoshizawa and Kennet 2002) or to stably measure phase-velocities from noisy records of ocean bottom seismometers (Takeo et al. 2018).
The new method allowed us to stably measure phase velocities with bootstrap uncertainties less than ∼0.1 km/s (Fig. 9a) without carefully preparing spline functions for each region. It also enabled us to avoided non-physical oscillations due to insufficient quality of data at certain frequencies (recognized as high residual vertical columns in Fig. 6c, 7c, and 8c). The phase unwrapping ambiguity is solved automatically in the frequency range, which can be recognized as the range with low bootstrap uncertainties. On the other hand, large uncertainties in higher frequencies partly correspond to the difficulty of unwrapping even with the assumptions of layered structures.
The method is applied to measure average phase velocities for several pairs of stations (Fig. 6b,   7b, and 8b) in this study. Although the same method is applicable to the phase-velocity measurement for each path, the uncertainty evaluation must be done differently. Moreover, the signal-tonoise ratios of each cross-correlation function (Fig. 6a, 7a, and 8a) are insufficient for the phasevelocity measurement of each path. We focused on averaging phase-velocity measurements at three regions, leaving the three-dimensional tomography for a separate study. The one-dimensional models still allow us to discuss the structure beneath the lava dome and the correspondence between S-wave velocities and densities, as described below.

S-wave Velocity Models
Our results show that the S-wave velocity at the summit area is significantly higher than the surrounding roof and base areas from the top to a depth of ∼200 m below sea level. The shallow highvelocity body above sea level represents the lava dome formed during the 1943-1945 eruption and then slowly cooled. This is consistent with insights gained from previous seismic explorations (Nemoto et al. 1957;Miyamachi et al. 1987), muon radiography ( The average S-wave velocity at the roof area is not as high as the summit area but higher than the base area. This velocity structure may indicate that the first few hundred meters beneath the roof area contain some portion of intruded magma as a high-velocity body, as fumarole temperature measurements suggest (Symonds et al. 1996;Hernández et al. 2006). These results suggest that the narrow root of the lava dome existing directly beneath the summit with a minor effect on the rocks beneath the roof area.
We should also note that our S-wave velocity structure at the base area is roughly consistent with the regional seismic structure in Japan Seismic Hazard Information Station (http://www.jshis.bosai.go.jp/map/?lang=en). In this regional structure model, the S-wave velocity around Showa-Shinzan at a depth of 30 m is 300-400 m/s, and it reaches 1 km/s at a depth of 300-400 m and 2 km/s at a depth of ∼800 m.

Density and S-wave Velocity
We compare our results with the latest three-dimensional density model in the Showa-Shinzan obtained from the combination of muon radiography and gravity measurements by Nishiyama  Figure 9b also have corresponding densities.
We first consider the density of the roof area's surface, where the value is estimated to be 1.5-1.7×10 3 kg/m 3 in N2017 and 1.7-2.0×10 3 kg/m 3 in our model (red squares in Fig. 10a). To evaluate whether models with a low-density body can fit phase velocity measurements, we conducted an inversion by constraining the density at the uppermost layer (0-20 m from the surface) to be 1.6×10 3 kg/m 3 , or the S-wave velocity of 0.27 km/s, and by introducing the smoothing term in equation (B.1) only from the second layer. The obtained structure can still fit phase velocity measurements (green triangles in Fig. 10a), indicating that the near-surface S-wave velocity and density in the roof area is extremely low.
We next consider density in the summit area, which is 2.0-2.1×10 3 kg/m 3 in our model (red in Fig. 10b) and ∼2.1×10 3 kg/m 3 at the surface to ∼2.3×10 3 kg/m 3 inside by N2017. Similar to the roof area, we fixed densities in the top three layers: 2.1×10 3 kg/m 3 at depths of 0-20 m, 2.2×10 3 kg/m 3 at depths of 20-50 m, and 2.3×10 3 kg/m 3 at depths of 50-100 m. We inverted for the S-wave velocities by introducing the smoothing term in equation (B.1) only to the 4 th layer or below. Green triangles in Figure 10b show that the obtained structure does not fit with phasevelocity measurements, indicating that the B2005 scaling between density and seismic velocities is inappropriate to reconcile with our phase-velocity measurements with the N2017 model. Therefore, we performed another inversion by fixing the density at the top three layers, freely searching S-wave velocities without scaling with density, applying the smoothing for all layers. Here the P-wave velocity is scaled with the S-wave velocity by B2005. Blue diamonds in Figure 10b show that the obtained model fits well with phase-velocity measurements.
These results indicate that the S-wave velocity at the summit is ∼1 km/s, ∼45% lower than that expected from the N2007 density model with the B2005 scaling as summarized in Figure   11. One possible cause of the S-wave velocity reduction is cracking in the dome that effectively reduces S-wave velocity compared to density. Hudson (1980Hudson ( , 1981 theoretically demonstrated that the inclusion of infinitely thin fluid-filled cracks reduces only S-wave velocities. Another possible cause is the different rock types. The main composition of Showa-Shinzan is dacite, whereas B2005 scaling is based on a variety of rocks. Gaunt et al. (2016) showed dacite sampled at Mt. St. Helens whose density is 2.5×10 3 kg/m 3 , and S-wave velocity is ∼1.5 km/s (Fig. 11) for a rock with connected porosity of 6%. The S-wave velocity is much lower than that expected from the density and B2005 scaling. This tendency is consistent with our results, although we cannot isolate the effects of rock type and crack density from this sample.
Since we expect high temperatures beneath the summit area from the ongoing fumarolic activities (Fig. 2), we also focus on S-wave velocity reduction of up to ∼15% at high temperature compared to that at room temperature for the dacite sample (Gaunt et al. 2016;Fig. 11). The interpretation of attributing S-wave velocity reduction to high temperature, however, conflicts with our result of sudden velocity increase at a depth of ∼100 m below the surface (red in Fig. 9b and Fig. 10b) because the velocity should not increase at this depth if the temperature is still high in the dome but is cooled from the surface. Higher crack density near the surface is more likely to explain the low S-wave velocity near the surface. We thus conclude that the cracks should explain the inconsistency between B2005 scaling, N2017 density model, and our result.

Prospects
Although we attempted to conduct ambient noise tomography, several technical problems arose.
To obtain a three-dimensional S-wave velocity model, we need to address the following things.
First, we need to improve the SN ratio of surface waves extracted in cross-correlation functions (Fig. 6a, 7a and 8a) to measure phase velocities for each pair of stations. Changing the normalization method during the calculation of cross spectrum or cross-correlation functions might work, but several approaches (e.g., Bensen et al. 2007) do not succeed at this moment. The strong lateral heterogeneity may cause strong scattering, and hence the low SN. Second, we need to overcome the fact that the available frequency range of phase-velocity measurements depends on areas ( Fig. 9a) and even on paths due to short interstation distances and strong lateral heterogeneity compared to the wavelengths. This situation makes the conventional phase-velocity mapping and ambient noise tomography difficult. Third, we need to consider the effect of rugged topography around Showa-Shinzan (Fig. 1c). The conventional correction of topography (Snieder 1986) is not applicable because it assumes a homogeneous S-wave velocity near the surface. The numerical or theoretical consideration for shallow structure with topography (e.g., Wang et al. 2012) is required. To overcome these difficulties, for example, we might need to construct an initial velocity structure with the one-dimensional models obtained in this study and rugged topography, and then iteratively update towards the final velocity model by waveform fitting.
While we observed three-component seismic motion in our seismic campaign, we only used vertical records to extract the propagation of the Rayleigh wave by taking cross-correlations of seismic records. By extracting the spectral ratio between horizontal and vertical components of the Rayleigh wave, the shallow structure might be possible to constrain. Moreover, utilizing horizontal records will allow us to extract the wavefield of the Love wave, although the analysis will be complicated by the ray bending and scattering of the Love wave stronger than that of the Rayleigh wave. Because the Rayleigh and Love wavefield lead to the velocity of SV and SH waves, respectively, utilizing both Rayleigh and Love waves enable us to infer anisotropic characteristics of the lava dome and surrounding materials assuming transverse isotropy. Such information gives us insights into the lava dome's finer structure, such as whether the lava is horizontally or vertically layered (e.g., Jaxybulatov et al. 2014).
The three-dimensional velocity structure thus obtained can be combined with three-dimensional density structures derived from muon radiography and gravity measurements, as we have discussed above. This also has the potential to infer elastic constants, gaining more insights into the material property of the lava dome and surroundings.

CONCLUSION
We conducted a seismic campaign with 21 1-Hz seismometers in the Showa-Shinzan lava dome, which emerged during the 1943-1945 eruption of Usu Volcano, Japan (Section 2). Before the campaign, we calibrated seismometers to be deployed on the field for their frequency response.
The calibration demonstrates that the data down to 0.1 Hz can be used for the analysis in this study. Also, we confirmed that the frequency response of the seismometers varies only slightly Showa-Shinzan S-wave Modelling 19 between them. For example, the phase response varies only less than 1-2 degrees between the seismometers (Section 3). The obtained vertical seismic record was cross-correlated to extract the seismic wavefield between seismic sites. The obtained cross-correlation functions gave us dispersion of the Rayleigh wave in three relatively homogeneous areas (Section 4), then obtained one-dimensional S-wave velocities down to a few hundred meters below sea level (Section 5). The new phase-velocity measurement method enabled us to automatically measure phase velocities in heterogeneous regions by assuming layered structures (Section 6.1). The S-wave velocity in the dome beneath the summit area is higher than surrounding areas by a few tens of percent down to a depth of ∼200 m below sea level, consistent with previous studies and indicating the root of the lava dome (Section 6.2). The obtained S-wave velocity inside the dome of ∼1 km/s is much slower than that expected from a density value of ∼2.3×10 3 kg/m 3 (Nishiyama et al. 2017) and a conventional scaling law (Brocher 2005), indicating the effect of cracking in the shallow lava dome (Section 6.3).
Further studies are required to gain more information on the lava dome's seismic structure from the observed seismic data. For example, taking rugged topography and lateral heterogeneity into account will extract the wavefield more accurately, leading to a delineation of an accurate three-dimensional seismic velocity structure. Also, utilizing both vertical and horizontal records will allow us to extract information on seismic anisotropy, thereby gaining information on a finer structure inside the lava dome and the roof area (Section 6.4).

DATA AVAILABILITY
Cross-correlation functions taken from the obtained seismic data are uploaded at the Zenodo repository (doi:10.5281/zenodo.4313960). . We used SAC2000 (Goldstein and Snoke 2005) during this study. The data analysis was carried out using ObsPy (Krischer et al. 2015).    1943-1945, 1977-1982, and 2000 eruptions. YH and FB denote Yanagihara and Futaba, where the initial and subsequent uplift, respectively, took place before the dome emergence during the 1943-1945 eruption.
Light gray area is the south end of Lake Toya corresponding to the Toya caldera. A rectangle marks the area of (c). (c) Spatial distribution of temporary seismic sites with a contour interval of 10 m.  (Fig. 1c). Ongoing steam can be seen near the summit. The panels plot vertical (U), east-west (E), and north-south (N) components. The reddish colors represent a higher probability. The spectra below 1 Hz are similar because the origin of the noise below 1 Hz is ocean swells. Between 1 and 10 Hz, two peaks for each probability distribution correspond to daytime and nighttime because the primary origin of the spectrum is the cultural noise. The high noise with a dominant frequency of 10 Hz at T.U06 could be attributed to the traffic noise. Thick grey lines are low and high noise models by Peterson (1993).

APPENDIX A: TENTATIVE STRUCTURES
In section 4, we assumed layered structures during phase-velocity measurement. In general, this kind of tentative structure fluctuates strongly due to tradeoff between adjacent layers (e.g., Yoshizawa and Kennet 2002;Takeo et al. 2018). Figure A1 shows tentative structures corresponding to phasevelocity measurements in Figure 9a. Since we assumed a monotonic increase of S-wave velocity with depth, those tentative models' uncertainties are small enough to recognize some structural differences between areas. However, the models are still nominal, and the model uncertainties depend on the choice of the allowed maximum S-wave velocity of a layer relative to that one layer above. In contrast, the modelled phase velocities are robust and fit well with the observations with small residuals (Fig. 6c, 7, and 8c). We, therefore, employ those modelled phase velocities to obtain final models as described in Section 5 and Appendix B.

APPENDIX B: MODEL REGULARIZATION
We assume layered structures to obtain final S-wave velocity models. The details of layer definitions and the scaling between S-wave velocity, P-wave velocity and density are as same as those used during phase-velocity measurements (Section 4). There are only two differences: the cost function and the vertical smoothing parameter, ε, which control and reduce the model uncertainties. The method follows Takeo et al. (2013) as described below.
We define cost function as where N is the number of phase-velocity measurements used in the inversion, and V i S is S-wave velocity in the ith layer. The measurements are taken by equal intervals in log frequency. The second term is only summed for layers in the crust. We obtained 100 models corresponding to 100 phase-velocity curves for each area, from which we calculated the average and standard deviation of S-wave velocities for each layer to define model uncertainties.