Spatiotemporal Variations and Postseismic Relaxation Process Around Mt. Fuji, Japan, During and After the 2011 Tohoku‐Oki Earthquake

To monitor the volcanoes at a high spatiotemporal resolution, we introduce the singular value decomposition‐based Wiener filter and the three‐component waveforms in ambient noise velocity monitoring. The continuous ambient noise data from 63 stations around Mt. Fuji and Mt. Hakone, Japan, during the January‐September 2011 were analyzed to estimate the seismic velocity variations at a 1‐day temporal resolution, allowing us to distinguish the velocity drops caused by the 2011 Mw 9.0 Tohoku‐oki and the Mw 6.0 East Shizuoka earthquake. The velocity drop during the Tohoku‐oki earthquake was large in volcanic areas and was larger around Mt. Hakone than Mt. Fuji. This difference is possibly due to the existence of fluid‐ and gas‐rich zones at shallower depths and a higher crack density around Mt. Hakone. In addition, the velocity drop at Mt. Fuji during the Tohoku‐oki and the East Shizuoka earthquake was the same level, despite larger static stress changes beneath Mt. Fuji during the East Shizuoka earthquake. We interpret this inconsistency between the velocity drops and static stress changes to arise from incomplete recovery of the generated cracks during the Tohoku‐oki earthquake when the East Shizuoka earthquake occurred. This study also investigates the spatial variations in recovery speed and recovery amount, finding slow recovery speeds in the volcanic areas and fault areas, possibly due to larger crack densities in the crust. Furthermore, we observe the lowest velocity recovery amount in the volcanic areas, which is likely attributed to the maintained increase in pore pressure due to the volcanic gas bubbles.


