The ability to detect relatively small quakes depends on the spatial density and capability of seismometers. Taiwan and Japan are both the most famous high-seismicity areas in the world. Dense seismometers evenly distributed throughout the whole area are beneficial for monitoring the earthquake occurrences close to and far away from fault zones (Chang, 2014). Earthquake catalogs retrieved from Taiwan and Japan were obtained from the Central Weather Bureau (CWB), Taiwan and the Japan Meteorological Agency (JMA). To distinguish dependencies from independent seismicity, the earthquake catalogs are declustered. Therefore, the ZMAP software package for MATLAB (Wiemer, 2001) was utilized to remove and/or omit influence from duplicate events, such as aftershocks. The declustering algorithm used in ZMAP is based on the algorithm developed by Reasenberg (Reasenberg, 1985). We classify clusters by using the standard input parameters (proposed in Reasenberg, 1985, and Uhrhammer, 1986) for the declustering algorithm. The aftershock clusters in a small area and in a short period of time do not conform to the Poisson distribution, which requires removing the aftershocks from the earthquake sequence. Therefore, some parameters can be set as follows: the look-ahead time for un-clustered events is 1 d, and the maximum look-ahead time for clustered events is 10 d. The measure of probability to detect the next event in the earthquake sequence is 0.95. The effective minimum magnitude cut-off for the catalog is given by 1.5, and the interaction radius of dependent events is given by 10 km (van Stiphout et al., 2012). Earthquakes with depth > 30 km were eliminated from the declustered catalogs to understand seismicity changes before mainshocks mainly in the crust.

Before the analytical processes in this study, we assumed that earthquakes with relatively small magnitude can be the cracks and potentially related to the far mainshocks based on the large seismogenic areas (Bedford et al., 2020). The minimum magnitudes of completeness M_c are 2.0 and 0.0 that can be determined by the declustered earthquake catalogs in Taiwan and Japan, respectively (also see Figs. S1–S4 in the Supplement). The earthquakes with M≥2 are selected and utilized in this study for fair comparison of the seismicity changes during earthquakes in Taiwan and Japan. We classified the selected earthquakes via their magnitudes into three groups (i.e., $\mathrm{3}\le M<\mathrm{4}$, $\mathrm{4}\le M<\mathrm{5}$ and $\mathrm{5}\le M<\mathrm{6}$). Note that the classified earthquakes in each group are determined to be the break events (i.e., the mainshocks). In contrast, the other selected earthquakes with magnitudes smaller than the minima of the classified magnitude are determined to be the crack events.

We construct a spatiotemporal distribution of the crack events for each break quake. The spatiotemporal distribution from 0 to 400 km away from the epicenter of the break quake during a period of 60 d before and after the break occurrence is constructed to illustrate the relationship between the crack events and the break quake in the spatial and temporal domain. Note that the spatial and temporal resolutions of the grids of the spatiotemporal distribution are 10 km and 1 d, respectively, based on the declustering parameters in the ZMAP software (Wiemer, 2001). We count the crack events in each spatiotemporal grid according to distance away from the epicenter and the differences in time before and after the occurrence of the break quake.

Figure 1Spatiotemporal seismic density distributions in Taiwan and Japan. The seismic densities constructed by using the declustered earthquake catalogs of Taiwan and Japan are shown in the left and right panels, respectively. The seismic density reveals changes in seismicity at distances from the epicenters ranging from 0 to 400 km up to 60 d before and after quakes in a particular magnitude group. The superimposed number in each grid is further normalized for a fair comparison by using the total number of quakes and their areas. Notably, the total number of quakes is shown in the title of each diagram.

Download

