Real-time crustal monitoring system of Japanese Islands based on spatio-temporal seismic velocity variation

Abstract

To continuously monitor crustal behavior associated with earthquakes, magmatic activities and other environmental effects (e.g., tides and rain precipitation), we have developed a continuous monitoring system of seismic velocity of the Japanese Islands. The system includes four main processing procedures to obtain spatio-temporal velocity changes: (1) preparing ambient-noise data; (2) creating virtual seismograms between pairs of seismometer stations by applying seismic interferometry; (3) estimating temporal velocity variations from virtual seismograms by stretching interpolation approach, and (4) mapping spatio-temporal velocity variations. We developed a data-processing scheme that removes unstable stretching interpolation results by using the median absolute deviation technique and a median filter. To map velocity changes with high stability and high temporal resolution during long-term (i.e., longer-term monitoring), we proposed the “sliding reference method”. We also developed evaluation method to select the optimum parameters related to stability and temporal resolution. To reduce computation time for continuous monitoring, we applied parallel computation methods, such as shared memory and hybrid distributed memory parallelization. Using our efficient and stable monitoring system, we succeeded to continuously monitor the spatio-temporal velocity variation of the whole Japanese Islands using ambient-noise data from 767 seismometers. Finally, we developed a web application that displays spatio-temporal velocity changes. In the monitoring results that we open through the website, we identified velocity variation (e.g., pore pressure variation) that could be related to earthquake, aftershock, magmatic activities and environmental effects in a stable manner.

Introduction

Changes in the subsurface environment have been often monitored through repeated seismic surveys using active sources (Chadwick et al. 2010), but the high cost of active-source seismic surveys makes it problematic for continuously monitoring subsurface behavior. The use of a permanent seismic source system such as the accurately controlled routinely operated signal system (ACROSS) (Kumazawa and Takei 1994; Yamaoka et al. 2008) is an effective approach to enhance temporal resolution and source repeatability. In ACROSS, repeatable signals can be continuously generated by rotating eccentric mass. Previous ACROSS-based monitoring succeeded in identifying temporal changes in seismic velocity associated with earthquakes (Ikuta et al. 2002; Ikuta and Yamaoka 2004), volcanic activity (Maeda et al. 2015), ground freezing (Ikeda et al. 2017), and secular velocity change (Tsuji et al. 2018). However, the number of ACROSS units is currently limited, resulting in low spatial resolution in the monitoring.

Seismic interferometry is an alternative method which can estimate variations in seismic velocity using ambient-noise data (e.g., Brenguier et al. 2008; Wegler et al. 2009). This method avoids the cost of active sources by relying on natural earth vibrations as the seismic source (e.g., Brenguier et al. 2008; Tsuji et al. 2016). Seismic interferometry can extract virtual seismograms, simulating the signal recorded at a station from a seismic event occurring at the location of another station, by computing the cross-correlation of ambient-noise data (Wapenaar et al. 2010). To estimate velocity changes, coda waves in the virtual seismograms are used rather than direct waves, because coda waves are more sensitive to velocity changes in geological materials (Meier et al. 2010) and are less affected by variations of noise sources (e.g., Chaves and Schwartz 2016). Several studies have successfully used this method to estimate temporal seismic velocity variations (Brenguier et al. 2008; Ikeda and Tsuji 2018; Meier et al. 2010; Obermann et al. 2013; Wegler et al. 2009). However, mapping of temporal velocity variation for longer-term (i.e., spatiotemporal monitoring; Wang et al. 2017) has not been well established, particularly when unstable velocity variations due to temporal variation of noise source characteristics (Zhan et al. 2013) and scatterers preclude us to map the monitoring results in a stable manner. Therefore, we need to develop a method to obtain stable monitoring results from many seismometer pairs for long-term. Furthermore, we need to design a high-performance computation scheme to continuously calculate velocity variation from huge ambient-noise data. Note that Nagaoka et al. (2010) monitored temporal change in Mt. Asama, Japan using S-coda for multiple natural earthquakes based on seismic interferometry. Their approach has a potential to exclude the biases introduced by temporal variation of noise source characteristics, but the temporal resolution depends on the occurrence of natural earthquakes.

