Bridging the gap between low-frequency and very-low-frequency earthquakes

Slow earthquakes that are observed in the > 1 Hz frequency band are called tectonic tremor or low-frequency earthquakes (LFEs) and those in the 0.01–0.10 Hz band are called very-low-frequency earthquakes (VLFEs). These two phenomena are separated by large microseismic noise at 0.1–1.0 Hz. However, recent observations of the signal in this microseismic frequency band accompanying LFEs and VLFEs in the shallow part of the Nankai subduction zone suggest that LFEs and VLFEs are parts of the same broadband phenomenon, “broadband slow earthquakes”. Here, we report the observation of slow earthquake signals in the microseismic frequency band in the western Shikoku region of the Nankai subduction zone, Japan, by stacking many seismograms relative to the timing of the high-frequency LFE signals. We relocate LFE events detected by the Japan Meteorological Agency, use these LFE waveforms to construct synthetic templates, perform a matched-filter event detection analysis using these templates, stack the seismograms recorded by broadband high-sensitivity accelerometers relative to the timing of the detections, and compare the amplitude of the stacked waveforms at different frequency bands. The stacked waveforms have a continuous signal in the 0.015625 Hz (64 s) to 8 Hz frequency band, and support the idea that LFEs are just a small part of the broadband slow earthquake spectrum, which extends to the VLFE frequency band. Furthermore, the frequency dependency of the maximum amplitudes in this study is similar to that of slow earthquakes in the Cascadia subduction zone, and this is also explained by a Brownian slow earthquake model. However, the frequency dependency is inconsistent with the omega-square model, which is a model for ordinary earthquakes.


Introduction
Slow episodic deformation events occur as marginally observable signals in subduction zones and other tectonically active environments. This phenomenon was first discovered in 2000s, with different names used to denote "events" observed in different frequency bands: low-frequency earthquakes (LFEs; e.g., Katsumata and Kamaya 2003) and tectonic tremors (e.g., Obara 2002) in the > 1 Hz band, very-low-frequency earthquakes (VLFEs; e.g., Ito et al. 2007) in the 0.01-0.10 Hz band, and slowslip events (SSEs; e.g., Hirose et al. 1999) that last for days or weeks. However, the definition of a single event is often ambiguous. For example, Shelly et al. (2007) showed that a tectonic tremor is composed of many swarm-like LFEs, and suggested that these are essentially the same process. Ide et al. (2007) further proposed that each of these events (tremor, LFEs, VLFEs, and SSEs) is part of a unified process and termed this phenomenon a slow earthquake, which may follow scaling relations that differ from those for ordinary earthquakes. However, large ambient noise between discrete observable frequency ranges has precluded a comprehensive understanding of the entire range of slow earthquake phenomena. Kaneko et al. (2018) detected broadband signals across all frequencies (0.01-10 Hz) accompanying tremors and VLFEs in the shallow part of the Nankai subduction zone using the Dense Ocean floor Network System Open Access *Correspondence: k.masuda@eps.s.u-tokyo.ac.jp 1 Department of Earth and Planetary Science, The University of Tokyo, Tokyo, Japan Full list of author information is available at the end of the article for Earthquakes and Tsunamis (DONET1), and demonstrated the existence of signals between the VLFE and tremor frequency bands, where microseismic noise is always dominant. Therefore, VLFE and tremor are both part of an inseparable broadband event. This observation is likely common for many other slow earthquake phenomena. However, the strategy of Kaneko et al. (2018) is not applicable to other regions since the microseismic noise level is usually high and such near-field observations are impossible for most slow earthquakes, especially deep events.
An alternative way to investigate slow earthquake phenomena is a statistical approach, such as the stacking of similar signals to improve the signal-to-noise ratio (SNR). Takeo et al. (2010) and Ide and Yabe (2014) stacked broadband signals relative to tremor timing and demonstrated that VLFE signals are accompanied by tremor signals. This procedure has been successfully applied to constrain the focal mechanisms in many regions (Ide and Yabe 2014;Ide 2016;Maury et al. 2016). Rousset et al. (2017) applied a similar stacking technique to geodetic observations, and successfully detected SSE signals during tremor bursts. These observations confirm the universality of the broadband signals that radiate from slow earthquakes, such that these events may be called "broadband slow earthquakes".
Several theoretical and numerical models of broadband slow earthquakes have been proposed (e.g., Ando et al. 2010;Nakata et al. 2011;Ben-Zion 2012;Gomberg et al. 2016;Chestler and Creager 2017;Hawthorne and Bartlow 2018). Ide (2008) proposed a simple Brownian stochastic process that produced broadband slow earthquake signals, known as the Brownian slow earthquake (BSE) model, and Ide and Yabe (2019) developed a corresponding two-dimensional model. The changes in the characteristic source length are based on a stochastic differential equation in the BSE model. The BSE model can explain the randomness of the seismic moment function, a constant seismic moment rate, the flat spectral amplitude of the seismic moment acceleration over a broad frequency range, and the proportionality between the seismic moment rate and seismic energy rate. While the BSE model explains the first observation of broadband slow earthquake spectra in the Kii-Oki region by Kaneko et al. (2018), there are few other examples based on observations, especially in the microseismic frequency band. Ide (2019) first applied a stacking scheme in the microseismic frequency band using broadband LFE signals that were detected in Cascadia (Bostock et al. 2015), and compared these stacked signals with synthetic spectra from the BSE model and omegasquare model (Aki 1967) for ordinary earthquakes. While the spectra of stacked waveforms were close to the BSE model prediction, the signal amplitude was insufficient for comparisons at very broadband frequencies. In the present study, we expected to obtain clearer signals in applying a similar method to the Nankai subduction zone, where tremor activity is high and a wealth of observational data are available. Here, we construct an LFE catalog via a matched-filter method (Gibbons and Ringdal 2006), which is an effective scheme for detecting LFEs from a continuous signal (e.g., Shelly et al. 2007;Bostock et al. 2012;Chamberlain et al. 2014;Frank 2016) and then stack the broadband records relative to the LFE timing to discuss the very broadband nature of slow earthquakes. We select western Shikoku as a study region because various seismic studies have already been conducted in the region (e.g., Shelly et al. 2006;Ohta and Ide 2017).

