Dynamic rupture process of the 1999 Chi-Chi, Taiwan, earthquake

1.   Introduction

It is expected that dynamic rupture processes on faults intersecting the ground surface with small dip angles would be strongly influenced by the ground surface. The 1999 M_W7.6 Chi-Chi, Taiwan, earthquake is a typical example. This earthquake is the largest event in Taiwan in the 20th century, and is well recorded by more than 400 three-component acceleration seismometers all around the island. Field survey after the mainshock revealed that this earthquake caused more than 100 km surface rupture along the Chelungpu fault, with a maximum vertical offset more than 8 m (Shin and Teng, 2001).

Abundant data collected during and after this earthquake provide unprecedented chance for seismologists to investigate the source of a large earthquake in detail. Numerous studies based on different inversion methods have been carried out on the detailed rupture process of this earthquake (e.g., Iwata et al, 2000; Ma et al, 2001; Wu et al, 2001; Xu et al, 2002; Wang, 2003). Although the detailed patterns of rupture process and slip distribution in these studies are slightly different, the overall patterns are consistent. For example, it is generally accepted that slip on the northern part of Chelungpu fault is much larger than that in southern part, and the rupture began to propagate outwards at some time (about 6 s later) after it appeared, and the whole rupture process last about 26 s. This kind of slip distribution is responsible for the difference in the amount of displacement on ground between hanging wall and footwall (Oglesby and Day, 2001a, b).

On the other hand, however, kinematic inversion studies cannot provide explanations to the phenomena they revealed. To this end, one must fall back on dynamic modeling. Since in dynamic studies, rupture processes and slip distributions are the results of certain friction law and physical state in the source region, they can provide explanations for the physical nature of source processes. Nevertheless, dynamic studies on the 1999 Chi-Chi earthquake are relatively rare. Some authors, such as Oglesby and Day(2001a, b) and Zhang et al (2003, 2004) among others, have investigated the dynamic process of this earthquake by finite element and finite difference method, respectively. Oglesby and Day(2001a, b) suggested that some near-field observation features, such as movement of hanging wall is much larger than footwall and mechanisms in southern and northern part of the fault are thrust and strike, respectively, can be explained as the simple result of the geometry of the fault without assuming any heterogeneous distribution of stress, complex friction law and site effects. Based on kinematic inversion results, Zhang et al(2003, 2004) performed dynamic modeling and concluded that the rupture process of the Chi-Chi earthquake roughly followed the slip-weakening friction law, and obtained complex distribution of slip and critical slip-weakening displacement. These results are helpful to an in-depth understanding of the source of this earthquake.

In this study, we investigate the dynamic process of the 1999 Chi-Chi earthquake by applying an extended boundary integral equation method (BIEM) we developed recently (Zhang and Chen, 2006a, b; Chen and Zhang, 2006). Unlike previous study based on traditional BIEM (e.g., Fukuyama and Madariaga, 1998; Aochi et al, 2000), the extended BIEM can naturally included the effect of free surface, which is achieved by applying the half-space Green's function for elastodynamic equation. It is therefore very suitable for modelling rupture processes on faults intersecting the ground surface with small dip angle, such as the Chelungpu fault, which is difficult to handle by using traditional BIEM. After a brief introduction to the extended BIEM, we first carefully assign the model parameters based on published literature, and then compare the numerical results with existing kinematic (Iwata et al, 2000; Ma et al, 2001; Wu et al, 2001; Xu et al, 2002; Wang, 2003) and dynamic (Oglesby and Day, 2001a; Zhang et al, 2004) results. Finally, we discuss possible further improvement on the current modeling.

2.   Method

Details of the extended BIEM are referred to Zhang and Chen (2006a) and Chen and Zhang (2006). We only briefly review the method here.

The extended BIEM is developed for dynamic ruptures on shallow faults in half space, especially those run up to the ground surface. For faults embedded more than 1 km or more, the effect of ground surface is so small that can be neglected (Oglesby et al, 2000; Zhang and Chen, 2006b), and a traditional BIEM based on full-space Green's function (e.g., Aochi et al, 2000) is an appropriate alternative. The main difference between the extended BIEM and the traditional one is the Green's function they use. In the extended BIEM, we use the Green's function for 3-D elastic half space, which can only be expressed in double integrals (e.g., Kennett, 1983; Chen, 1999) rather than in an analytical closed form as in Green's function for full space. It is therefore more difficult in handling the dynamic ruptures on faults embedded in half space. Besides, hypersingularities in integrals cannot be avoided as in Aochi et al (2000) and must be carefully treated.