Introduction
Accurate estimations of future volcanic activity are important for disaster mitigation.Continuous volcanic monitoring at a high temporal resolution is a key requirement for volcano warning systems.Several monitoring tools, such as geodetic and seismic stations, are regularly used to continuously monitor active volcanoes worldwide.Ambient noise monitoring is a relatively new monitoring approach that has been widely used for crustal-scale monitoring (Brenguier et al., 2014;Hobiger et al., 2016;Ikeda & Tsuji, 2018;Meier et al., 2010), whereby the spatiotemporal seismic velocity variations are estimated without active seismic sources.If this monitoring system can provide deep crustal constraints in active volcanic areas (e.g., the magmatic area) at a high temporal resolution, it can provide additional information for the continuous monitoring of volcanoes.
Ambient noise cross-correlation provides an estimate of the Green's function, which is the impulse response between a given station pair (Campillo & Paul, 2003;Sabra et al., 2005;Shapiro & Campillo, 2004).The seismic waves that are extracted via ambient noise cross-correlation are generally considered to be dominant surface waves.The surface-wave velocity exhibits a frequency dependence based on its penetration depth.Therefore, the physical changes at specific depths can be studied by estimating the velocity changes from the extracted surface waves that are obtained for different time periods.We can use the information from the cross-correlation functions (CCFs) of coda waves (i.e., the scattered waves in the subsurface due to crustal heterogeneities) to obtain stable estimates of the velocity changes (e.g., Meier et al., 2010;Sens-Schönfelder & Wegler, 2006).The ambient noise-derived monitoring data can provide continuous information that is relevant to volcanic activity because the observed velocity changes are related to the crustal stress state, pore pressure changes, crack generation, and/or the presence of volcanic and hydrothermal fluids (Andajani et al., 2020(Andajani et al., , 2022;;Brenguier et al., 2014;Wang et al., 2017).However, most studies of ambient noise monitoring have used a longer period of data (i.e., longer data window) for the velocity monitoring, such that the monitoring results generally possess a lower temporal resolution (e.g., several to 10 days; Hutapea et al., 2020).This lower temporal resolution makes it difficult to detect sudden events, such as rapid volcanic eruptions.A continuous ambient noise monitoring system with a high temporal resolution is therefore required to better detect such events.Volcanic complexes in seismogenic regions are commonly active following large earthquakes (Manga & Brodsky, 2006;Nishimura, 2017); therefore, it is important to evaluate the effects of earthquakes on this volcanic activity based on monitoring with high temporal resolution.For example, Karpinsky volcano erupted the day after the 4 November 1952 M w 9.0 Kamchatka earthquake.Furthermore, the 1701 Hoei earthquake (M w 8.7) in the Nankai subduction zone triggered one of the largest historical eruptions of Mt.Fuji (e.g., Chesley et al., 2012).The 2011 M w 9.0 Tohoku-oki earthquake (hereafter, the Tohoku-oki earthquake) triggered seismicity globally (Gonzalez-Huizar et al., 2012), including in volcanic areas (Yukutake et al., 2011).The volcanic areas around Mt. Fuji and Mt.Hakone, which are within the collision zone of the Eurasian, North American, and Philippine Sea plates, experienced a 0.1%-1.0%decrease in velocity (Brenguier et al., 2014;Yukutake et al., 2013) after the Tohoku-oki earthquake, and subsequent seismicity was triggered around Mt. Hakone.Furthermore, the 15 March 2011 M w 6.0 East Shizuoka earthquake occurred 4 days after the Tohoku-oki earthquake at the southern foot of Mt.Fuji, and a subsequent increase in seismic activity was observed around Mt. Fuji for several weeks (Aizawa et al., 2016).However, these previous studies had a low temporal resolution, with average velocity changes captured over a 5day window, and were unable to separate the velocity changes caused by the Tohoku-oki earthquake from those caused by the East Shizuoka earthquake.Furthermore, the detailed spatial distribution of the velocity changes around Mt. Fuji and Mt.Hakone were not clarified due to the low spatial resolution of these previous studies.
Here we improve the temporal resolution of the monitoring results using three-component ambient noise records, the singular value decomposition-based Wiener filter (Moreau et al., 2017), and a denser observation data set.This approach allows us to monitor velocity changes with a high temporal resolution of one day.This increased temporal resolution means that the temporal velocity changes associated with earthquakes can now be discussed in more detail and the earthquake mechanisms are more discernible.Furthermore, the high temporal resolution of this approach can provide important information for volcano warning systems.

Data Preparation and Cross-Correlation
We used the continuous seismic waveforms that were recorded during the January-September 2011 period at 63 three-component seismic stations from Hi-net (NIED, 2019a), V-net (NIED, 2019b), and volcano seismometers of Japan Meteorological Agency and Hot Springs Research Institute of Kanagawa around Mt. Fuji and Mt.Hakone (Figure 1) to monitor the spatiotemporal velocity variations around Mt. Fuji and Mt.Hakone during the Tohoku-oki earthquake.
The estimated velocity changes via ambient noise seismic interferometry are strongly affected by earthquakes because earthquakes happen irregularly, and their amplitudes are much larger than those of ambient noise.Many earthquakes occurred around Mt. Fuji and Mt.Hakone after 11 March 2011.The continuous waveform data were therefore normalized in the time domain, following Bensen et al. (2007), to reduce the influence of these earthquakes and obtain stable monitoring results from the ambient noise analysis.The raw waveform data were divided by the moving average of the amplitudes of the 0.1-0.9-Hzbandpass-filtered data to reduce the earthquake amplitudes, and these normalized data were then filtered in the 0.1-0.9-Hzfrequency band.Furthermore, the seismometer orientations at the stations were corrected and transformed from EW and NS components to radial and transverse components, respectively.
The daily waveform data were then divided into 95 30-min segments with 50% overlap.We defined the threshold to be five times of the root-mean-square amplitudes for 1-day seismogram.If the amplitude of the waveform is above the threshold, we discarded the segment.A Hamming window and spectral whitening were applied to each segment, and the CCF was calculated using the three-component waveforms (Figure 2).Cross-correlation with spectral whitening is performed in the frequency domain as follows: where CCF(ω) is a CCF, F(ω) and G(ω) are the waveforms of the same time segment at different stations in the frequency domain, and G* is the complex conjugate of G.We obtained an 18-component CCF that consisted of causal and acausal 3 × 3 cross-correlation tensors.The causal CCF is the seismogram with an impulse source at X and receiver at Y, where X is the station recording F(ω) and Y is the station recording G(ω).The acausal CCF is the seismogram with an impulse source at Y and receiver at X.The CCF was stacked for one day to construct the 1-day waveforms, and also for the entire observation period to construct a reference waveform.