The superimposition process, a statistical tool utilized in data analysis, is capable of either detecting periodicities within a time sequence or revealing a correlation between more than two data sequences (Chree, 1913). The process is known as the superposed epoch analysis (Adams et al., 2003; Hocke, 2008). In practice, the superimposition is a process to stack numerous datasets that can migrate unique features for a few datasets and enhance common characteristics for the most datasets. The count in each grid of the spatiotemporal distributions for all the break quakes is superimposed as a total one based on the occurrence time and epicentral distance of the break quakes. The total count of the superimposed distribution in each spatiotemporal grid is normalized to seismic density (count per square kilometer) for comparing to the total number of the break quakes and the related spatial area. Moreover, we compute the average values every distance grid using the seismic densities 60 d before and after the quake. The average values are subtracted from the seismic densities, and the obtained differences are divided by the average values in each distance grid to obtain the normalized variation, clarifying changes of the seismic density in the spatiotemporal domain.

The earthquakes with magnitude ≥ 2 listed in the declustered catalogs of Taiwan from January 1991 to June 2017 are utilized to construct a spatiotemporal distribution of foreshocks and aftershocks corresponding to the quakes with $\mathrm{3}\le M<\mathrm{4}$. We superimposed all the crack events corresponding to the 15 625 quakes ($\mathrm{3}\le M<\mathrm{4}$). The seismic density is more than 1000 times greater in a hot region at a distance of 10 km away from an epicenter (which is generally considered to be the gestation area of foreshocks) than it is in areas located >200 km from the epicenter (Fig. 1a). The sudden increase in seismic density suggests that earthquake-related stress accumulates mainly around the hot region, triggering many foreshocks a few days before the earthquakes with $\mathrm{3}\le M<\mathrm{4}$. This partial agreement of the numerous recent studies reported that the seismicity migrates toward the fault rupture zone within tens of kilometers from epicenters a couple of days before earthquakes (Kato et al., 2012; Kato and Obara, 2014; Liu et al., 2019). Meanwhile, the events mainly occur 0–1 d after the quakes, which is irrelevant to the smaller distribution 0–1 d before the quakes (also see Fig. 1). The seismic density close to epicenters (Fig. 1) suddenly increases before and gradually decreases after the quakes. The irrelevance and the differences of changes rates with epicentral distance smaller than 20 km before and after the quakes reveal that the increase in seismicity before the quakes is not contributed by the seismicity after due to the analytical processes in this study. In addition, these analytical results of the seismic activity are also in agreement with the studies in Lippiello et al. (2012, 2017, 2019) and de Arcangelis et al. (2016) with regards to distinct methods.

Figure 2Changes of the normalized spatiotemporal variations in Taiwan and Japan. The normalized variations corresponding to the seismic density in Taiwan and Japan (in Fig. 1) are shown in the left and right panels, respectively. The colors reveal changes of the normalized variations at distances from the epicenters ranging from 0 to 400 km up to 60 d before and after quakes in a particular magnitude group.

Download

On the other hand, the increase in seismic density is not only always limited within the hot region, but also extends outward to a distance of over 50 km away from the epicenters about 0–40 d leading up to the occurrence of the quakes (Fig. 1a). We further examine the spatiotemporal changes in the seismic density up to the M≥4 quakes utilizing the same superimposition process (Fig. 1b and c). The expansion of the increased seismic density about 0–40 d leading up to the occurrence of the quakes and the sharp increases in seismic density a few days before the quakes can be consistently observed using the M≥4 quakes in Fig. 1b and c. Similar results (i.e., the sharp increases in seismic density a few days before the quakes and areas where the increase in the seismicity density is much larger than that of the hot region) can also be obtained using the earthquake catalogs between 2001 and 2010 from the Japan Meteorological Agency (JMA) in Japan (Fig. 1d–f). Note that the earthquakes that occurred on the northern side of the latitude of 32^∘ N were selected from the Japan catalogs. The selection is based on the fact that the earthquakes occurred in the area monitored by the dense seismometer network and to avoid the double count of events in the Taiwan catalogs. The normalized variations corresponding to seismic density in Fig. 1 are shown in Fig. 2. The radii of the positive normalized variations are approximately 50 km while earthquake magnitude increases from 3 to 6 in Taiwan (Fig. 2a–c). The land area of Taiwan is approximately 250 km by 400 km, which causes underestimation of the seismic density in the spatial domain. In contrast, the positive normalized variations roughly expand along the radii ranging from 50 to 150 km, while earthquake magnitude increases from 3 to 6 in Japan (Fig. 2d–f). However, variations in the lead time mostly range from 40 to 20 d, and relationships between the positive normalized variations and the earthquake magnitude can be found in neither Taiwan nor Japan (Fig. 2).

