Exploring the Rheology of a Seismogenic Zone by Applying Seismic Variation

Featured Application

A 4D tomographic technique is utilized to investigate the variation of rheology around the fault zone.

Abstract

Although the study of spatiotemporal variation of a subsurface velocity structure is a challenging task, it can provide a description of the fault geometry as well as important information on the rheological changes caused by fault rupture. Our main objective is to investigate whether rheological changes of faults can be associated with the seismogenic process before a strong earthquake. For this purpose, a 3D tomographic technique is applied to obtain P- and S-wave velocity structures in central Taiwan using travel time data. The results show that temporal variations in the Vs structure in the source area demonstrate significant spatiotemporal variation before and after the Chi-Chi earthquake. We infer that, before the mainshock, Vs began to decrease (and Vp/Vs increased) at the hanging wall of the Chelungpu fault, which may be induced by the increasing density of microcracks and fluid. However, in the vicinity of the Chi-Chi earthquake’s source area, Vs increased (and Vp/Vs decreased), which may be attributed to the closing of cracks or migration of fluid. The different physical characteristics at the junctional zone may easily generate strong earthquakes. Therefore, seismic velocity changes are found to be associated with a subsurface evolution around the source area in Taiwan. Our findings suggest that monitoring the Vp and Vs (or Vp/Vs) structures in high seismic potential zones is an important ongoing task, which may minimize the damage caused by future large earthquakes.

1. Introduction

Changes in the physical properties of a fault zone before an earthquake rupture have become an important issue that scientists attempt to explore. Generally, the most popular methods for the topic of earthquake precursors are ionospheric anomalies [1,2,3] and seismic activity [4,5]. They can include variations in rock properties or pore pressures during the inter-seismic period, but these characteristics still cannot be directly observed around the seismogenic zone. Results of seismic velocity measurements reveal some information about the potential role of crustal rock strength and pore pressure in crustal dynamics. From a wave propagation simulation, it is known that Vp is inversely proportional to the crack density and Vs is proportional to the fractional volume of the fluid content in fractured rock [6]. In addition, for the velocity structures of Vp and Vs, the Vp/Vs ratio is correlated to Poisson’s ratio; therefore, they are regarded as important parameters in research into crustal rock characteristics [6,7]. Poisson’s ratio reflects the porosity, degree of crushing, and hydraulic pressure of rock formation.

Many factors such as the composition of rock, orientation of cracks, pore pressure, temperature, and fluid saturation would affect the Vp/Vs ratio. Since most seismic events occur in active fault zones, the Vp/Vs ratio can be used to determine if pressure changes will cause cracks to affect the velocity structure, which may result in gas or fluid migration and further large changes in pore pressure [8,9]. Therefore, the relationship between the characteristics of a fault zone and the distribution of the Vp/Vs ratio can be explored through three-dimensional (3D) tomography inversion. Furthermore, velocity structures can also be utilized to detect the attitude of subsurface formation and a brief overview of interface properties and features. Based on the above advantages, the results obtained in this study can be used to further understand the spatiotemporal evolution of fault materials.

Recently, some studies about seismic velocity variation have been conducted on the basis of the rupture of rocks caused by faulting [10], migration of magmatic fluids accompanying volcanic activity [11,12], and dilatancy and fluid migration before earthquakes [13], among others. The authors of [6] also described the variations in Vp before and after the 1999 Chi-Chi earthquake and the corresponding variations in Vp/Vs. Although there were no significant changes in Vp, after the occurrence of the Chi-Chi earthquake, deeper regions of the low-velocity anomaly appeared to have further expanded. The Vp/Vs model showed a substantial change reflecting an outward shift of the fluid-filled fractured source region. These changes are attributed to a localized stress change caused by the slip of the Chi-Chi rupture. The authors of [14] stated that the relationship between the velocity structure and time can be applied in earthquake forecasting and prediction, volcano warning systems, geothermal exploration, oil and gas reservoir assessment, and carbon dioxide capture and storage. However, the study of spatiotemporal variations of subsurface velocity structures is still a challenging work.