In Zhang and Chen (2006a), we systematically presented the theoretical development, and obtain the extended boundary integral equations as follows:

where τ'_3α^ij0 and τ'_3α^(abs)ijk are initial and absolute shear stresses on the fault, respectively. V'_α^lmn is slip rate on the fault, and C_α^{ijk; lmn} is the integral kernel. α=1 or 2 indicate the direction of slip. Integer pairs (i, j) and (l, m) indicate the locations of observation point and source point on the fault, respectively, and k and n are the serial number of instants of observation point and source point, respectively. According to equation (1), shear stress at an arbitrary kth instant on the fault can be expressed as a weighted summation of slip rate on entire fault at present and all previous instants. In equation (1), it is assumed that the CFL (Courant-Friedrichs-Lewy) ratio (see, e.g., Madariaga et al, 1998) w=v_pΔt/Δs≤1/2 (here, v_p, Δs and Δt are the P-wave velocity, space step and time step, respectively), therefore contribution to the stress comes only from the current element within the same time step (second term on the right hand side, which is called instantaneous term). Explicit separation of the instantaneous term greatly reduces the computation, because only an algebraic equation, rather than complicated simultaneous equations, must be solved. The form of equation (1) is exactly the same as that in traditional BIEM (e.g., Aochi et al, 2000). However, the integral kernels, which are the most important factor, are completely different.

For a spontaneous rupture problem, the BIEs alone are not enough. An appropriate friction law must be incorporated. As in most dynamic modeling studies (e.g., Andrews, 1976, 1985; Harris and Day, 1993, 1999; Olsen et al, 1997; Fukuyama and Madariaga, 1998; Oglesby and Day, 2001a, b; Aochi et al, 2000; Zhang et al, 2004), we adopt the slip-weakening friction law, which is proved to be the dominant friction law in the 1999 Chi-Chi, Taiwan, earthquake (Zhang et al, 2003).

3.   Parameters setting

In numerical modeling, parameters setting are very important, because the results rely strongly on the setting. In this section, we carefully assign computation parameters based on the published literature (Iwata et al, 2000; Oglesby and Day, 2001a; Shin and Teng, 2001; Ma et al, 2001; Xu et al, 2002; Wang, 2003; Zhang et al, 2003, 2004). Since the current extended BIEM can only deal with planar fault, while the Chelungpu fault is a typical nonplanar one with the northern part turns eastwards (Ouchi et al, 2001), equivalent treatment on initial stress and shear strength is necessary. No other heterogeneity is considered other than that.

3.1   Geometrical parameters

Geometrical parameters include the length (L) and width (W) of the fault, dip angle of the fault (δ), space element size (Δs), location of focus (d and H denote the distance between the focus and south boundary of the fault and focal depth, respectively), and the size of asperity (for simplicity, we assume the initial region of the rupture is a circle with a radius R_asp).

Geometry of the fault is different in different studies. A nonplanar fault is adopted in Iwata et al (2000) and Zhang et al (2004). The size of projection of the fault on the ground is 78 km×39 km, and the dip angle of the main part of the fault is 29°. A similar spade fault is assumed in Wang (2003), and the size of projection of the fault on the ground is 70 km×40 km. Wu et al (2001) constructed a fault model consisting of three segments. The size of the main subfault is 84 km×44 km, and the dip angle is 30°. On the other hand, a planar fault with a size 170 km×90 km and a dip angle 30° is assumed in Xu et al (2002). In this study, we adopted a fault model suggested by Shin and Teng (2001). The size of the planar fault we consider is 100 km×40 km and the dip angle is 30°. Space element size Δs is set to be 2 km (hence the fault is divided into 50×20 elements). Although this value seems a bit large, considering the size of the Chelungpu fault, overall patterns of the rupture process can be revealed according to the resolution test (Appendix D in Chen and Zhang, 2006).