Monitoring of seismic velocity variations has a wide range of applications. Velocity changes can indicate subsurface changes caused by natural activities such as earthquakes, volcanic eruptions, and seasonal cycles. Continuous monitoring may lead to the use of seismic velocity changes in forecasting volcanic eruptions or characterizing aftershock sequences (Brenguier et al. 2008; DeVries et al. 2018; Nimiya et al. 2017). Furthermore, monitoring velocities can be used to track subsurface behaviors due to human activities such as CO₂ geological storage, oil/gas production, and geothermal reservoir managements. If seismic velocities can be used to monitor pore pressure variations in resource reservoirs due to fluid injection or extraction, the information derived from monitoring could help avoid induced earthquakes. Indeed, Taira et al. (2018) estimated the reservoir response due to fluid extraction at the Salton Sea geothermal field.

Permanent seismic networks now provide continuous seismic data in the United States, Japan and elsewhere (Kaneda et al. 2015; Kawaguchi et al. 2015; Okada et al. 2004). Ambient-noise data (i.e., big data) derived from such networks can be used to conduct real-time monitoring of seismic velocity. For example, Taira and Brenguier (2016) constructed a quasi-real-time monitoring system at the Lassen volcanic center in California using ambient-noise data. Recently, seismic velocity variation beneath the seafloor due to strain accumulation process has been monitored by ocean bottom seismometers (Ikeda and Tsuji 2018).

In ambient-noise seismic monitoring, we calculate velocity variations between two seismic traces acquired at different times. One trace, called the current trace, is defined by stacking cross-correlations over a short period (e.g., 10 days), and the other is a reference trace, typically obtained by stacking cross-correlation over a longer period (e.g., 1 year) (Nimiya et al. 2017). Because the velocity change is estimated with respect to the reference trace, the definition of the reference trace is crucial. In previous studies (e.g., Hobiger et al. 2016; Ikeda and Tsuji 2018; Meier et al. 2010; Nimiya et al. 2017; Obermann et al. 2013), reference traces were defined by stacking all the cross-correlations for the respective study periods. However, when compiling monitoring results for longer terms, it is problematic to use such a fixed reference trace. Because the characteristics of cross-correlations vary with time due to earthquakes and magmatic activities, the quality of a reference trace may decrease even when we increase the number of stacking. The results derived from ambient-noise monitoring are influenced by temporal variation of scatters and noise sources (Nimiya et al. 2017). Furthermore, there is a trade-off with current traces in that longer periods of data yield more stable results, but shorter periods offer better temporal resolution.

In this study, we built a system for continuous velocity monitoring using ambient noise and applied the system to the seismometers distributed on the whole Japanese Islands. To evaluate processing parameters in details, we mainly used seismometers in Kyushu Island. In this area, velocity changes associated with the 2016 Kumamoto earthquake (Mw 7.0) and eruption of Aso volcano could be identified (Nimiya et al. 2017). We considered two different methods for calculating reference traces: (1) the absolute reference method (ARM) and (2) the sliding reference method (SRM) that is a new approach for defining reference and current traces. To select an optimum time window of current traces for these two methods, we performed a stability evaluation. Furthermore, we developed an outlier removal method based on median absolute deviation (MAD) (Leys et al. 2013) and a median filter to remove unstable results of temporal velocity changes. Because our continuous monitoring system requires fast computation with proper memory management, we designed the system with a shared memory model and hybrid distributed memory model (cluster parallelization). This system was constructed with Python 3.6.4 (Python Software Foundation 2017) from the Anaconda Distribution to support any Linux distribution. Our parallel computation design used seven server nodes to achieve greater than sixfold increases in computation speed. To display the monitoring results, we installed a web application on a virtual private server to give free access to the updated 1 year monitoring results.

Methods

Our monitoring system retrieved velocity information from ambient-noise data in four processing steps (Fig. 1), described in this section.

Fig. 1

Data processing workflow diagram for estimating spatio-temporal velocity variations

Preparing data