For the Taiwan region, from studies using 3D velocity models, similar results were discovered [15,16,17,18]; however, temporal variation was not considered in these 3D models. In 1999, the disastrous Chi-Chi earthquake M_L 7.3 occurred in central Taiwan and greatly affected Taiwanese society. In addition, although considerable research has been conducted to better understand the changes in the velocity structure in response to seismic events in central Taiwan, few studies have emphasized the significance of the velocity structure in relation to time. Therefore, this study aims to determine the impact of the Chi-Chi earthquake in central Taiwan before and after the event in order to understand the behavior of subsurface velocity structures. The main objectives of this research can be divided into three parts. First verifying the resolution of the model in each period. Next, detecting the uncertainty of the difference in different periods. Finally, the time evolution of the imaging velocity structure.

2. Geological Setting

Taiwan is situated on the collision boundary between the Eurasian and Philippine Sea Plates. Since the late Miocene, the southeast of the rifted continental margin, which was part of the Eurasian Plate, has been colliding with the Philippine Sea Plate, leading to the Taiwan orogeny [19]. The Philippine Sea Plate moves in a NW direction with a convergence rate of 8.2 cm/year [20]. The ongoing orogenic process causes a series of imbricate fold–thrust belts in the Western Foothills of Taiwan. Therefore, central Taiwan is composed of a complex system of belts, including the Changhua Fault (CHF), Chelungpu Fault (CLF), Shuangtung Fault (STF), Shuilikeng Fault (SKF), Lishan Fault (LF), and Longitudinal Valley Fault (LVF) (see Figure 1).

The great 1999 Chi-Chi earthquake (M_L = 7.3; depth = 8 km) that occurred at the CLF has received the most attention. The SKF is on the east of the foreland basin, and to its west are the STF, CLF, and CHF. They represent the major geological structures of central Taiwan. In addition, the continental margin of the Eurasian Plate contains several Eocene rift basins [21,22,23], of which the Hsuehshan Basin holds an important part of the Taiwan orogenic belt [21]. The sedimentary stratum beneath the foreland basin can reach more than 6 km in thickness, and the deformation front is also included [24].

The Eocene rocks of the Hsuehshan Basin now topographically occupy the higher ground and structurally overlie Pleistocene rocks of the foreland basin along the SKF [25,26]. Among them, the Peikang Basement High is located in western Taiwan, which impedes the orogenic process. Moreover, central Taiwan is the transition zone of the Peikang Basement High and the Central Range.

According to thin-skinned thrust theory [27], a basal detachment, composed of brittle structures, may exist beneath the fold thrust belts in the Western Foothills of Taiwan. During the collision of plates, structures above the basal detachment form a series of imbricate fold–thrust belts in the shallow layer. However, when the rift-related basin bounding faults were reactivated as a result of the collision process, the crust deformation proceeded toward the deeper parts and caused the underlying basement to be uplifted and exhumed [28,29]. This orogenic process can lead to destructive deformity on the surface, which leads to heavy casualties and property losses. In addition, the seismicity can occur along the reactivating rift-related faults, and steeply inclined seismic clusters extending into the lower crust are exhibited [28,30,31].

3. Data and Methods

In this research, we utilized the massive volume of the seismic events compiled by the Central Weather Bureau Seismic Network (CWBSN) from January 1991 to December 2003. We divided this time span into three segments to investigate velocity structures in different time periods, namely, 1 January 1991–31 December 1995; 1 January 1996–19 September 1999; 1 January 2001–31 December 2003 (hereinafter, period I, period II, and period III, respectively) (see Figure 1b). We excluded the days immediately before and after the Chi-Chi earthquake because of complicated source effects that make some seismic events unsuitable to be adopted as a point source; therefore, these data could not be used for tomography [32]. This study suggests that a large earthquake caused by stress adjustment makes the subsurface structure unstable, so the data for the year after the Chi-Chi earthquake was eliminated.

In this study, we used the travel times of P- and S-waves from earthquakes located inside the study area. Our study area stretched from 120.0° E to 122.0° E in longitude and from 23.5° N to 24.5° N in latitude. Earthquakes with magnitudes equal to or greater than 1.0 were considered, and data with poor location quality were excluded. Twenty-five stations were chosen in the study area, and at least 10 readings of each earthquake event were recorded, as the amount of observation of the matched earthquake event arrival affects the quality of the location of the earthquake; fewer observations result in poor location quality.