The distance between the focus and south boundary of the fault (d) is generally considered as 30 km (Oglesby and Day, 2001a, b; Xu et al, 2002), or some other value around 30 km (Iwata et al, 2000; Ma et al, 2001; Wang, 2003; Zhang et al, 2003, 2004). In the study, we set d=30 km. The focal depth H is a quantity difficult to determine accurately. For example, H determined by USGS, Harvard and CWB are 5 km, 21 km and 8 km, respectively. Most studies adopt a value between 5 km and 15 km (Iwata et al, 2000; Ma et al, 2001; Oglesby and Day, 2001a, b; Wu et al, 2001; Wang et al, 2003; Zhang et al, 2003, 2004). Some authors (e.g., Xu et al, 2002), however, used a larger value (H=21 km). In this study, we adopt the value suggested by CWB, that is, H=8 km.

As for the radius of asperity R_asp, since it is a quantity involved in dynamic modeling, studies that can be referred to are rare. Zhang et al (2004) assumed that the rupture of Chi-Chi earthquake started from a 9 km×9 km region around the focus. Considering the size of the fault, we assume here that the asperity is a circle with a radius R_asp=4 km.

Based on the above setting, a fault model is constructed as in Figure 1. The 100 km×40 km fault intersects the ground surface with a dip angle 30° (Figure 1a), and the rupture is assumed to start from the grey region (Figure 1b). Geometrical parameters setting in published literature and this study discussed above are summarized in Table 1.

Figure  1.  Geometry of the fault in this study. (a) 3D view of the fault. Ground surface is in the plane Ox₁x₂. The fault in grey color intersects the ground surface, and intersection is along axis x₁. Dip angle of the fault is 30°. The black dot indicates the focus; (b) 2D view of the fault. The 100 km×40 km fault is divided into 50×20 elements. The region in grey around the black dot (focus) is the asperity, from which the rupture initiates.

DownLoad: Full-Size Img PowerPoint

Table  1.  Geometrical parameters setting in published literature and this study

 | Show Table

DownLoad: CSV

3.2   Physical parameters

Physical parameters include the shear strength on the fault (T_u), initial stress inside and outside of the asperity (denoted as T_i and T_e, respectively), and the critical slip-weakening displacement (D_c).

A simple planar fault model may not well interpret the near-field coseismic data of Chi-Chi earthquake, because field survey revealed that the northern part of the Chelungpu fault turns eastwards (Ouchi et al, 2001), implying that the fault consists of several segments. A non-planar fault model is therefore adopted in many kinematic inversion studies (e.g., Iwata et al, 2000; Wu et al, 2001; Wang, 2003). In the current extended BIEM, however, only a planar fault model is considered due to the complicated Green's function. In order to fit the kinematic inversion results, an equivalent treatment on the geometry of fault is necessary. Aochi et al (2000) studied the dynamic rupture problem in a complex fault system in full space by using a traditional BIEM, and found that even in a homogeneous triaxial compressional field, rupture will automatically slow down or stop due to the heterogeneous distribution of initial stress, which is a result of variation of the fault strike. It is therefore reasonable to replace the effect of fault bend with a higher shear strength and initial shear stress (assuming the increase in shear strength is larger than that in initial shear stress) in a planar fault model, as also did in Oglesby and Day(2001a, b). Based on the results of kinematic inversion, Zhang et al(2003, 2004) obtained a peak stress drop of 35 MPa, and found that the tectonic shear stress before the mainshock is close to the shear strength. According to the above analysis, we construct a simple distribution of shear strength and initial shear stress as shown in Figure 2, and assume that the initial stress inside the asperity just exceeds the shear strength (T_i=1.001T_u), such that the rupture is able to start from this region. Shear strength around the fault (except the upper bound, which intersects the ground) is assumed to be very large (10T_u), so that the rupture can stop automatically after it reaches the boundary of the fault.

Figure  2.  Distributions of shear strength (a) and initial shear stress (b). Vertical and horizontal axes are the down-dip distance and the distance along strike, respectively (see Figure 1b).

DownLoad: Full-Size Img PowerPoint