Overview of method
We employed the method described by Ohta and Ide (2017) to construct an LFE catalog. We first prepared templates using the LFEs that were detected and located by the Japan Meteorological Agency (JMA) (Katsumata and Kamaya 2003). The LFE catalog events that occurred in the open rectangle area (Fig. 1) were relocated via the network correlation coefficient (NCC) method Ide 2008, 2011), and then used to construct synthetic template waveforms in the studied region, where we assumed a 16 km × 16 km fault plane represented the seismogenic portion of the plate interface (Fig. 1). The depth and orientation of this fault plane were determined by minimizing the depth difference between the fault and relocated LFEs in a least-squares sense. We computed the synthetic template waveforms at each grid point (250-m interval) on the fault plane, and then performed a matched-filter analysis to detect additional LFEs in the continuous records. Finally, we stacked the waveforms of all the detected LFEs to assess the nature of the slow earthquake signals over a very broadband frequency range.

Data
The data consist of continuous records from shortperiod velocity seismometers and high-sensitivity accelerometers (tiltmeters) at nine Hi-net stations [Hi-net is the nationwide Japanese high-sensitivity borehole seismometer network, which is maintained by the National Research Institute for Earth Science and Disaster Resilience (Obara et al. 2005)]. These three-component velocity seismometers measure the up-down, north-south, and east-west ground motions with a pendulum frequency of 1 Hz. The original seismograms, which were recorded at 100 samples per second (sps), were down sampled to 20 sps after applying an anti-aliasing filter.
The records were also bandpass filtered at 2-8 Hz for the hypocenter relocations, synthetic template computation, and matched-filter analysis. The tiltmeters are accelerometers that measure the two horizontal components (north-south and east-west) of ground motion. These original data, which are recorded at 20 sps, were used in this study. Tiltmeters are considered to have a flat response at low frequencies when the records are integrated after applying a high-pass filter, similar to the response observed in broadband seismometers (Tonegawa et al. 2006). These records were only used to produce the final stacked records after the matched-filter analysis. The study period covered 15 years, from April 2004 to March 2019, which is a marked improvement over the previous study of Ohta and Ide (2017), which analyzed data over a 7-year period, from April 2004 to March 2011, for the relocations, and a 3-year period, from January 2005 to December 2008, for the matchedfilter analysis.

Relocation by NCC method
The NCC method determines the time difference, t ij , and relative hypocenter location, x ij , between a reference event i and target event j by maximizing the relocated NCC, NCC reloc , which is defined as follows: where n is the station index, l the component index, u ln i the observed waveform for event i at station n and component l , u ln j the observed waveform for event j at station n and component l , and T x ij the travel time difference determined by x ij when the hypocenter of event i is fixed. We used the JMA2001 one-dimensional velocity structure (Ueno et al. 2002) to calculate the travel times. We determined t ij and x ij for all of the event pairs and eliminated outlier pairs, for which the maximum NCC reloc was smaller than 6 times the standard deviation σ NCC or the difference between x ij and − x ji was larger than 1 km. As an exception to the latter condition, if the maximum NCC reloc was larger than 7σ NCC and the maximum NCC reloc of the opposite pair was smaller than 5σ NCC , the pair was not eliminated. Finally, we calculated (1) The black and red dots represent the LFE epicenters before and after relocation, respectively. The gray square shows the fault plane, which is a 16 km × 16 km square that is centered at the centroid of the relocated LFEs. The triangles indicate Hi-net stations used in this study. The inset map shows the location of the study region in western Shikoku relative to Japan the self-consistent hypocenters for the events which had at least one time difference and relative hypocenter location to another event using a least-squares inversion.