We applied seismic travel time data for tomography inversion by local tomography software [33]. The local tomography software was modified according to the bending-ray tracing algorithm by [34]. First, we utilized the Goal function [35] to detect the initial locations of the earthquakes. A straight path was drawn from the hypocenter to the station to locate the turning point at the midpoint and create a new path in sequence. Then, the corresponding travel time residue was calculated and the one with the minimum value was set as the ray path [28]. Finally, travel-time inversion was performed using the LSQR algorithm [36,37] to measure the P- and S-wave velocity anomalies and for obtaining the source parameters (four parameters for each seismic event). Furthermore, these results are fed back on a 3D grid for performing iterative inversion calculations until the iteration number or the travel-time residual value was less than the threshold value.

The study area was divided into numerous grids to invert the velocity variation. The grid was centered (121° E, 24° N) on the positive and negative sides of the x-, y-, and z-axes to 200 km and had dimensions of 5 × 5 × 5 km³. As a single-grid azimuth may cause artifacts [38], this study averaged the results from four angles—0°, 22°, 45°, and 67°—to reduce any grid dependency. We integrated different combinations of Vp and Vs damping values to obtain the optimal model through the trial and error process. Finally, we selected the appropriate Vp and Vs damping parameters that would give results with the stable travel-time residual (RMS). The different combinations of the Vp and Vs damping parameters for each stage are 15 and 15, respectively. In addition, the initial velocity model often plays a decisive role in the problem of velocity inversion. Therefore, for the selection of the velocity structure in this study area, we chose the one that closely approximates the real velocity structure model to facilitate the inversion. We applied the one-dimensional (1D) velocity model adopted by [17] as the initial 1D velocity structure in this study.

This formula refers to the method of [32], and the description is simplified. In the subsequent study, the following formula was applied on the temporal variation in the velocity structure:

$$\left[\raisebox{1ex}{$\left(V_n^{i+1}-V_n^i\right)$}\!\left/ \!\raisebox{-1ex}{$V_n^i$}\right.\right]\times 100\%$$

(1)

where V is the velocity, i is the number of periods, and n is the number of grids. Finally, the velocity perturbation model was obtained to explore the temporal variations in the velocity structure in central Taiwan.

4. Model Examination

An important objective of this study was to investigate the local velocity structural changes before and after the occurrence of strong earthquakes in central Taiwan. The evolution of velocity models in different time periods before and after the Chi-Chi earthquake is discussed. Through this analysis, we can further understand the seismogenic process before the earthquake, which is of great importance for future strong earthquake monitoring in central Taiwan.

In this study, the travel time data of the earthquake events compiled by the CWBSN were used, and the whole study time interval was divided into three periods for the inversion of 3D temporal velocity models for central Taiwan. To understand the reliability of the inversion results, a checkerboard test was used to verify the resolution at different depths.

We utilized the 1D velocity model adopted by [17] and added a ±8% perturbation to Vp and Vs to generate a new 3D velocity model in the form of a checkerboard (grid dimensions: 10 × 10 km² in horizontal space). Utilizing this 3D velocity model, the assumed travel times of P- and S-waves are calculated. The results yielded by the checkerboard test for each period are shown in Figure 2, Figure 3 and Figure 4.

The checkerboard test reflected the density of the ray paths at different depths, which could explain the resolution at different depths and confirm the reliability of the velocity model. The results show that both Vp and Vs have poor resolution at shallow depths of 0–5 km; this is because, near the surface, the ray path is nearly perpendicular to the station, resulting in an incomplete reflection of the lateral variation in the velocity structure. At depths of 10–25 km, excellent resolution was obtained because of a high ray density. These results suggest that the 3D velocity model of this area has a very high level of reliability for different time periods. The 3D velocity models of each period are shown in the Appendix A.

To verify the impact of artifacts, this study designed another test [32] restoration by identical synthetic models using two data subsets corresponding to different periods (Figure 5, Figure 6 and Figure 7).

The test results using data of periods I and II in the proposed model with the variation in the synthetic manner are shown in Figure 8.

The grid settings for the two cases were the same. After restoring models in these two time periods, we computed the differences between the recovered models, representing the Vp and Vs variations, at different depths. First, Figure 8a,b, respectively, show that the Vp and Vs results are similar, except for the border area, where the variation was 4%. In other areas, the variations were less than 0.40% and 0.54% on average. Next, Figure 9a,b, respectively, show that the Vp and Vs results are similar, except for the border area, where the variation was 4%. In other areas, the variations were less than 0.54% and 0.64% on average.