The direction of the initial stress is not along the strike or dip. For instance, in Oglesby and Day(2001a, b), the angle between initial stress and strike is 67°. However, the current extended BIEM cannot deal with an initial stress with arbitrary direction, and only a pure strike- or dip-slip is considered instead. Therefore, an equivalent treatment on initial stress distribution is also needed. We decompose the initial stress in the direction of dip- and strike-slip, and discuss the pure dip-slip and strike-slip problems. Although the slips obtained in such way are not the components of actual slips, it is expected that they represent the corresponding components of the latter to a great extent. Because according to Madariaga et al (1998), slip occurs mainly in the direction of stress drop, implying that a scalar friction law is a good simplification of a vector friction law.

An appropriate friction law is crucial for dynamic modeling. As stated in Fukuyama and Madariaga (1998), it is friction that controls the initiation and development of rupture and the healing of the fault. As for the Chi-Chi earthquake, Zhang et al (2003) analyzed the relationship among stress, slip and slip-rate, and found that behaviors of stress and slip during the earthquake obey the slip-weakening friction law, and the observation data do not support the rate-weakening friction law. Therefore, we assume that the rupture process is controlled by a slip-weakening friction law throughout this study. Critical slip-weakening displacement D_c is the most important parameter in slip-weakening friction law. According to Zhang et al (2003), distribution of D_c is quite heterogeneous, and the peak value is about 10 m. In the dynamic modeling of Oglesby and Day(2001a, b), however, D_c is assumed to be 5 cm. It is obvious that the freedom in choice of D_c is quite large. Taking this into account, we do not fix the value of D_c in advance, instead, we adjust D_c by trial and error such that the computed dynamic rupture process is roughly consistent with the result of kinematic inversion. The D_c obtained in this way is expected to be close to the actual value, so long as other parameters are assigned properly.

3.3   Other parameters

Other parameters required in the numerical computation include shear modulus μ, density ρ, Possion ratio v and CFL ratio w_α.

For medium parameters, we take those in Oglesby and Day(2001a, b) as a reference, and assume that μ=3×10¹⁰ Pa, ρ=3.0 g/cm³, and ν=0.25. And therefore the P wave and S wave velocities α and β are therefore 5.48 km/s and 3.16 km/s, respectively. Equation (1) is valid when w_α < 0.5 (Fukuyama and Madariaga, 1998; Aochi et al, 2000), and we set w_α=0.3. Since the space step Δs=2 km, the time step is obtained as Δt=0.11 s. The maximum frequency is f_max=20 Hz. It is noted that such a high frequency is impossible for other numerical method, such as the finite difference and the finite element method. It is, however, can be reached by using the BIEM, implying that more details of the rupture process can be revealed.

4.   Numerical results

In this section, we will first show the results obtained by numerical modeling with parameters assigned in section 3, and then compare them with those of kinematic inversion and dynamic modeling in some published studies.

4.1   Numerical results in this study

As stated in section 3.2, we have to calculate pure dip- and strike-slip problems by using the current extended BIEM.

Figures 3a and 3b show the snapshots of slip-rate and slip for dip-slip, respectively. Since T_i just exceeds the local shear strength, some time is required to accumulate energy for rupture to propagate outwards. Apparent propagation of rupture occurs around t=6.0 s. It is noted that this is an important feature inferred by kinematic inversion (e.g., Shin and Teng, 2001). In order to fit this feature, we adjusted D_c by trial and error and found that with other parameters fixed, D_c must locate in a very narrow range (results in Figures 3 and 4 are obtained with D_c=0.64 m), otherwise, the above feature cannot be satisfied. For instance, rupture will soon begin to propagate (at about t=3.0 s) when D_c=0.50 m, while it takes a long time to propagation (about 10.0 s) when D_c=0.80 m. And rupture will die away soon after it begin to propagate when D_c=0.96 m.

Figure  3.  Snapshots of slip rate (a) and slip (b) for dip-slip. The black dots indicate the focus. Vertical and horizontal axes are the down-dip distance and the distance along strike, respectively (see Figure 1b).

DownLoad: Full-Size Img PowerPoint