We obtained ambient-noise data recorded by the Hi-net seismic network from the National Research Institute for Earth Science and Disaster Prevention (NIED) server (NIED 2019). In the monitoring of whole Japanese Islands, we used the ambient-noise data of 767 Hi-net stations. We selected the vertical component of ambient-noise data, then we applied de-meaning, de-trending, and instrument correction for the data from each station. We also applied a band-pass filter with the frequency range of 0.1–0.9 Hz, and divided the daily data into 30-min-long segments with 50% overlap. To reduce the data volume, the output from this process was saved in the frequency domain. Data segments with 1 s or more of missing waveforms were rejected. Because earthquakes degrade the quality of ambient noise, we applied one-bit normalization to increase the ambient-noise signal and reduce other noise (e.g., Bensen et al. 2007; Hobiger et al. 2016; Meier et al. 2010; Minato et al. 2012). One-bit normalization has been shown to enhance the stability of monitoring results for our receiver and noise characteristics. Our previous study (Hutapea et al. 2019) demonstrated that one-bit normalization allows us to obtain stable results in Hi-net ambient-noise monitoring with similar processing parameters.

Creating virtual seismograms

By computing daily cross-correlations between the sites of two seismometers, virtual seismograms (i.e., Green’s functions) between the two seismometers can be retrieved (Wapenaar et al. 2010). We computed the power-normalized cross-correlations (cross-coherence) in the frequency domain (e.g., Nakata et al. 2011) from

$$CC_{\text{AB}} \left( f \right) = \frac{{F_{\text{A}} \left( f \right) F_{\text{B}}^{ *} \left( f \right)}}{{\left| {F_{\text{A}} \left( f \right)} \right|\left| {F_{\text{B}} \left( f \right)} \right|}},$$

(1)

where CC_AB is the power-normalized cross-correlation at frequency f between seismometers A and B, F_A and F_B are the Fourier transforms of the seismograms at each seismometer, and $F_{\text{B}}^{ *}$ is the complex conjugate version of $F_{\text{B}}$. The resulting cross-correlation consists of the positive time (causal) and negative time (anticausal) parts, corresponding to the wave propagation from seismometer B to A and from seismometer A to B, respectively. In this study, we use the trace side that the waveform propagates to Mount Aso (positive time) in order to reduce strong time-variant tremor noise from Mount Aso (Kawakatsu et al. 2000; Nimiya et al. 2017; Hendriyana and Tsuji 2019) (Fig. 2). Cross-correlations were calculated for station pairs whose distance apart, as computed by using the Haversine equation (Van Brummelen et al. 2013), was less than 40 km. Whereas, we added several seismometer pairs whose distance is > 40 km where the station pair is sparse. In total, we made 7235 pairs of seismic stations distributed on the Japanese Islands.

Fig. 2

Example of the current trace (blue) and the reference trace (red). The rectangles represent the time windows of the coda part of the seismograms, used to estimate seismic velocity changes by the stretching interpolation method. In order to reduce strong noise from Mt Aso, we use the trace side that the waveform propagates to Mt. Aso (positive time)

Estimating temporal velocity variation

To obtain temporal velocity variations from cross-correlations between two seismometers, we can use either the stretching interpolation method (e.g., Meier et al. 2010; Minato et al. 2012; Nimiya et al. 2017; Yukutake et al. 2016) or the moving-window cross-spectral (MWCS) analysis (e.g., Ratdomopurbo and Poupinet 1995; Clarke et al. 2011) to estimate the velocity changes between a current trace and a reference trace (Fig. 2). Nimiya et al. (2017) compared velocity changes estimated by the stretching method and the MWCS for Hi-net data in Kyushu Island, and demonstrated that when using a single side of cross-correlation to reduce the influence of noise from Mount Aso, the stretching method produced more reliable results. Furthermore, the MWCS method underestimates large changes in seismic velocity due to period skipping in the coda of the cross-correlation (Oliver et al. 2017). Therefore, to obtain reliable velocity variation including large velocity changes, we used the stretching interpolation technique.

Assuming a value of velocity changes, we stretched (shrink and expand) a current trace $f^{\text{cur}}$ and computed the similarity between the stretched current trace $f_{E}^{\text{cur}}$ and the reference trace $f^{\text{ref}}$ as follows:

$$C\left( E \right) = \frac{{\smallint f_{E}^{\text{cur}} \left( t \right)f^{\text{ref}} \left( t \right){\text{d}}t }}{{\left( {\smallint (f_{E}^{\text{cur}} (t ) )^{ 2} {\text{d}}t\smallint (f^{\text{ref}} \left( t \right) )^{ 2} {\text{d}}t} \right)^{1/2} }},$$