SVD-Based Wiener Filter for Denoising
The SVD-based Wiener filter (SVDWF) is the denoising approach introduced by Moreau et al. ( 2017) to improve the quality of noisy CCFs.The SVDWF  improves the SNR by first decomposing the CCF into a singular value that is associated with a level of energy and then applying the Wiener filter, which minimizes the incoherent signals via a statistical approach.The daily CCF matrix C, with dimension M × N, can be transformed into three matrices, U, S, and V, via singular value decomposition: where U is a unitary matrix of dimension M × M that contains the left singular vectors, S is a diagonal matrix of dimension M × N that contains the singular values, V is a unitary matrix of dimension N × N that contains right singular vectors, and M and N are the date and time dimensions, respectively.We can define to rewrite Equation 2 as follows: where s i is the ith singular value, and u i and v i are the ith left and right singular vectors, respectively.Each singular value represents a level of energy.The large singular values contain energetic signals, and the small ones contain less energetic signals.In other words, the coherent signal is dominant for the subspaces with the first several large singular values, and the incoherent signal and noise are dominant for the other subspaces with small singular values (Figure S1 in Supporting Information S1).We can improve the quality of the CCF by rejecting the subspaces with small singular values.We used 50 singular values to reconstruct the CCFs (Figure S2 in Supporting Information S1).
A Wiener filter is then applied to each ith singular vector s i u i v i T to minimize the incoherent signal and noise.CCFs are also reconstructed using the remaining smaller singular values, and the same Wiener filter is applied to the reconstructed CCFs to denoise the overall incoherent signal and noise.
The window size of the Wiener filter affects the velocity change estimation.Increasing the window size along the date axis suppresses the fluctuations in the estimated velocity changes, although at the expense of weakening the velocity drop and making it more gradual.Increasing the window size along the time axis removes highfrequency incoherent signal and noise, although at the expense of also removing some coherent signal and decreasing the correlation coefficient (CC).We tested several window sizes and selected [1 12], which corresponds to [1 day 1.2 s], to retain a high temporal resolution and remove incoherent noise.Figure 3 shows the raw and denoised CCF between stations FJSV and FJHV (Figure 1).We applied the SVDWF to the current CCFs, but not to the reference CCF, to preserve the coherent signals in the reference CCF.