It can be seen from Figure 3 that the rupture propagates outwards after t=7.67 s. Rupture front on the in-plane direction soon reaches the ground surface, and a new rupture front due to the reflection of ground surface occurs (t=10.96 s and t=14.24 s). At t=17.53 s, the rupture front propagating southwards has already reached the fault boundary, and a stopping phase is generated. After that, the rupture can be regarded as a unilateral rupture. It is noted that shear strength and initial shear stress on the right lower part of the fault are both larger, and the increase in shear strength is larger than initial stress, therefore the velocity of rupture front on the lower part of the fault is lower than that of the upper part, while the slip-rate on lower part is much larger. This results in a larger slip on the right lower part of the fault after t=24.1 s (Figure 3b). And the final slip on the northern part of the fault is much larger than that on the southern part (t=37.25 s). Especially, maximum slip on the ground surface is larger than 15 m on the northern part of the fault, while the corresponding value on the southern part is much smaller.

Figures 4a and 4b show the snapshots of slip-rate and slip for strike-slip, respectively. The overall pattern is similar to that for dip-slip. Since the in-plane direction is parallel to the ground surface, rupture speed is markedly larger than that for dip-slip. If the initial shear stress is close to the shear strength (Zhang et al, 2003), it is no wonder that the rupture speed of front is supershear (Zhang and Chen, 2007). In fact, in region without equivalent treatment, we assume that T_e/T_u=0.9. Therefore, the rupture speed of the front in in-plane direction soon reaches the supershear speed after it begins to propagate (about t=6.13 s). Similar to the case of dip-slip, final slip on the northern part of the fault is apparently larger than that on the southern part (t=28.81 s).

Figure  4.  Snapshots of slip rate (a) and slip (b) for strike-slip. The black dots indicate the focus. Vertical and horizontal axes are the down-dip distance and the distance along strike, respectively (see Figure 1b).

DownLoad: Full-Size Img PowerPoint

It is mentioned previously that the results of pure dip- and strike-slip are not strictly the dip and strike components of the actual problem, but dip and strike components can be replaced by corresponding pure slip respectively to some extent. From the above analysis, we can see that the evolution patterns of dip- and strike-slip are very similar, and therefore it is expected that the actual rupture process is also similar. The evolution of total slip should be faster than the dip-slip case and slower than the strike-slip case considered above.

4.2   Comparison with results of other studies

In studies based on kinematic inversion, slip on the fault is inferred by directly fit the recorded data, it can therefore be regarded as a restriction of dynamic modeling. However, due to the non-uniqueness of inversion problem, it is not infrequent that different authors obtained very different results. Results in a single study can hardly be regarded as an absolute standard. We take the common feature of different kinematic inversion studies, which are based on different methods and data, as a restriction of our dynamic modeling.

Final slip distributions from four different studies of kinematic inversion (Iwata et al, 2000; Ma et al, 2001; Wu et al, 2001; Xu et al, 2002) are shown in Figure 5. Iwata et al (2000) used a non-planar fault model and obtained a very heterogeneous final slip based on data from 31 near-field stations (Figure 5a). Ma et al (2001) and Wu et al (2001) performed a hybrid inversion of near-field data, far-field data and GPS data, and obtained quite similar results (Figure 5b and 5c). A relatively rough result is obtained by Xu et al (2002) based on long-period far-field data, taking the aftershocks as empirical Green's functions (Figure 5d). Although the details of these results are different due to different data and different inversion method, the overall pattern are the same, that is, slip concentrates mainly on the northern and lower part of the fault. We obtained the same results in this study.

Figure  5.  Distribution of final slip obtained in several studies based on kinematic inversion. (a) Iwata et al (2000); (b) Ma et al (2001); (c) Wu et al (2001); (d) Xu et al (2002).

DownLoad: Full-Size Img PowerPoint

Dynamic modelings of the Chi-Chi earthquake are relatively rare. We take the results of two typical studies (Oglesby and Day, 2001a; Zhang et al, 2004) as an example. Oglesby and Day (2001a) studied the dynamic process of the Chi-Chi earthquake based on half-space medium and planar fault model by using finite element method. Figure 6a shows the displacements of strike and dip component on hanging wall and footwall, respectively. Difference of displacements of hanging wall and footwall is the slip on the fault. They assumed a circular region with a stress drop 10 times higher than other region to fit the results of kinematic inversion. It can be inferred from Figure 6a that the slip distribution is very similar to ours, that is, the slip of the fault is much larger on the northern part. Zhang et al (2004) obtained a quite complex dynamic process based on "thick" fault and layered media model by using an improved finite difference method. Figure 6b shows the snapshots of slip. Since their results were obtained based on heterogeneous dynamic parameters inferred from kinematic inversion, the slip distribution is strongly heterogeneous and roughly consistent with that of Iwata et al (2000). Nevertheless, the main feature is still the same as our results.