In short, the expansion of the increase in seismic density becomes mitigation and may no longer impact a place at distances > 200 km away from the epicenters for the earthquakes with a magnitude < 6. The increase in seismicity density before the quakes suggests that the accumulation of the earthquake-related stress in the crust originates from the hot region and gradually extends to an external place before earthquakes occur. The area of this external place is several times that of a fault rupture zone that is concluded based on the sparse seismic arrays of the past. If a quake can excite seismicity changes over a wide area (i.e., over 50 km by 50 km), any crustal vibration related to stress accumulation before earthquakes can be too small to be identified from continuous seismic waveforms at one station. In contrast, crustal vibrations can be a common characteristic of continuous seismic waveforms at most stations around fault zones due to the fact that seismicity changes dominated by earthquake-related stress accumulation are distributed in a wide area.

Seismic waveforms obtained from 33 broadband seismometers operated by the National Center for Research on Earthquake Engineering (NCREE) of Taiwan, within a temporal span of approximately 1 year (from June 2015 to June 2016) are utilized in this study. Note that two seismometers of them are eliminated from following the analytical processes due to long data gaps. The principal component analysis (PCA) method (Jolliffe, 2002) is utilized to retrieve the possible stress-related common signals from continuous seismic waveforms on the vertical component at 31 seismic stations over a wide area and to mitigate local noise simultaneously. Figure 3a shows the energy and the cumulative energy of the principal components derived from the continuous seismic waveforms at the 31 stations. The energy of the first principal component is about 12 %, which is more than 3 times that of the following ones. Thus, we determined the first principal component to be the common signals of the ground vibrations before earthquakes. Figure 3b reveals changes in the common signals during the study period along the time. However, no obvious changes can be observed in the temporal domain.

Figure 3The energy and the first principal component derived from vertical seismic velocity data from the 31 stations. The energy and the cumulative energy of the principal components are shown in (a). Bars denote the energy of each principal component. The blue line shows the variation in the cumulative energy from distinct used principal components. The variations in the first principal component during the period (i.e., from June 2015 to June 2016) are revealed in (b). The red vertical line indicates the occurrence time of the M6.6 Meinong earthquake (on 2 February 2016).

Download

Thus, we sliced the common signals into several time spans using a 5 d moving window with 1 d steps to show time-varying changes. The common signals in each time span are transferred into the frequency domain using the Fourier transform to investigate frequency characteristics of ground vibrations before earthquakes. The amplitudes are normalized using the frequency-dependent average values computed from the amplitude 30 d before and after earthquakes via the temporal division. Here, we take the M6.6 Meinong earthquake (Wen and Chen, 2017; Chen et al., 2020c) as an example to understand the changes of the amplitude of the common signals in the spatiotemporal domain (Fig. 4a). Distinct patterns in the amplitude–frequency distributions can obviously be observed before and after the earthquake at a frequency higher than $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$ Hz (also see Fig. 4e and f). The amplitude at the frequency close to $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$ Hz was obviously enhanced approximately 20–40 d before the earthquake. Hereafter, the enhancements were significantly reduced and reached a relatively small value a few days after the earthquake. Meanwhile, the frequency is close to $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$ Hz approximately 60 d before the earthquake and tends to be high near 10⁻³ Hz a few days before the event (also see Fig. 4e and f). We next superimpose the amplitude based on the occurrence time of the 17 earthquakes with $\mathrm{4}\le M<\mathrm{5}$ and the 109 earthquakes with $\mathrm{3}\le M<\mathrm{4}$ during the 1-year temporal span shown in Fig. 4b and c, respectively. The consistent variations (i.e., the frequency is close to $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$ Hz approximately some days before the quakes and tends to be high near 10⁻³ Hz a few days before the quakes) can be observed in Fig. 4b and c.