Stretching Interpolation for Velocity Change Estimation
The stretching interpolation method (Meier et al., 2010) estimates slight velocity changes by stretching and shrinking the current waveform and then comparing it with the reference waveform.We used the denoised 1-day stacked CCF as the current waveform and the 242-day (entire observation period) stacked CCF as the reference waveform to estimate the velocity changes between the current and reference waveforms.We focused on the coda-wave part of the CCF, which is composed of the scattered waves after the direct surface wave.Coda waves propagate along a longer path than direct waves, thereby making it easier to estimate slight velocity changes.The stretched current waveform is defined as follows: where ε is the stretching parameter and f cur (t) is the current waveform.We then searched for ε that maximizes the CC cc(ε): Journal of Geophysical Research: Solid Earth 10.1029/2023JB027978 where f ref (t) is the reference waveform.The seismic velocity change dv/v is simply described using ε max , which maximizes cc(ε) in Equation 5: We defined the coda window as a 100-s window with time zero at 15 s after the surface-wave arrival time.The surface-wave arrival time was determined by dividing the surface-wave velocity by the distance between the observation station pairs.The surface-wave velocity was obtained from the slope of the best-fit line to the maximum absolute value of the amplitude across the observed arrival times.The solid and dashed red lines in Figure 3c represent the estimated surface-wave arrival times and start of the coda window, respectively; the estimated surface wave velocity is 2.98 km/s.We calculated the root mean square error of dv/v using an equation derived by Weaver et al. (2011) to evaluate the precision of the velocity changes.Equation 7gives the root mean square error of the estimated velocity change: where cc max is the maximum CC via stretching, T is the inverse of the frequency bandwidth, ω c is the center frequency, and the coda window extends from time t 1 to time t 2 .The velocity change was estimated for each of the 18 components of the CCF.

Average Velocity Change of Multi-Component
The averaged velocity change was estimated by combining the 18 components of the velocity changes, following Hobiger et al. (2012), to reduce the noise in the velocity changes.The velocity change Δ(t) and cc(t) are: and: respectively, where Δ k (t) and cc k (t) are the velocity change and CC of component pair k = 1,2,…,18 at time t.We defined the error of the averaged velocity change σ c as follows: where σ k is the root mean square error of dv/v in component pair k, and Δ(t) is the velocity change between the observation stations (i.e., the velocity change along the ray path between the stations).

Mapping of Velocity Change
To visualize the spatial distribution of the velocity changes, we calculated the velocity change at the observation stations based on the method in Hobiger et al. (2012).The basic idea of the Hobiger et al. (2012) approach can be expressed as follows: where ∆ obs,i is the velocity change for the station pair, and ∆ station,i 1 and ∆ station,i 2 are the velocity changes at the two stations.This equation can be formulated for all of the observation station pairs that are within 60 km of each other and incorporated into a single matrix to obtain the following equation: Journal of Geophysical Research: Solid Earth where n and m are the numbers of station pairs and stations, respectively.We can rewrite Equation 12as The spatial distribution of the velocity changes was obtained by spatially interpolating the velocity changes at each observation point.The velocity changes on the first day of the observation period were set to zero to determine the relative velocity changes.

Curve Fitting for the Observed Velocity Change
We performed curve fitting to the velocity changes between station pairs to evaluate the detected velocity changes.Seismic velocity changes can be influenced by external perturbations, such as the tides and precipitation (e.g., Sens-Schönfelder & Wegler, 2006), which also include seasonal variations related to the environment (Wang et al., 2017).However, we did not adopt the seasonal variation term in the fitting model because Wang et al. (2017) demonstrated that the velocity changes around the Fuji-Hakone area exhibit a weak seasonality.We fitted the velocity changes based on Hobiger et al. (2014Hobiger et al. ( , 2016) ) and Ikeda and Tsuji (2018).
The fitting curve for the observed velocity change f(t) is formulated as follows: where t is the day; A, B, C, and D are fitting parameters; H(t) is the Heaviside function; and t eq is the timing of the Tohoku-oki earthquake.A is a constant offset that represents the average velocity before a given earthquake, B is the non-recovered coseismic velocity change, and C is the recovered coseismic velocity change, which is exponentially scaled using time constant D. We chose the best-fit B, C, and D parameters as those that minimized the square of the standard deviation in Equation 14.We then plotted the parameters in the middle of each station pair and compared those within 60 km of each other to investigate the spatial distribution of each parameter.