Synthetic template waveform
Synthetic template waveforms were computed for each component and station by stacking the waveforms of the LFE events that were observed at a given station and component: where g ln is the synthetic template waveform at station n and component l , x a grid point on the fault plane, C ln the coefficient of normalization for a given component and station, N e the number of LFE events in the relocated catalog, W template a weighting coefficient that is dependent on the distance between two points, x k the event hypocenter, t k the origin time of the event, u ln the observed waveform after removing the mean and trend, dt the sampling interval (0.05 s since the sampling rate is 20 sps), x n the station location, and T the travel time between x n and x k . C ln is given as follows: W template is given as follows: where q is a coefficient that is assumed to be 0.7 km −1 , which is the empirically estimated value of Ohta and Ide (2017). We stacked the data 1.5 s before and 2.5 s after the S-wave arrival time of each event.

Matched-filter analysis
New LFEs were detected during the matched-filter analysis by calculating the correlation between these template waveforms and the continuous observed waveforms. The matched-filter analysis NCC, NCC MFA , used to detect the LFE events is defined for time t and location x as follows: (2) NCC MFA is calculated in 0.1-s increments throughout the entire continuous waveform dataset. We detect a new event when NCC MFA > 7σ , where σ is the standard deviation of NCC MFA for the entire 15-year study period. This threshold indicates the number of false detections is 10-11 for the entire study period if both the template and noise have independent Gaussian distributions. Although the actual number of false detections is greater than this estimation, it is sufficiently small compared with the total detection. We utilized GPU computing to reduce the massive computational cost of this matched-filter analysis for the entire study period.
We stack the tiltmeter records based on the catalog of newly detected LFEs from the continuous velocity seismograms as follows: where h ln is the stacked waveform for station n and component l , N E is the number of waveforms detected via the matched-filter analysis, v ln is the waveform observed by tiltmeter after removing the mean and linear trend, and W stack is the weighting coefficient, which is the L1 norm of the bandpassed waveforms at lower frequencies. W stack is defined as follows: where BP 20−100s is a bandpass filter at 20-100 s. We used 2000s ( M = 20, 000 ) records that were centered at the S-wave arrival time, and measured the L1 norm for this period, excluding the first and last 100 s. Fig. 1), with 1461 events relocated via the NCC method. The rests fail to be relocated due to the lack of relative hypocenter information. The strike and dip of the fault plane, as estimated via least-squares plane fitting of the relocated event distribution, are N216.83° E and 10.80°, respectively. The centroid of the relocated LFE hypocenters is at 33.601° N, 132.867° E, and 32.0 km depth.