Figure 4The amplitude ratio of the superimposed time–frequency–amplitude distribution associated with earthquakes with distinct magnitudes. The superimposed results 120 d before and after quakes with the M6.6 Meinong earthquake, $\mathrm{4}\le M<\mathrm{5}$ and $\mathrm{3}\le M<\mathrm{4}$ are shown in (a)–(c), respectively. The distribution is normalized for comparison by using the average amplitude in each frequency band of 30 d before and after the quakes. The total number of earthquakes in each magnitude group is shown in the title of each diagram. Variations in the amplitude ratios in (a)–(c) at frequencies of about $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}$, $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$ and $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{3}}$ Hz during the same period are shown in (d)–(f), respectively.

Download

Here, we retrieve the ratios at three frequencies of approximately $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}$, $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$ and $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{3}}$ Hz to reveal the relationship between the enhancements and earthquake magnitudes (Fig. 4d–f). For the Meinong earthquake, the enhancements could be identified at the low frequency of approximately $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}$ Hz. The ratios exhibit a relatively large value of ∼1.2 about 90 d earlier than the earthquake (Fig. 4d). The ratios rapidly decrease to a relatively small value of ∼0.5 near 60 d before the earthquake. The enhancements with the maxima reach ∼1.6 about ∼30 d before the earthquake. After the earthquake, the ratios fluctuate and recover with a relatively large value of ∼1.2 about 100 d after the earthquake. Regarding earthquakes with relatively small magnitude, the enhancements at $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}$ Hz are ∼1.2 for the group of $\mathrm{4}\le M<\mathrm{5}$ and ∼1.1 for the group of $\mathrm{3}\le M<\mathrm{4}$ between 30 and 50 d before the earthquake occurrence (Fig. 4d). Similarly, the enhancements at $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$ Hz are ∼1.4 for the Meinong earthquake, ∼1.15 for the group of $\mathrm{4}\le M<\mathrm{5}$ and ∼1.05 for the group of $\mathrm{3}\le M<\mathrm{4}$ between 5 and 30 d before the earthquake occurrence (Fig. 4e). The enhancements at $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{3}}$ Hz are ∼1.15 for the Meinong earthquake, ∼1.15 for the group of $\mathrm{4}\le M<\mathrm{5}$ and ∼1.05 for the group of $\mathrm{3}\le M<\mathrm{4}$ between 2 and 30 d before the earthquake occurrence (Fig. 4f). The ratios at the three frequencies in Fig. 4d–f suggest that the amplitude ratios of the enhancements and earthquake magnitudes generally show a proportional relationship. However, the ratios at $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{3}}$ Hz with a relatively large value of ∼1.6 can be observed during the period of 60–45 d before the Meinong earthquake due to unknown disturbances (Fig. 4f).

The findings suggest that the common-mode ground vibrations exist in a wide area before earthquakes due to the signals being retrieved from most stations distributed around all of Taiwan through the PCA method. In short, the common-mode vibrations are very difficult to be identified from the time series data but become significant in the frequency domain. If the expansion of the seismogeneric areas and the existence of the common-mode ground vibrations are true, the next step is to determine the potential mechanism hidden behind this nature.

Walczak et al. (2017) repeatedly observed stressed rocks exciting long-period vibrations during rock mechanics experiments. Leissa (1969) reported that the resonance frequency of an object is proportional to its Young modulus and exhibits an inverse relationship to its mass. Based on the crust, the outermost of the Earth, being lamellar, we assume that the earthquake-related stress accumulates in the volume of a square sheet with a width of 100 km, which is determined by using a distance of 50 km away from an earthquake due to the significant increase in the seismic density (Figs. 1 and 2). The resonance frequency near $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{4}}$ Hz (Fig. 4) can be derived from the square sheet once the thickness of the volume is ranged between 500 and 1000 m (Fig. S5). Although we do not fully understand the causal mechanism of the thickness, the agreement with the spatiotemporal domain of the relatively small quakes from the earthquake catalogs, the superimposition results of continuous seismic waveforms and the resonance frequency models suggests that the phenomenon of variable frequency may exist tens of days before earthquake occurrence and can be retrieved by broadband seismometers.