(2)

$$f_{E}^{\text{cur}} \left( t \right) = f^{\text{cur}} (t\left( {1 + E)} \right) ,$$

(3)

where E ($= - {\text{d}}V/V$) is the relative velocity change with respect to the reference trace and C(E) is the cross-correlation coefficient between the reference trace and the stretched version of the current trace. We estimated the relative velocity change E that maximizes C(E) by applying a grid search algorithm that searched E in the range − 0.025 < E < 0.025 with a step size of 0.0005. We further applied a ternary search algorithm to increase the resolution of E values. We used 100 s of coda waves to obtain velocity changes by the stretching interpolation. The starting time of coda waves was defined as t_coda = d/v_min, where d is the distance between a station pair and v_min is the minimum apparent velocity between the stations. As v_min, we used a relatively slow velocity of 1 km/s to select coda parts.

Define current and reference traces

The choice of periods of current and reference traces is important to obtain stable velocity changes with high temporal resolution. To estimate daily velocity changes, the time window (stacking period) for the current trace is fixed (N days) and the window is slid forward each day (Fig. 3). We used the latest day of the time window for the current trace as the date of the monitoring result. To define the reference traces, we evaluated two methods: the absolute reference method (ARM) and the sliding reference method (SRM).

Fig. 3

Schematic representation of estimating temporal velocity variation by a ARM and b SRM. N and M are the stacking periods, in days, for the current and reference trace, respectively. In ARM (panel a), the data period for the new current trace is moved by 1 day while the data period of the reference trace is fixed. In SRM (panel b), the data periods for the reference and current traces are both moved by 1 day, then stretching interpolation is applied for the entire length of the reference trace. In SRM, the current trace data always overlap with the reference trace

In the ARM, the reference trace is defined by stacking M days of data fixed at a particular time period (Fig. 3a). We then slide the current trace, with its shorter time window (N days), and estimate the velocity change between reference and current traces by stretching interpolation. This approach (ARM) was often used in previous studies (e.g., Meier et al. 2010; Nimiya et al. 2017).

On the other hand, SRM is a new approach we proposed in this study. In the SRM, the reference trace is also defined by stacking M days of data, but it slides each day (Fig. 3b). The current trace (N days) slides within the period of the reference trace (M days), and the latest days of both traces are always the same. Because velocity variations estimated from different reference traces are difficult to compare in SRM, we calculated temporal velocity changes within the fixed reference period each day. To do this, we applied stretching interpolation M−(N + 1) times for each station pair, using the fixed reference trace and sliding the current trace. Therefore, SRM needs a longer computation time than ARM. In this study, we subtracted the average value of E for the first 30 days of the monitoring period (E₀) from the estimated value of E as follows:

$$E^{\prime} = E - E_{ 0} ,$$

(4)

where E′ is the relative velocity changes with respect to E₀. In our monitoring system, the estimated relative velocity change dV/V′ = –E′ was used in both the ARM and SRM.

Detect and remove outliers

When the value of C(E) is low, the corresponding temporal velocity change is usually unstable. To remove these outliers, we marked all temporal velocity changes for which C(E) was below a threshold value. In this study, we used 0.5 as the threshold because a cross-correlation coefficient smaller than 0.5 indicates a weak relationship between two variables, in this case the reference and current traces (Asuero et al. 2006; Mukaka 2012). We also marked velocity changes if there were multiple peaks in C(E). We removed outliers with the MAD algorithm and a median filter. We defined the range of acceptable data in the following equations:

$${\text{MAD}} = {\text{Median}} _{i} \left( {\left| { s_{i} - {\text{Median}}_{j} \left( {s_{j} } \right)} \right|} \right),$$

(5)

$${\text{Median}}_{j} \left( {s_{j} } \right) - TC *{\text{MAD}} < s_{i} < {\text{Median}}_{j} \left( {s_{j} } \right) + TC *{\text{MAD}},$$

(6)

where s_i is the ith velocity change in the time series of velocity changes at a station pair, and TC is the tolerance coefficient for MAD. For this study we used TC = 3. Velocity changes outside the MAD range were removed. We further applied a median filter with a time window of 3 days. The procedure is diagrammed in Fig. 4.

Fig. 4