Results
We estimated the seismic velocity changes between all of the station pairs that were within 60 km of each other during the January-September 2011 period, which included the Tohoku-oki and East Shizuoka earthquakes (Figures 4 and 5).The velocity changes that were estimated using only single-component data are too noisy to discern the daily variations (gray lines in Figures 4a and 4b), whereas the averaged velocity changes from the 18component CCFs are stable (red lines in Figures 4a and 4b).There was a significant reduction in the noise levels of the averaged velocity changes, such that we could clearly distinguish the velocity drop due to the Tohoku-oki and East Shizuoka earthquakes.Even though the CC is less than about 0.6 for the two presented station pairs, the observed velocity drop can be clearly attributed to the two earthquakes, because the velocity drop is larger than the error σ (about 0.01%; Figures 4c and 4d).The velocity changes for the FJHV-FJSV station pair, which cross the summit of Mt.Fuji, highlight seismic reductions after both the Tohoku-oki (∼0.2%) and East Shizuoka (∼0.2%) earthquakes (right panel of Figure 4a).The velocity changes for the ASGH-MNZH station pair, which are in the Mt.Hakone area, highlight a strong velocity reduction (∼0.4%) during the Tohoku-oki earthquake (right panel of Figure 4b); however, no velocity drop due to the East Shizuoka earthquake was observed, despite the epicenter being only ∼30 km from Mt. Hakone (Figure 1).The CC was stable during the observation period for the FJHV-FJSV station pair (Figure 4c), and the CC decreased during the time period between the Tohoku-oki and East Shizuoka earthquakes for the ASGH-MNZH station pair (Figure 4d).
We estimated the local velocity changes at each station using the least-squares method (Equation 13) to investigate the spatial distribution of the velocity changes caused by the Tohoku-oki and East Shizuoka earthquakes.We then mapped the velocity changes by spatially interpolating the station-dependent velocity changes (Figure 5).The velocity changes during the 10-13 March 2011 period (i.e., before and after the Tohoku-oki earthquake) were mapped by setting the velocity changes on 10 March to zero to only extract the influence from the Tohoku-oki earthquake (Figure 5a).A distinct coseismic velocity drop was observed in the area surrounding Mt.Fuji and Mt.Hakone, and the largest velocity drop was observed around Mt. Hakone.The velocity changes during the 13-15 March 2011 period (i.e., between the Tohoku-oki and East Shizuoka earthquakes) revealed minimal velocity changes across the study region (Figure 5b).The velocity changes during the 15-16 March 2011 period (i.e., before and after the East Shizuoka earthquake) showed that the coseismic velocity drop only occurred around the epicenter of the East Shizuoka earthquake and possessed an amplitude that was considerably smaller than that of the Tohoku-oki earthquake (Figure 5c).The East Shizuoka earthquake did not induce seismic velocity changes around Mt. Hakone, which is only 30 km from the summit of Mt.Fuji.We note that the surface-wave velocity in our frequency range (0.1-0.9 Hz) is sensitive to the regional structure from the surface to about 10 km depth (Brenguier et al., 2014;Nimiya et al., 2017).We fitted the velocity changes using the curve-fitting approach in Equation 14to investigate the characteristics of the velocity changes in a more quantitative sense.Figure 6 shows the seismic velocity changes and their best-fit curves for the FJHV-FJSV station pair, which crosses the summit of Mt.Fuji, and the ASGH-MNZH station pair, which crosses near the summit of Mt.Hakone.The recovered coseismic velocity drop for ASGH-MNZH is notably larger than that for FJHV-FJSV.
We also plotted the spatial distributions of the fitting parameters in Equation 14 (Figure 7).The dot sizes in Figure 7 represent the residual between the velocity changes and models, whereby a large dot implies small residuals and vice versa.B, which is the non-recovered coseismic velocity drop, is large in the volcanic areas, especially Mt.Fuji and Mt.Hakone.C, which is the recovered velocity drop, is also large in the volcanic areas, although it exhibits greater spatial variability than B. D, which is the recovery time constant, is relatively large in the volcanic areas and notably large to the east of Mt.Fuji and across the Izu collision zone (dashed and solid black circles, respectively, in Figure 7c).Uncertainties of these velocity fitting parameters are provided in Supporting Information S1.In Figure S3 in Supporting Information S1, we also showed the relationship between these fitting parameters.