The JMA catalog lists 1958 LFEs in the study area from April 2004 to March 2019 (open rectangle in
The number of LFEs ( x and t pairs) detected by our matched-filter analysis is 1,067,710 during the 15-year study period (Additional file 1). A histogram of the (v(t)), temporal variations in detected event frequency is shown in Fig. 2a, where a clear episodicity in LFE activity is observed, which is in agreement with previous findings (e.g., Obara 2002). The spatial distribution of LFEs is heterogeneous, with the LFEs being concentrated in two major patchy regions in the central and northern parts of the study region (Fig. 2b), which is consistent with the findings of Ohta and Ide (2017). Some of the detected LFEs are close in both space and time. This suggests, for example, one event is detected at several neighbor grid points simultaneously. However, we did not eliminate such detections because the purpose of this study was not to investigate the spatial and temporal LFE patterns, but rather improve the SNR of the overall broadband signal. Figure 3 shows an example of the stacked waveforms, calculated using Eq. (7) and the detected LFE catalog. This example is from the east-west component of the tiltmeter at TBEH. Examples of stacked tiltmeter waveforms for each station-component pair are provided in Additional file 2: Fig. S1. The waveforms in the top two boxes are the stacked waveform, which has been converted to velocity data via zero-phase high-pass filtering above 0.005 Hz (200 s) and integration. The ten traces below are the zero-phase bandpass filtered waveforms at ten frequency (period) bands: 0. 125-0.25, 0.25-0.5, 0.5-1, 1-2, 2-4, 4-8, 8-16, 16-32, 32-64, and 64-128 s. The amplitudes are normalized by the maximum amplitude of the trace in the top two boxes. Ide (2019) explained that it is natural to have large amplitude pulses in the 2-8 Hz band, which is the frequency band used in the matched-filter analysis, although the existence of signal outside of this frequency band is not always guaranteed. However, the stacked waveforms show clear coherent pulses at the timing of detection in all frequency bands from 0.125 to 64 s. The pulses often overcome the microseismic noise, which is usually dominant in the 1-20 s band. Most of the components at all stations possess such a lower signal, with the exception of some noisy frequency bands at some stations (gray panels in Additional file 2: Fig. S1).
The present study clearly shows that the stacked LFE waveforms include continuous signals that range from the LFE and tremor frequency bands (> 1 Hz) to the VLFE frequency band (0.01-0.10 Hz). This is consistent with results for LFE families in the Cascadia subduction zone for periods of < 16 s (Ide 2019). The broader frequency range in the present study is due to the large number of stacked signals, the high sensitivity of the tiltmeters at lower frequencies, and quiet observational conditions at the bottom of the Hi-net boreholes. These observations suggest that the slow deformation radiating these signals is a broadband phenomenon, which is similar to the "broadband slow earthquakes" detected in the shallow part of the Nankai subduction zone (Kaneko et al. 2018).
It is not appropriate to measure the signal amplitudes from the stacked waveform spectra (Additional file 2: Fig. S2) in this study, even though clear signals are observed in many of the frequency bands, as shown in Fig. 3 and Additional file 2: Fig. S1. This problem arises because the noise contribution may be significant in these long waveforms. Here, we focus on the variation The maximum amplitudes of the signals in each frequency band are averaged for each component of each station, except for the results where noise dominated the signal, as shown in Fig. 4. The x-and y-axis error bars represent the frequency band of the applied bandpass filter and an unbiased estimation of the standard deviation, respectively. We set two criteria based on the maximum amplitude and root-mean-squared (RMS) amplitude to select the results for which the signal is coherent and exceeds the noise levels. One criterion is that the maximum amplitude between −4T and 4T must be identical to that between −T and T . T is 4 s in the 0.125-0.25, 0.25-0.5, 0.5-1, 1-2, 2-4, and 4-8 s cases, 8 s in the 8-16 s case, 16 s in the 16-32 s case, 32 s in the 32-64 s case, and 64 s in the 64-128 s case. The other criterion is that the RMS amplitude during the −T to T period must be at least twice as large as that during the −4T to 4T period, excluding the period from −T to T . The definition of T is the same for both criteria.
The maximum amplitude increases with the frequency band of the applied bandpass filter until the 2-4 s band. Its slope, ω, is almost constant, and the deviation from this constant slope can be explained by the error. The results from LFEs in Cascadia (Ide 2019) are also plotted in Fig. 4 to highlight the high degree of consistency between these two studies.

Discussion and conclusion
The same maximum amplitude operation was performed on synthetic waveforms that obey the omega-square model (Aki 1967), where normal earthquakes with a 2 Hz corner frequency were modeled in this study, and synthetic waveforms that obey the BSE model (Fig. 4). The large amplitudes of the BSE model at lower frequencies are largely consistent with the observations of the present study and those from Cascadia, with the exception of a deviation slightly below 1 Hz. This deviation reflects the spectral shape of the template wave that was used to make the stacked wave for the BSE model, as explained in detail by Ide (2019). However, the omega-square model is proportional to ω 2 at frequencies below its corner frequency as described in detail in Additional file 2: Text S1, such that the deviation from the observations is significant at very low frequencies. The slope of this deviation does not change, even though the corner frequency of the omega-square model can be arbitrarily changed. Therefore, it is difficult to explain the spectral characteristics of slow earthquake signals via the superposition of a finite number of ordinary earthquakes.
In summary, we performed a matched-filter analysis using 15 years of continuous short-period seismometer and tiltmeter data following the method of Ohta and Ide (2017), and stacked the detected waveforms. The stacked waveforms have coherent signals in the LFE to VLFE frequency bands, and even in the microseismic frequency band, with these results supporting the idea of "broadband slow earthquakes", which represent broadband phenomena including LFEs and VLFEs. The results of this study can be considered robust because Kaneko et al. (2018) obtained the same results under a specific set of conditions. A comparison of the maximum amplitudes in the time domain for each frequency band shows that the results of this study are consistent with LFEs in Cascadia in that they have a constant slope, ω, which is shallower than that of the omega-square model, ω 2 , and that the slow earthquake model must have a spectrum rich in low-frequency components, as observed in the BSE model.

Additional file 1. List of matched-filter analysis detections.
Additional file 2. Additional Figures S1, S2 and Text S1.   Fig. 4 Average maximum amplitude of the stacked waveforms for each frequency band. The average was taken for each of the analyzed frequency bands, except for the case that noise was dominant. Blue inverted triangles, red diamonds, and yellow squares represent the results of this study, the BSE model (Ide 2019), and the omega-square model (Aki 1967), respectively. Green pentagons and hexagons represent data from Cascadia (Bostock et al. 2015), as presented by Ide (2019)