Flowchart for detection and removal of unstable temporal velocity variations (outliers). Input data consist of the stretching coefficient C(E) and the temporal velocity changes E for each pair of stations. Given the threshold of C(E) and the tolerance coefficient TC of MAD, we can estimate the range of acceptable velocity variations by using Eqs. 5 and 6. If C(E) is outside that range, the velocity change is rejected. The final step is applying a median filter to obtain stable E values

Map spatio-temporal velocity variations

After estimating the temporal velocity variations between station pairs, we mapped the spatio-temporal velocity changes. To obtain seismic velocity changes at each station, we decided to apply the simple averaging technique (Brenguier et al. 2014). For each station, velocity changes were defined by averaging the velocity changes for all pairs that are separated by 40 km or less including a station (e.g., Brenguier et al. 2014; Nimiya et al. 2017). In order to make 2D velocity map, we make a grid image data using N-dimension linear interpolation (Jones et al. 2019). As we described above, because seismometer pairs are limited due to sparse seismometers in a few areas, we added several seismometer pairs whose distance is > 40 km.

In the monitoring results derived from SRM (Fig. 5 and see Additional file 1: Movie S1 online), we can observe stable spatio-temporal velocity variations of the whole Japanese Islands, and identify velocity variations due to the earthquake, volcanic eruptions and weathering (and tiding), as discussed later. In the next section, we describe how we optimized parameters to define temporal resolution (i.e., how many days of data we need to stack for current traces).

Fig. 5

Spatio-temporal variation of seismic velocity of whole Japanese Islands. This velocity variation was derived by applying SRM to ambient-noise data recorded by 767 Hi-net stations (blue symbols). Each panel shows the latest date within the window of current trace: a 2 February 2016, b 26 April 2016, c 2 August 2016, and d 27 November 2016. Warm colors indicate regions where seismic velocity was decreased. The Kumamoto earthquakes (Mw7) occurred on 16 April, and the Mt Aso eruptions occurred on 7 and 8 October

Optimum parameters for ARM and SRM

There is a trade-off between stability and temporal resolution, depending on how many days of data are stacked for current traces and how reference traces are defined (i.e., ARM or SRM). To find optimum parameters for defining current traces, we assessed the stability of the stretching interpolation result by evaluating different time windows of current traces for ARM and SRM. C(E) for the estimated velocity change (Eq. 2) is commonly used for this stability evaluation (e.g., Budi-Santoso and Lesage 2016; Yukutake et al. 2016). We evaluated stability according to the number of data points for which the correlation coefficients for the estimated velocity changes exceeded the threshold value (0.5 in this study). The input data for this evaluation consisted of the velocity variations before outlier removal.

For our stability evaluation, we used ambient-noise data (vertical component) recorded by Hi-net stations (~ 96 stations) in Kyushu from 2015 to 2016, a period including the 2016 Kumamoto earthquake, and tested different stacking periods to create the current trace. The results (Fig. 6) showed that as expected, increasing the time window yielded more stable monitoring results for the ARM and SRM. In the SRM, the incremental improvement from each added day was less than 5% when we stacked more than 11 days of data to make the current trace, therefore we chose 11 days of stacking as the best compromise between stability and temporal resolution. Likewise, we could choose 17 days of data to make the current trace for the ARM. The chosen time window depends on several parameters, such as the window size for reference traces (here 1 year), seismometer density and station interval, ambient-noise characteristics, lithology, and the frequency range for analysis. Therefore, the stability evaluation (Fig. 6) should be done for any application of our monitoring method.

Fig. 6

The result of stability evaluation for the ARM and SRM in Kyushu. The study period includes the 2016 Kumamoto earthquake. The time windows chosen to define the current trace are 17 and 11 days for the ARM and SRM, respectively

We found that the seismic velocities clearly decreased around the hypocenter of the 2016 Kumamoto earthquake (Mw7) and Mt Aso eruption after the earthquake by using both the ARM (Fig. 7a–e) and the SRM (Fig. 7f–j). We assumed that surface waves dominated the coda wave and that the surface-wave velocity was close to the S-wave velocity. Because we analyzed the frequency range from 0.1 to 0.9 Hz, our results were sensitive to S-wave velocity variations shallower than ~ 10 km depth (see Fig. 3 in Nimiya et al. 2017). A different choice of frequency range should reveal velocity variations for a different depth range.