Interpretation and Discussion
We do not observe large velocity changes in the study region before the Tohoku-oki earthquake (e.g., Figures 4a  and 4b), whereas significant velocity drops are detected around Mt. Fuji and Mt.Hakone, with maximum velocity drops of 0.25% and 0.5%, respectively, immediately after the Tohoku-oki earthquake (Figures 5a and 6).Brenguier et al. (2014) reported a similar result, whereby the area that experienced a large velocity decrease correlated with the main volcanic zone on the eastern part of Honshu Island and determined that volcanic areas are sensitive to stress perturbations due to the presence of pressurized hydrothermal and/or magmatic fluids at depth.These large velocity drops around Mt. Fuji and Mt.Hakone could be caused by the generation and/or opening of cracks in the crust due to pressurized volcanic fluids.Furthermore, the velocity drop around Mt. Hakone is larger than that around Mt. Fuji.This observation can be explained by two interpretations.One interpretation is that the depth difference between the fluid-and gas-rich regions could induce a variable velocity drop.Yukutake et al. (2015) reported that low V p /V s anomaly exist at 3-10 km depth beneath Mt.Hakone, which could be crack- filled fluids or gases that have been supplied from the deep magma chamber of Mt.Hakone.However, such lowvelocity anomalies exist at 7-17 km depth beneath Mt.Fuji (Nakamichi et al., 2007).Therefore, pressurized magmatic fluids exist at 3-10 km depth beneath Mt.Hakone and 7-17 km depth beneath Mt.Fuji.These depth differences in pressurization would indicate that the velocity changes beneath Mt.Hakone are more likely to occur at shallower depths than those beneath Mt.Fuji.As we noted, we performed our analysis in the 0.1-0.9-Hzfrequency range, such that the observed surface-wave velocity penetrates to 10 km depth at most (Brenguier et al., 2014;Nimiya et al., 2017).Our depth constraints should therefore yield a larger velocity reduction around Mt. Hakone than around Mt. Fuji, as observed.The second interpretation is that the difference in velocity drops can be due to spatial variations in the density of crustal fractures.Mt.Hakone is in the Izu collision zone, which is where the Izu-Bonin arc is colliding with the Honshu arc.Furthermore, a large change in the velocity field has been observed around and to the south of Mt.Hakone and has been considered to represent a shear deformation zone (Doke et al., 2021).Therefore, the larger velocity drop around Mt. Hakone seems to be related to the heavily deformed fault zones associated with the Izu collision and shear deformation zone, where cracks are well developed.This second interpretation is consistent with previous research, whereby a decrease in seismic velocity can be attributed to an increase in crack density (O'Connell & Budiansky, 1974).
We could not detect velocity changes during the 13-15 March period, which spans the time interval after the Tohoku-oki earthquake and before the East Shizuoka earthquake (Figure 5b).A local velocity drop of ∼0.2% was observed around Mt. Fuji after the East Shizuoka earthquake, with this velocity drop being as large as that during the Tohoku-oki earthquake (Figure 5c).The seismic velocity generally depends on the static stress conditions (e. g., Christensen & Wang, 1985;Toksöz et al., 1976); however, the static stress changes beneath Mt.Fuji during the East Shizuoka earthquake were 10-100 times larger than those during the Tohoku-oki earthquake (Fujita et al., 2013).We could explain this inconsistency between the velocity and static stress changes from the perspective of the opening and closing of cracks.The seismic velocity generally decreases during earthquakes in response to the opening of cracks caused by ground shaking and an increase in pore pressure before gradually recovering due to the closure of cracks.However, the pore pressure had already increased due to volcanic gas migration during the Tohoku-oki earthquake (Aizawa et al., 2016), thereby making it difficult to open new cracks further during the East Shizuoka earthquake.This result indicates that any already damaged crust could be less susceptible to further damage.
The spatial distributions of the fitting parameter values highlight relationships between the local geological features and the fitting parameters (Figures 7 and 8).D, the recovery time constant, was larger in the volcanic areas compared with other areas, especially to the east of Mt.Fuji and to the south of Mt.Hakone (dashed and solid circles, respectively, in Figure 7c).Saade et al. (2019) observed changes in the polarization anomaly and anisotropy change to the east of Mt.Fuji after the Tohoku-oki earthquake and suggested that these were caused by changes in crack alignment in response to the reactivation of hydrothermal fluids.To the south of Mt.Hakone, many low-velocity anomalies have been observed between the microplates due to the Izu collision (Nimiya et al., 2020).Therefore, these highly fractured areas likely undergo a slow postseismic recovery.
B, the non-recovered velocity drop, was large around Mt. Fuji and Mt.Hakone, whereas C, the recovered velocity drop, had large variations.Therefore, we calculated the recovery amount (or recovery rate), which is the ratio between C and the total coseismic velocity drop, B + C (Figure 7d).The recovery amount was large for the station pairs that were far from Mt. Fuji and Mt.Hakone, and small (i.e., significant non-recovered coseismic velocity drop) for the station pairs that were around the volcanoes, suggesting inelastic deformation continued in the volcanic areas.We focused on the two areas where the recovery speed (Figure 7c) was slow: the recovery amount (Figure 7d) was relatively large to the south of Mt.Hakone, whereas it was small to the east of Mt.Fuji (solid and dashed black circles, respectively, in Figure 7d).
The relationship between the fitting parameters and geological characteristics of the study region is summarized in Figure 8.The recovery speed in both the faulted (i.e., the Izu collision zone) and volcanic areas, where cracks are well developed, is slower than that in the normal areas.This occurs because the number of cracks that opened during the earthquake could be larger.The lowest recovery amount are observed in the volcanic areas due to the presence of volcanic fluid or gas.Ground shaking during earthquakes causes the upwelling of gas-rich fluids and/ or gas bubbles (Aizawa et al., 2016;Hill & Prejean, 2005), which increases the pore pressure causing crack opening (Tsuji et al., 2008).The increased pore pressure may then induce seismicity that is triggered for several months (Aizawa et al., 2016;Hill & Prejean, 2005;Yukutake et al., 2013).Therefore, the increase in pore pressure in the volcanic areas can potentially be either prolonged owing to the upwelling of gas-rich fluids and gases or maintained indefinitely due to a continued fluid and gas supply from depth.

Conclusions
We estimated the velocity changes around Mt. Fuji and the surrounding region at a 1-day temporal resolution using the SVDWF and 18-component CCFs.The high temporal resolution of this monitoring approach highlights its ability to contribute to volcano monitoring systems, particularly through providing additional information on the state of the volcanic system.Furthermore, our ability to resolve the temporal variations in the velocity changes (i.e., pore pressure and stress changes) in an extensive area means this approach could contribute to our understanding of how earthquakes can trigger volcanic eruptions (e.g., Andajani et al., 2020).
Our findings are summarized as follows.
1.The high temporal resolution of our monitoring scheme allowed us to separately detect the seismic velocity changes caused by the Tohoku-oki and East Shizuoka earthquakes.Large velocity drops were observed around both volcanic area regions after the Tohoku-oki earthquake, whereas a local velocity drop was only observed around Mt. Fuji after the East Shizuoka earthquake.2. The large velocity drops that occurred around the volcanic areas after the Tohoku-oki earthquake were likely to be caused by the generation of cracks due to pressurized volcanic fluids.The velocity drops were larger around Mt. Hakone than around Mt. Fuji because the fluid-and gas-rich zones are at shallower depths beneath Mt.
Hakone, and the crack density around Mt. Hakone is higher due to the Izu collision.These results indicate that highly fractured areas (volcanic and faulted areas) are sensitive to strong ground shaking, such as the Tohokuoki earthquake, and are more likely to generate cracks.Therefore, the pore pressure in these areas could significantly increase after a large earthquake, with such a pressure and/or stress variation potentially generating local aftershocks.3. The velocity drop due to the East Shizuoka earthquake only occurred around Mt. Fuji and was as large as the velocity drop due to the Tohoku-oki earthquake, even though the static stress change at the magma reservoir of Mt.Fuji was larger during the East Shizuoka earthquake than during the Tohoku-oki earthquake.We interpret this inconsistency between the velocity changes and stress changes to be the result of the generated cracks and stress changes due to the Tohoku-oki earthquake not being fully recovered before the East Shizuoka earthquake occurred.4. The observed slow recoveries in the volcanic and faulted areas may be due to larger crack densities in these areas.The small recovery amounts in the volcanic areas suggesting intensive inelastic deformation are attributed to the prolonged (or indefinitely maintained) increase in pore pressure owing to the upwelling of volcanic gas-rich fluids or gas bubbles from depth.

Figure 1 .
Figure 1.Distribution of the seismic observation stations.The Hi-net and V-net (red circles) stations and volcano seismometers of the Japan Meteorological Agency (white circles) and Hot Springs Research Institute (green circles) are shown.Red triangles denote the summits of Mt.Fuji and Mt.Hakone.The orange star is the epicenter of the East Shizuoka earthquake.Doubleheaded arrows denote the approximate ray paths between the two station pairs that are discussed in the text.

Figure 3 .
Figure 3. (a) Raw and (b) denoised cross-correlations between stations FJSV and FJHV (Figure 1).(c) Shot gather derived from the denoised cross-correlation.Solid and dashed red lines represent denote the surface-wave arrivals and the start of the coda window, respectively.The estimated surface-wave velocity is 2.98 km/s.

Figure 4 .
Figure 4. Single-component (gray) and averaged (red) velocity changes for the (a) FJHV-FJSV and (b) ASGH-MNZH station pairs.The right panels show the enlarged view of the averaged velocity changes around the timing of the Tohoku-oki earthquake.Shaded pink regions denote the standard deviation of the averaged velocity changes.Average correlation coefficient for the (c) FJHV-FJSV and (d) ASGH-MNZH station pairs.Dotted blue lines represent the timings of the Tohokuoki (left) and the East Shizuoka (right) earthquakes.

Figure 5 .
Figure 5. Spatial distributions of the velocity changes during the (a) 10-13 March (i.e., during the Tohoku-oki earthquake), (b) 13-15 March (i.e., between the Tohokuoki and East Shizuoka earthquakes), and (c) 15-16 March (i.e., during the East Shizuoka earthquake) time periods.Black dots represent observation stations, and red triangles denote the summits of Mt.Fuji and Mt.Hakone.The orange star is the epicenter of the East Shizuoka earthquake.

Figure 6 .
Figure 6.Observed velocity changes (black) and best-fit model curves (orange) for the (a) FJSV-FJHV and (b) ASGH-MNZH station pairs.Best-fit parameters (A-D) are provided at the top of each plot.

Figure 7 .
Figure 7. Spatial distributions of (a)-(c) the fitting parameters in Equation 14 and (d) the recovery amount (or recovery rate) C/(B + C).Dashed and solid black circles represent the areas to the east of Mt.Fuji and across the Izu collision zone, respectively.Dot sizes represent the residuals between the velocity changes and models.Dashed and solid black circles in (c) and (d) represent the areas to the east of Mt.Fuji and across the Izu collision zone, respectively, and highlight where D is especially large.

Figure 8 .
Figure 8. Schematic of the temporal variations in velocity changes for the different geological structures across the study region.The qualitative meaning of the fitting parameters B, C and D in Figure 7 are included in the images.