We set the white area as resolvable, where the average variation was less than the actual temporal evolution of the velocity structure. These results revealed that the change over time was reliable. The blue and red areas represented greater variation, which implies greater effect of artifacts. In contrast, no temporal resolution capability could be shown in the black area. Furthermore, the variation in the synthetic models in periods II and III (Figure 9) was similar to that in periods I and II (Figure 8). However, the variations in Vp and Vs increase, which is attributed to the aftershocks of the Chi-Chi earthquake. Therefore, special attention must be given to the subsequent interpretation of the results.

As will be explained in the following section, the seismic velocity variations obtained by the model were less than those obtained from actual observation in the data analysis. This implies that the temporal variation in our model is not significantly affected by the artifacts associated with changes in the ray configuration for the selected data subset. However, a much wider range of artifact effects must be considered when exploring local changes in the study area. The presented test results demonstrate that, for the study region, the proposed model robustly revealed the variations in the crust with a scale of more than 2%. The results of the experimental data analysis for stronger variations, presented in the next section, thus represent variations that occur in the crust as a result of tectonic processes.

5. Results and Discussion

5.1. Temporal Evolution of a Velocity Structure

The pore pressure is inversely proportional to the seismic wave velocity, evidenced by oil field detection and laboratory experiments [39]. Moreover, many experiments have examined the effects of fluid-filled cracks on seismic velocity [40,41,42,43]. Therefore, it is generally accepted that Vp is inversely proportional to the crack density and that Vp/Vs is proportional to the fluid saturation in the crack. In addition, its correlation with the rock crack and pore pressure was discussed. Fortunately, there have been studies using more realistic observational data to build fracture models [44] and there are detailed discussions on the differences between various fracture models.

First, variations in Vp, Vs, and Vp/Vs in the α zone (source area, see Figure 10b) before the occurrence of the Chi-Chi earthquake sequence in periods I and II are discussed and are shown in Figure 10.

The results indicate that Vp decreased by approximately 2%, Vs decreased by approximately 4%, and Vp/Vs increased by approximately 4%. In the β zone (around the source area, see Figure 10b), Vp increased by approximately 4%, Vs increased by approximately 10%, and Vp/Vs decreased by approximately 6%. The results demonstrate that the subsurface structure had already changed before the impending Chi-Chi earthquake, and the reason for the significant decrease in Vs in the α zone compared with Vp is possibly the proliferation of cracks and fluid intrusion. The authors of [45] also speculated using a magnetotelluric technique that the seismogenic mechanism of the Chi-Chi earthquake may be associated with fluid intrusion. The velocity structure variation in the β zone was opposite to that in the α zone, presumably due to the closure of the crack, resulting in an increase in Vs. Figure 10b,c clearly show that, at depths of 5–10 km, the hypocenter of the Chi-Chi earthquake occurred in the boundaries of the areas of high and low Vp/Vs. However, no significant difference was noticed at depths greater than 15 km.

From periods II and III, the variations in Vp, Vs, and Vp/Vs after the earthquake sequence are shown in Figure 11.

The results demonstrate that Vp and Vs varied within ±10% in the vicinity of the CLF. These variations were more complicated in the shallow crust (equal to or less than 15 km). Note that in Figure 11a, the result shows a significant decrease in Vp at a depth of 5–15 km in the northern segment of the CLF. The distribution of the low Vp anomaly is in good agreement with the result of maximum slip distribution in the GPS calculation [46,47], because the slip vector in the GPS calculation is the combined effect from the surface observation. However, from the temporal variations in the velocity structure, we can estimate the possible slip pattern beneath the subsurface. When the depth was less than 15 km, the lateral variation in Vp was more prominent than in the lower crust, which reflected changes in rock cracks after the earthquake occurred. No more lateral change was noted for depths greater than 25 km, which indicated that the lower crust was less affected by the Chi-Chi earthquake. At a depth of more than 15 km, the variations in Vs and Vp/Vs near the hypocenter reduced after the earthquake, and the area with obvious velocity variation was located in the area north of the CLF to LF.

5.2. Tectonic Implications