Fig. 7

Spatio-temporal velocity maps determined by using a–e the ARM and f–j the SRM. The ARM and SRM have different temporal resolutions because the current trace was made with stacking periods of 17 and 11 days, respectively. Both methods identified spatio-temporal velocity changes caused by the Kumamoto earthquake. But, the velocity changes derived from both methods are largely different

The characteristics of velocity variations differed in the results of the ARM and SRM (Fig. 7). Thus, we mapped C(E) derived by both methods to evaluate the stability of the monitoring results (Fig. 8). In the results using the ARM, C(E) suddenly decreased at the time of the earthquake, particularly around the hypocenter and Aso volcano (Fig. 8a, b). This indicates that the characteristics of the coda wave changed after the earthquake, probably the effect of a change in scattering due to the earthquake (Obermann et al. 2013) and change in ambient-noise characteristics. The C(E) value for the ARM significantly decreased in October 2016 (Fig. 8e), which could be related to multiple eruptions occurred at Aso volcano from mid-September 2016 to mid-November 2016. On the other hand, in the SRM results, C(E) remained high even after the earthquake and eruption (Fig. 8f–j) because the sliding reference period included time after their occurrence. These results indicate that the velocity variations derived from the ARM are less accurate when the current trace is affected by dynamic events (volcanic activity and a change in scattering) and also when the current trace date is too far separated from the reference trace. The advantage of the SRM is that the reference trace always reflects geological conditions at times close to the current trace, and this method can detect velocity variations with greater stability (Fig. 8) and higher temporal resolution (Fig. 6). The weakness of the SRM is the limited period of velocity variation results (no more than M − (N + 1) days; 1 year in present case) due to its intensive computation. Given these considerations, we decided to use the SRM for monitoring daily velocity variations, using an 11-day window for the current trace and a 365-day window for the reference trace (Fig. 5).

Fig. 8

Stretching coefficient C(E) maps determined by using a–e the ARM and f–j the SRM. The SRM results show higher C(E) values than the ARM results because in the SRM the current trace always overlaps the reference trace. Panels d and e show that C(E) decreased because the current trace has continued to separate in time from the fixed reference trace and also is affected by multiple eruptions, such that the velocity variation could become unstable. The ARM appears to be useful to detect changes in the scatter wave, whereas the SRM determines velocity variations with high stability and high temporal resolution

In order to analyze the error using the SRM, we calculated temporal velocity variations by using shorter time window for coda waves (50 s) and shifting the starting time of the window every 10 s interval (Meier et al. 2010; Nimiya et al. 2017). Using those velocity changes, we estimated the standard deviation for each station pair. We then estimated the error for each station using the standard deviations between all pairs that included a station within a distance of 40 km as follows:

$$\sigma = \sqrt {\frac{1}{N}\frac{{\mathop \sum \nolimits_{i = 1} n_{i} \sigma_{i}^{2} }}{{\mathop \sum \nolimits_{i = 1} n_{i} }}} ,$$

(7)

where $\sigma_{i}$ and $n_{i}$ are the standard deviation and the number of measurements (maximum is 6) for the $i{\text{th}}$ station pair, respectively, and $N$ is the number of station pairs. When only one measurement data was accepted for a station pair after detection and removal of outlier data, we did not include the pair in the error evaluation (Eq. 7), because we cannot define the standard deviation. In the error map using the conventional approach (without detection and removal of outlier data) (Fig. 9a–e), we locally observed large values of errors (~ 0.2%), which are comparable to velocity variation due to the earthquake (Fig. 7g). On the other hand, by considering outlier data using MAD and a median filter, we did not observe such larger values of errors and most of the errors are smaller than 0.02% (Fig. 9f–j), indicating high stability of our monitoring results.

Fig. 9

Error maps determined by using the SRM. a–e Without MAD and a median filter and f–j with MAD and a median filter. Note that color scales for panels a–e and f–j are different. By considering outlier data using our approach, the errors of velocity variations can be significantly suppressed

Develop automated monitoring system