Figure  6.  Results of dynamic studies. (a) Distribution of particle displacement on footwall and hanging wall obtained in Oglesby and Day (2001a); (b) Snapshots of slip obtained in Zhang et al (2004).

DownLoad: Full-Size Img PowerPoint

5.   Discussion

Although the effect of ground surface can be exactly included in the dynamic modeling, the model in this study is much too simple for a real earthquake.

We used a planar fault model in this study. However, outcrop of the Chenglungpu fault indicates that this fault is no way a planar fault. Taking this into account, we perform an equivalent treatment as stated in section 3.2. Figures 7a and 7b show distributions of final slip for dip-slip (left in each panel) and strike-slip (right in each panel) without and with equivalent treatment, respectively. According to Figure 7, a nearly symmetric slip distribution is obtained without equivalent treatment, while the pattern with equivalent treatment is completely different, implying that fault geometry plays an important role in the dynamic rupture process. It must be point out that such equivalent treatment is a pis aller under the current scheme. A complex fault model can be included without difficulty in principle, and this is a project to be studied in a near future.

Figure  7.  Distributions of final slip for dip-slip (left in each panel) and strike-slip (right in each panel) without equivalent treatment (a) and with equivalent treatment (b). The black dots indicate the focus.

DownLoad: Full-Size Img PowerPoint

In this study, we use a half-space medium model. Although it is an important progress to a real modeling compared with the traditional BIEM (e.g., Fukuyama and Madariaga, 1998; Aochi et al, 2000), it is still too simple for an actual problem. For example, slip on the ground inferred by kinematic inversion (Figure 5) is not as large as expected by dynamic modeling (Figures 3 and 4). Because when the rupture runs up to the ground, the effect of low-speed layer close to the surface cannot be neglected. A more real medium model might be a half-space covered by a low-speed layer. It is expected, however, the Green's function is even more complex in such case.

Due to the simplification assumption (slip occurs only in the direction of stress drop, Zhang and Chen, 2006), we cannot deal with the case of stress drop with an arbitrary direction by the current extended BIEM. Instead, we consider pure strike-slip and dip-slip respectively in this study. In fact, when the stress drop in a direction is considered, the slip on other directions cannot be supposed to be zero, though the value is in fact very small. A further improvement is to adopt a vector friction law rather than a scalar one.

In the slip-weakening friction law, we assume the residual stress to be zero for simplicity. This, however, is not the actual case, because this cannot explain the occurrence of aftershocks. Furthermore, we do not consider the effect of normal stress, as in Aochi et al (2000), which may also play a role in the dynamic process.

In a word, this study on such a complex large earthquake as the Chi-Chi earthquake is quite preliminary, and many problems are still open to further study.

6.   Conclusions

We preliminary studied the dynamic rupture process of the 1999 Chi-Chi, Taiwan, earthquake by applying the extended BIEM we recently developed (Zhang and Chen, 2006a; Chen and Zhang, 2006). By carefully assigning the parameters based on published studies, we took the common feature inferred by kinematic inversion as a restriction and adjusted D_c by trial and error, and found that: ① To fit the results of kinematic inversion (i.e., rupture begins to propagate outwards after about 6 s, and the initial shear stress is close to the local shear strength), D_c must located in a narrow region [60cm, 70cm]; ② With the equivalent treatment on the fault geometry, the final slip distribution is quite heterogeneous, i.e., slip on the northern part of the fault (especially the northern lower part) is much larger than other regions on the fault; ③ If the conclusion "initial shear stress before the mainshock is close to the local shear strength" comes into existence, supershear rupture should occur during the rupture process according to the results of dynamic modeling.

We neglect all the heterogeneities of parameters in this study, except some necessary equivalent treatment. Even with such simple assumptions, the main feature of rupture is still consistent with that inferred by kinematic inversions. It is therefore can be expected that geometry of the fault, such as the intersection with the ground surface and the non-planar geometry, plays a very important role in the dynamic process of the Chi-Chi earthquake. This conclusion supports that in Oglesby and Day (2001a).