In this study, we discussed the relationship between the velocity structures of Vs and Vp/Vs over time and the geometry of faults in order to determine the variations in Vs and Vp/Vs before and after the occurrence of the Chi-Chi earthquake sequences along the velocity profiles. Figure 12 shows the time-dependent perturbation profile in periods I and II. It can be seen from the section of CC*–EE* connection that in the α area (hypocenter) during the Chi-Chi earthquake, Vp/Vs increased by approximately 4% and Vs decreased by approximately 4%.

In the β area (near the hypocenter), Vp/Vs decreased by approximately 6% and Vs increased by approximately 10%. The fluid-filled fractured region at the α area (hypocenter) was assumed to lead to fissure expansion; therefore, the decline of Vs was greater than Vp. Furthermore, in the south to north sections of Vs, Vs changed along an eastward-tilted inclined plane. Moreover, Vs changed more significantly south of the hypocenter compared with the north. This phenomenon occurred before the Chi-Chi earthquake and was concentrated south of the CLF. However, the variations in the β area were opposite to those in the α area.

In contrast, in this study, we discuss the velocity variation over a long time scale before the occurrence of the earthquake sequence; we found remarkable changes in periods I and II, such as the decrease in Vs in the α region and the increase in Vs in the β region, while the hypocenter was located at the junction of the α and β areas. This phenomenon indicated that the direction of the tectonic stress was from the south-east (approximately 100°, [48]), resulting in the closure of the pores in the β area, which may allow the fluid to migrate toward the α area. The decrease in the frictional strength of the fault in this area might be attributed to the vicinity of the α area to the hypocenter, the increase in the pore density, and the possibility of fluid intrusion. These factors may have led to the occurrence of the Chi-Chi earthquake, suggesting it is a precursor phenomenon.

From periods II and III, the time-dependent perturbation profile of Vs is shown in Figure 13; it can be seen from the section of AA*–EE* connection in the low anomalous regions that Vs significantly expanded from south to north, revealing that the northern section of the CLF has a larger slip displacement after the rupture, which leads to an east-dipping, low-Vs anomalous region.

Furthermore, the slip displacement of the southern section was smaller than that of the northern section, and Vp/Vs exhibited fewer anomalies compared with those in the northern section. This anomalous region is speculated to be the Chi-Chi earthquake rupture zone, which is consistent with the location of the CLF. Further analysis is required to explore the spatial variation of the earthquake rupture, the results of which would also be valuable for confirming the larger rupture zone in the northern section of the CLF compared to the southern section.

In addition, the authors of [49] independently analyzed the Qs of earthquakes and found that Qs changed with time. Among them, the region where Qs is higher before the Chi-Chi earthquake corresponds to the β region where Vs is increased in Figure 12a of the present study. The lower area of Qs after the Chi-Chi earthquake is corresponding to the Vs falling area in Figure 13a of this study. They speculate that the decrease in Qs (increased S-wave attenuation) may indicate an increase in fluid content in the ruptured area after the Chi-Chi earthquake which resulting in increased pore fluid saturation.

6. Conclusions

This research used travel time data for earthquake events compiled by the Central Weather Bureau Seismic Network (CWBSN). Using these data, tomography inversion was conducted for different periods of structures in central Taiwan. The results were utilized to deduce the relationship between seismicity and fault structure with spatiotemporal evolution. Based on the aforementioned information, the following conclusions can be drawn.

• The results showed that the temporal variations of Vp were less than those of Vs, and Vp was not affected significantly before the occurrence of the Chi-Chi earthquake.
• The velocity profile indicates that the temporal variation of the subsurface velocity structure may have been affected before the occurrence of the Chi-Chi earthquake sequences. In the α area, the decreased Vs and increased pore pressure are attributed to the increasing density of microcracks and the possibility of fluid or gas intrusion. Moreover, we observed increasing Vs in the β area; this was because the β area was located on the eastern side of the α area. It was presumed that the direction of the tectonic stress was from the southeast, resulting in the closure of the pores and fluid migration in the β area. This phenomenon may be attributed to the different mechanism of action in the α and β areas.
• Based on the present results, the temporal variations reflect precursory phenomena in the source area before the Chi-Chi earthquake and the range of the CLF affected by the Chi-Chi earthquake.