To realize daily monitoring using SRM, we need to consider an effective system for the intensive computation. Here we describe our continuous monitoring system in terms of its three main tasks: (1) downloading ambient-noise data from the NIED server; (2) processing the ambient-noise data to estimate spatio-temporal velocity variations, and (3) displaying the results in a web application on a virtual private server (Fig. 10). We used CentOS 6.10-x64 and Red Hat 6.4-x64 as our main operating system platform to build and test the system. Python 3.6.4 from the Anaconda Distribution was used to ensure that this monitoring system works on various Linux platforms. We used the NumPy library (Van der Walt et al. 2011) as the input–output data format for processing the ambient-noise data.

Fig. 10

Schematic diagram of the ambient-noise monitoring system. The system comprises three operations: downloading ambient-noise data, processing the data, and displaying the spatio-temporal velocity variations on a website. A crontab module starts the update cycle every day at midnight. Presently we update the monitoring results every week on the website

Update spatio-temporal velocity variation

To update spatio-temporal velocity changes automatically, we download the ambient-noise data of previous day (in win32 format) from the NIED server every midnight. Before starting the processing step, we check whether the number of downloaded files is complete or not. The analysis of the ambient-noise data then begins (Figs. 1 and 10). After obtaining both temporal and spatio-temporal velocity variations, the results are transformed into a tabular data format by NumPy. A map of velocity variations is compiled by linearly interpolating the spatio-temporal velocity variations at each station and then converting the results into a grid data format by using Netcdf4 (Unidata 2015). These files are used as input data to make objects for the web application. The web application is updated by a subroutine that checks the modification time of the input files and restarts the application whenever the input files are updated. A crontab module in the Linux operating system automatically starts all of these processes daily.

Parallel computation design

We adopted a parallel computation architecture with a shared memory model and hybrid distributed memory (or cluster) model (Barney 2018). The shared memory model used multi-core central processing unit (CPU) computation inside a single node of workstations, and the hybrid distributed memory model used multi-core CPU and several nodes of workstations connected through a network. We used Python’s Multiprocessing module (Python Software Foundation 2017) to handle the shared memory model and the combination of Multiprocessing and MPI4PY (Dalcin et al. 2011) modules to handle the hybrid distributed memory model.

The ambient-noise processing system consisted of parent, child, and storage nodes. The storage for the parent and child nodes was connected to the storage node through the network file system. During computation in the hybrid distributed memory model, only the child node (rank ≠ 0) ran the computation job and the input–output data were in the storage node. In terms of Flynn’s taxonomy (Flynn 1972), our monitoring system can be categorized as single instruction stream, multiple data streams (SIMD) because ambient-noise data can be processed independently by every seismometer station or every station pair. Although we could use both the shared memory and hybrid distributed memory models in all processing steps in Fig. 1, we did not use the latter in mapping velocity changes because this processing step is faster than the others. During the computation process, the volume of data loaded caused computation crashes when it exceeded the available RAM. To prevent this problem, we wrote a subroutine to divide the input data into portions smaller than 70% of the free RAM before starting the computation.

To benchmark the computation performance, we used a cluster system with eight server nodes (one parent and seven child nodes) and estimated the velocity variation using the SRM. We recorded the total processing time for a system that used hybrid distributed memory parallelization (cluster parallelization) for preparing data, creating virtual seismograms, and estimating temporal velocity variations. The computation process ran only on the child nodes. In theory, this system would yield a sevenfold increase in processing speed, but because the hardware specification of every node was not uniform (the CPU clocks and memory sizes were not uniform), we could not attain that level of performance. Table 1 lists the computation performance for different numbers of child nodes. Using all seven child nodes led to the fastest computation time (33.3 min), whereas using the fastest single child node (with dual Intel Xeon E5-2680 CPUs and 96 GB RAM) required 202.4 min. The use of parallel computation thus sped up processing by 6.09 times and consumed almost 87% of the cluster nodes resource.

Table 1 Parallel computation performance by using hybrid-distribution memory model (cluster) parallelization

Although the monitoring scheme using the SRM is computationally expensive, the high efficiency of the designed system enables us to update monitoring results daily with different processing parameters (e.g., different stacking windows). Our monitoring system has flexibility in computation environments, thus this system could be used in cloud computing.

Display results on website