In this study, we determined the seismogenic areas using the relatively small earthquakes in the spatiotemporal distribution and found that the areas are significantly larger than the fault rupture zone (Figs. 1 and 2). Meanwhile, the ground vibrations can exhibit frequency-dependent characteristics at about 10⁻⁴ Hz (Fig. 4) that could relate to the large seismogenic areas due to the resonance model (Fig. S5). If these are true, the seismo-TEC (total electron content) anomalies in the ionosphere, which is generally observed in a large-scale area with more than 10 000 km² (Liu et al., 2009), have a high potential to be driven by upward propagation of acoustic waves before earthquakes (Molchanov et al., 1998, 2011; Korepanov et al., 2009; Hayakawa et al., 2010, 2011; Sun et al., 2011; Oyama et al., 2016). The existence of the ground vibrations can generate the acoustic-gravity waves that have been reported (Liu et al., 2016, 2017). However, the acoustic-gravity waves in a period of <300 s are difficult to propagate upward into the atmosphere and the ionosphere (Yeh and Liu, 1974; Azeem et al., 2018). The wide seismogenic areas observed in this study can contribute the larger-scale ground vibrations at approximately $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$–10⁻³ Hz that cover the frequency channel ($<\mathrm{1}/\mathrm{300}$ Hz) for the acoustic-gravity waves propagating into the atmosphere and changing the TEC in the ionosphere. Meanwhile, the seismo-atmospheric and the seismo-ionospheric anomalies in a large-scale area can also be supported by the acoustic-gravity waves due to the wide seismogenic areas. While partial aforementioned relationships cannot be quickly proven, the ground vibrations at a low frequency ($<\mathrm{1}/\mathrm{300}$ Hz) in a wide area assist our understanding of the essence of the seismo-anomalies in the atmosphere and the ionosphere.

The earthquake catalogs of Taiwan and Japan were obtained from the Central Weather Bureau (https://www.cwb.gov.tw/, last access: December 2020) and the Japan Meteorological Agency (JMA; https://www.data.jma.go.jp/svd/eqev/data/bulletin/index_e.html last access: December 2020) (JMA, 2020), respectively. Seismic waveform data in Taiwan were provided by the Seismic Array of NCREE in Taiwan (SANTA; http://santa.ncree.org/, last access: December 2020) (NCREE, 2020). The downsampled seismic waveforms with the temporal interval of 10 s can be utilized to reproduce the analytical results in this study through the MATLAB software that can be download at https://doi.org/10.5061/dryad.1jwstqjqq (Chen, 2019).

The supplement related to this article is available online at: https://doi.org/10.5194/nhess-20-3333-2020-supplement.

YYS contributed to discussion and revision. SW contributed to discussion and revision. PH contributed to data collection. LCL contributed to discussion and revision. HZY contributed to discussion. XZ contributed to discussion. YG contributed to discussion. CCT contributed to discussion and revision. CHL contributed to discussion and revision. JYL contributed to discussion and revision.

The authors declare that they have no conflict of interest.

This research was funded by the Joint Funds of the National Natural Science Foundation of China (Grant no. U2039205), the National Natural Science Foundation of China (Grant nos. 41474038 and 41774048), the Spark Program of Earthquake Science of China (Grant No. xh17045), the Ministry of Science and Technology of Taiwan (Grant nos. MOST 106-2116-M-194-016- and MOST 106-2628-M-008-002), and the Sichuan earthquake Agency-Research Team of GNSS-based geodetic tectonophysics and mantle–crust dynamics of the Chuan-Dian region (Grant no. 201803). This work was also supported by the Center for Astronautical Physics and Engineering (CAPE) from the Featured Area Research Center program within the framework of the Higher Education Sprout Project by the Ministry of Education (MOE) in Taiwan.

This paper was edited by Maria Ana Baptista and reviewed by two anonymous referees.