To make daily monitoring results available to the public, we developed a web application by using Holoviews, Geoviews, Panes (Stevens et al. 2015), and Bokeh (Bokeh Development Team 2018). For the data structure, we used Pandas to manage the tabular data (McKinney 2010) and X-array (Hoyer and Hamman 2017) to manage the grid data file. The application produces three main products: a map of the spatio-temporal velocity variations derived from the SRM (Fig. 11a), the temporal velocity variations for all station pairs (Fig. 11b), and a graph of spatio-temporal velocity variations for each station (Fig. 11c). We included some interaction tools such as zoom, hover tools, pan, and selecting box (Fig. 11d). The newest update for spatio-temporal velocity variation for Kyushu area is shown on http://geo.mine.kyushu-u.ac.jp/tsuji/monitoring.html. We see room for continued improvement in the design and content of the website.

Fig. 11

Template of the web application. a Spatio-temporal velocity changes in Kyushu; b temporal velocity variations for all station pairs; c spatio-temporal velocity variations for individual stations, and d attributes of the data used to calculate velocity variations. This template of web application design is not fixed, because we will update the design in the future in order to improve the feature and performance

Discussion

We developed a continuous spatio-temporal S-wave velocity monitoring system, based on cross-correlation of ambient-noise data from seismic stations, and we report here on its performance during 2015 and 2016, a time period that included the Kumamoto earthquake (Mw 7.0) and eruption of Aso volcano (Figs. 5 and 7). Because the S-wave velocity reflects the shear modulus of the geological formation, the velocity decrease we detected near the earthquake hypocenter and Aso volcano could have been caused by formation damage, increased pore pressure, or pressurized volcanic fluids (Nimiya et al. 2017). In areas far from the hypocenter, velocity variations could be caused by pore pressure changes due to earthquake vibrations. Increased pore pressure could reduce effective stress as well as the friction on the seismogenic fault (Tsuji et al. 2014), thus our velocity variation could be related to the aftershock sequence.

Moreover, before volcanic eruptions, increased pore pressure associated with the intrusion of magma (e.g., Budi-Santoso and Lesage 2016; Obermann et al. 2013) could decrease seismic velocity. Indeed, we detected a decrease in velocity before the Aso eruption and a velocity increase after the eruption, consistent with these pore pressure changes (Figs. 5 and 7).

In addition, our monitoring result (Fig. 5; Additional file 1: Movie S1) contains signals of other processes that may affect earthquake generation. Using ambient-noise monitoring, we can identify the velocity variation due to rain precipitation and snowfall (Wang et al. 2017), because the overburden or diffusion change the pressure conditions of the crust. The tides and snowfall are known to trigger small earthquakes (e.g., Heki 2003; Ide et al. 2016) through their effect on pressure and stress conditions around seismogenic faults. In Kyushu Island, the influence of snowing can be neglected due to little snowfall. However, we identified the velocity variation caused by rain precipitation and ocean tiding (Wang et al. 2017; Fig. 5; Additional file 1: Movie S1). Therefore, our monitoring system can detect dynamic crustal behaviors associated with earthquake triggering processes on a daily basis that may inform disaster warnings and related applications.

Summary

We summarize our monitoring system as follows:

• 1. We built the continuous monitoring system of the whole Japanese Islands. Geodetical approaches such as GNSS or InSAR are typical to identify dynamic crustal behavior, but our monitoring system uses the wave propagating within the crust.

• 2. The monitoring system we developed includes a module to automatically remove unstable velocity results by using the MAD algorithm and a median filter.

• 3. We compared the ARM and SRM for treatment of reference and current traces for stretching interpolation. Because the SRM generated a reference trace that reflects geological conditions close to the current trace, the SRM produced results with greater stability and higher temporal resolution.

• 4. To find optimum parameters for defining current traces, we assessed the stability of the stretching interpolation results by trying different time windows for current traces. On the basis of cross-correlation coefficient C(E), we found that an 11-day stacking window offered the best compromise between stability and temporal resolution in the SRM.

• 5. Our system used parallel computation on seven server nodes to calculate spatio-temporal velocity variations, achieving a greater than sixfold gain in computation speed. We were thus able to continuously monitor seismic velocity using ambient-noise data.

• 6. We developed a web application to enable public access to spatio-temporal velocity variations detected by the monitoring network. The monitoring results can be seen from the following site: http://geo.mine.kyushu-u.ac.jp/tsuji/monitoring.html.