# Slip distribution on a thrust fault at a plate boundary: the 2003 Chengkung earthquake, Taiwan

L. Mozziconacci, B. Delouis, J. Angelier, J.-C. Hu, B.-S. Huang

► To cite this version:

L. Mozziconacci, B. Delouis, J. Angelier, J.-C. Hu, B.-S. Huang. Slip distribution on a thrust fault at a plate boundary: the 2003 Chengkung earthquake, Taiwan. Geophysical Journal International, 2009, 177 (2), pp.609 - 623. 10.1111/j.1365-246X.2009.04097.x. hal-00497741

HAL Id: hal-00497741

https://hal.science/hal-00497741v1

Submitted on 15 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Slip distribution on a thrust fault at a plate boundary: the 2003 Chengkung earthquake, Taiwan

Laetitia Mozziconacci,<sup>1,2</sup> Bertrand Delouis,<sup>1,2</sup> Jacques Angelier,<sup>1,2,3</sup> Jyr-Ching Hu<sup>2,4</sup> and Bor-Shouh Huang<sup>2,5</sup>

Accepted 2008 December 21. Received 2008 October 12; in original form 2008 May 14

## SUMMARY

The 2003 Chengkung earthquake ($M_w$ 6.8) occurred on the east-dipping Chihshang fault of eastern Taiwan, a listric fault with dominant reverse motion. A joint inversion taking into account geodetic (GPS) and seismological (accelerometric and teleseismic) data is performed. This modelling highlights the rupture process behaviour on the fault plane at depth. The coseismic rupture developed essentially between 12 and 26 km depth, expanding laterally by about 30 km towards the SSW. Two slip patches dominated, with slip values exceeding 1.5 m. During the 30 s of the whole process, the rupture propagation decelerated with time, starting at $3.4 \pm 0.3$ km s$^{-1}$ in the two main slip patches and ending with velocities below $1.6 \pm 0.1$ km s$^{-1}$. The 48 hr of aftershocks that followed the main shock are distributed essentially inside the coseismic slip patches on the fault surface, except in the deepest, southern part of the fault where slip is less resolved. The results of our modelling illustrate the particular behaviour of the 2003 Chengkung earthquake where the coseismic slip rapidly diminished towards the surface, in good agreement with the observation of fast post-seismic slip followed by interseismic creep along the Chihshang fault.

**Key words:** Time-series analysis; Earthquake dynamics; Earthquake ground motions; Earthquake source observations; Dynamics and mechanics of faulting; Asia.

## 1 INTRODUCTION

The Taiwanese orogen is the consequence of the collision between the Philippine Sea Plate and the Eurasian continental margin (Fig. 1a). The suture zone between these two plates follows the Longitudinal Valley not far from the eastern coast of Taiwan (Ho 1986; Teng 1990). There, the Coastal Range is thrusted over the Central Range along the Longitudinal Valley Fault (LVF) (Figs 1a and b). The LVF shows active thrusting with a minor component of left-lateral slip (Fig. 1a), as shown by surface studies of its ongoing creep displacements along the fault trace (Angelier et al. 1997, 2000). The Chihshang Fault (Ho 1986; Tsai 1986), which was ruptured during the 2003 Chengkung earthquake (Fig. 1b), is the main segment of the LVF in the South.

The Chengkung earthquake occurred on 2003 December 10, 04h38 (UTC) near the coastal town of Chengkung (Fig. 1a). The epicentre was relocated at 121.324°E–23.106°N (Kuochen et al. 2007), only 2 km northwest from the initial location by the Central

Weather Bureau (CWB) (Fig. 2a). Hypocentral depth given by different sources [Broadband Array in Taiwan for Seismology (BATS: http://bats.earth.sinica.edu.tw); CWB; Kuochen et al. 2007] varies from 10 to 25 km (Table 1). Published focal mechanisms for this event indicate a reverse fault motion with a possible small strike-slip component [BATS, Global Centroïd Moment Tensor catalogue (GCMTc), Kuochen et al. 2007, Fig. 2a].

The causative fault of the Chengkung earthquake, the 20-km-long Chihshang fault, is the best-known segment of the LVF (Ho 1986; Tsai 1986). This is an east-dipping thrust segment striking $N20°E$; the slip is oblique with a minor left-lateral component (Angelier et al. 1997, 2000; Chang et al. 2000; Lee et al. 2006). At its southernmost tip, the Chihshang fault splits into two branches separated by the Pinanshan Massif, the Lichi Fault on its East side and the Luyeh Fault on its west side (Yu & Kuo 2001, Figs 1b and 2a). Mainly based on geodetic triangulation surveys, Lee et al. (1998) demonstrated that the Lichi Fault is a left-lateral strike-slip fault whereas the Luyeh Fault is pure reverse, and explained this

<footnote n="1"><sup>1</sup>Geosciences Azur, Observatoire Océanologique de Villefranche, BP. 48, 06235, Villefranche-sur-mer and CNRS-University of Nice Sophia-Antipolis–250, 7 Rue Albert Einstein Sophia Antipolis, 06560 Valbonne, France. E-mail: mozzi@geoazur.obs-vlfr.fr</footnote>

<footnote n="2"><sup>2</sup>International Associated Laboratory ADEPT France-Taiwan, CNRS-NSC</footnote>

<footnote n="3"><sup>3</sup>Institut Universitaire de France, France</footnote>

<footnote n="4"><sup>4</sup>Department of Geosciences, National Taiwan University, Taipei, Taiwan</footnote>

<footnote n="5"><sup>5</sup>Institute of Earth Sciences, Academia Sinica, P.O. Box 1-55, Nankang, Taipei, Taiwan</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** (a) Tectonic setting of Taiwan with the convergence of the Philippine Sea Plate and the Eurasian Plate. Open arrow: convergence direction and velocity from Yu *et al.* (1997). Dashed rectangle: study area. LVF: Longitudinal Valley Fault. (b) Perspective view of study area. Black star: epicentre of the 2003 Chengkung earthquake. The Central Range (CER) is separated from the Coastal Range (COR) by the Longitudinal Valley (LV). The southern segment of the Longitudinal Valley Fault is the Chihshang Fault (CF). PM: Pinanshan Massif; LiF: the Lichi Fault; LuF: the Luyeh Fault and YF: the Yuli Fault.</figcaption></figure>

difference in behaviour between these parallel active faults in terms of strain partitioning due to the obliquity of plate collision with respect to the LVF. Farther North, the Yuli Fault runs parallel to the LVF. It occupies a more central position within the Longitudinal Valley (Fig. 1b) and exhibits a left lateral movement (Shyu *et al.* 2007).

Both the Yuli Fault and the Chihshang Fault were activated during one of the most destructive seismic crises of the 20th century in Taiwan that affected the entire Longitudinal Valley on 1951 November. This seismic sequence was dominated by two main shocks with magnitudes greater than 6.5 (Hsu 1962; Shyu *et al.* 2007). After the 1951 major crisis, the Chihshang Fault, that we are interested in, did not generate large earthquakes until the 2003 Chengkung earthquake sequence. However, fast surface creep had been observed since 1986 (Angelier *et al.* 1997, 2000; Lee *et al.* 2006; Hu *et al.* 2007). From the background seismicity in the Chihshang Fault area between 1998 and 2001, a clear listric geometry was evidenced in cross-section, with the fault dip varying from 60° to 70° at depths shallower that 10 km to less than 40° at 25–30 km depth (Chen & Rau 2002; Kuochen *et al.* 2004).

The link between the 2003 Chengkung earthquake and the Chihshang fault is attested by the location of the main surface ruptures that took place along the pre-existing creeping segment of the fault (Lee *et al.* 2004), also by the listric shape of the after-shocks distribution in cross-section (Kuochen *et al.* 2007) (Fig. 2b), and by the location of the epicentre itself, about 12 km east of the Chihshang Fault.

Previous studies of the Chengkung event determined the coseismic displacement from the geodetic (GPS) data or from the static part of strong motion records (Wu *et al.* 2006; Ching *et al.* 2007; Hu *et al.* 2007). These analyses produced a fixed image of the global slip distribution and could not illustrate the evolution of the rupture with time. In this paper, we aim at combining the resolving power of the GPS data and of the seismological time-series (strong motion and teleseismic) to produce a joint inversion modelling of the coseismic rupture event. We thus intend to better constrain the distribution of the coseismic slip, and through consideration of time-series, to establish the rupture evolution throughout time.

## 2 GEODETIC AND SEISMIC DATA

Coseismic GPS data is issued from 89 stations located at distance between 3 and 62 km from the Chengkung earthquake epicentre. The GPS data processing has been described by Chen *et al.* (2006). The resulting horizontal displacement vectors (Fig. 3, left-hand panels) display two symmetrical fans with respect to the Chihshang Fault. Remarkably, the vertical component (Fig. 3, right-hand panels) of all stations within the two blocks recorded an uplift (Chen *et al.* 2006; Wu *et al.* 2006; Hu *et al.* 2007). The maximum displacement is located in the hanging wall (in the Coastal Range) near the epicentre with 126 and 263 mm for the horizontal and vertical component, respectively (Chen *et al.* 2006). However, we should weight some GPS stations according to two criteria. First, less weight was given to stations where measurement uncertainties are large. Second, we
<figure><img src="image_3.png" /><figcaption>**Figure 2.** (a) Study area with main faults (black lines). Epicentre location (black circle) and focal mechanisms from different sources (Kuochen et al. (2007), BATS and GCMTc). Grey open dots: aftershocks of the Chengkung earthquake. The cross-section location (line AB) is perpendicular to the Chihshang Fault. Earthquake focal mechanisms shown as usual ‘beachball’ stereoplots, with tension and pressure dihedral, respectively black and white. (b) Cross-section AB of the 2003 Chengkung earthquake aftershocks (grey circles). Same focal mechanisms as in (a) top but shown in lateral projection.</figcaption></figure>

considered some structural peculiarities that could not be taken into account with the fault model used in the inversion. These structural complex zones are located in the northern part of the study area, where the coseismic process induced some deformation on the Yuli fault, and also in its southern part, where the LVF splits into two faults so that strain partitioning is expected to occur (Figs 1b and 2a). In these two areas, GPS data are downweighted by 50 per cent. In addition stations located in the Central Range have been affected by a more or less continuous uplift that contaminates the coseismic data. Vertical components from Central Range stations are hence not taken into consideration.

Accelerometric data were recorded at five stations of the CWB network located at distances between 3 and 34 km from the epicentre of the Chengkung earthquake (Fig. 4a). Stations located near the surface rupture of the Chihshang fault were removed from the analysis because their waveforms, likely to be affected by travel paths within the fault zone, could not be reproduced properly. Seismograms, initially in acceleration, were processed as follows. They were integrated twice in order to get displacement and a bandpass filter was applied between 0.05 and 0.25 Hz. The low cut-off frequency (0.05 Hz) was chosen in order to reduce the influence of the long-period noise that results from the base line shift in acceleration and its double integration. The high cut-off frequency (0.25 Hz) was defined to limit high frequencies related to waveform complexities that could not be reproduced with the use of simple 1-D velocity models. For each station, a specific 1-D model was established by using a grid search that explored more than a thousand different models with varying P-wave velocity ($Vp$) and layer thickness as in Örgülü et al. (2005). The $Vp/Vs$ ratio was adjusted for each station and a unique value was allowed for each model. In this preliminary analysis the rupture was assumed to be circular, with a constant rupture velocity on a simple rectangular fault plane. The criterion for selecting the optimal velocity model was the improvement in waveform modelling. Fig. 4(b) displays the resulting velocity model for each station. Only two stations, both located south of the Chengkung Earthquake area, slightly differ from a half-space model by the addition of a thin low velocity layer in the uppermost three kilometres. In particular, TTN046 is clearly located in an intramountain basin, the Taiyan basin (Figs 4a and b). $Vp/Vs$ ratios vary from 1.74 to 1.91. Stations requiring a high ratio are all located in the Coastal Range. Such a variation in the $Vp/Vs$ ratio was also pointed out by Kim et al. (2006) who related this change to the plate boundary between the continental crust (Central Range), characterized by a low ratio (Christensen 1996) and the oceanic crust (Coastal Range), characterized by both high ratio and high $Vp$ (Fowler 1990).

The teleseismic data are retrieved from the IRIS data centre (http://www.iris.edu/cgi-bin/wilberII_page1.pl) and correspond to 17 broad-band stations located worldwide, at angular distances between 30° and 90° from the main shock epicentre (Fig. 4c). The records were deconvoluted from the instrumental response,

<table><caption>**Table 1.** 2003 Chengkung earthquake parameters from different sources.</caption><thead><tr><th>Source</th><th>Epicentre location</th><th>Hypocentre depth</th><th>M<sub>w</sub></th><th>Focal mechanism (strike/dip/rake)</th></tr></thead><tbody><tr><td>Kuochen et al. (2007)</td><td>121.324°E / 23.106°N</td><td>16 km</td><td>6.8</td><td>37/50/94</td></tr><tr><td>CWB</td><td>121.34°E / 23.10°N</td><td>10 km</td><td>M<sub>L</sub> = 6.6</td><td></td></tr><tr><td>BATS</td><td>CWB</td><td>20 km</td><td>6.5</td><td>23/42/104</td></tr><tr><td>GCMTc</td><td>121.43°E / 22.94°N</td><td>25 km</td><td>6.8</td><td>10/51/69</td></tr><tr><td>This study</td><td></td><td>18 km</td><td>6.8</td><td>24/45/74</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** GPS data from Chen *et al.* (2006) On left: map of horizontal component (open arrow) of the coseismic displacement from GPS data. Line AB, CD and EF in rectangles: location of thick cross-sections, used to show the vertical component. Black star: Chengkung earthquake epicentre. Faults as thick black lines. On right: Vertical component of the GPS data (black dots) projected along lines AB, CD and EF. Small open rectangles: intersection between the cross-section and the Longitudinal Valley Fault (LVF).</figcaption></figure>

integrated to displacement and equalized to the same epicentral
distance (40°), according to Nabelek's method (1984).

## 3 COSEISMIC SLIP DISTRIBUTION: THE METHOD

In order to recover the coseismic slip distribution, we performed separate and simultaneous inversions of the three data subsets (teleseismic, strong-motion and GPS) following the approach described by Delouis et al. (2002). In this method, the fault plane is discretized in subfaults of fixed dimensions, and a simulated annealing algorithm is used to determine the inversion free parameters, which are for each subfault:

(1) The rake, allowed to vary within the interval 65° ± 35°,

(2) The amplitude of two mutually overlapping isosceles triangular functions (with a fixed width of 1.5 s) defining the subfault source time function,

(3) The rupture onset time that may vary according to two bounding rupture velocities 1.0 and 3.4 km s⁻¹, respectively.

Throughout this paper, rupture velocity refers to the average rupture velocity from the hypocentre, not the local rupture velocity.

The convergence criterion adopted in the inversion method is based on the minimization of a cost function defined as the sum of the rms (root mean square) misfits between observed and computed data and a function minimizing the total seismic moment. Syn-

thetic seismograms at strong motion stations are computed using the discrete wavenumber method of Bouchon (1981) which built all the waveform at regional distance in a given range of frequencies (between 0.05 and 0.25 Hz in this case) inside a stratified 1-D medium. At teleseismic distances, the ray-theory approximation of the approach of Nabelek (1984) is used to compute synthetics in the range of frequencies between 0.01 and 0.8 Hz for the *P* wave and between 0.01 and 0.4 Hz for the *SH* wave. For GPS data, the synthetic near-field static displacement is generated using the formulation of Savage (1980) that considers each subfault as a slip surface embedded in an elastic half-space. In our final result, all data sets are equally weighted.

To model the coseismic slip distribution of the Chengkung earthquake, the fault geometry used was retrieved from the aftershock seismicity at depth. Since the aftershocks display a listric shape, a single planar surface could not be adopted. Our fault model consists of four segments with identical strike but decreasing dip as depth increases. We checked different dip and width for these segments to obtain the best geometrical fit with the aftershocks distribution. At the surface, the shallower segment coincides with the surface trace of the Chihshang fault and with the coseismic surface breaks (Fig. 5). The four segments were subdivided into a total of 120 subfaults of equal length but with different widths depending on the fault segment (Table 2). As the uppermost segment never displayed any significant slip in all preliminary inversions, a relatively large width was allocated to this segment. This absence of large slip was
<figure><img src="image_5.png" /><figcaption>**Figure 4.** (a) Map of the accelerometric stations (black triangles) used in this study. Black star: Chengkung earthquake epicentre. (b) Velocity models selected for each accelerometric stations. H (km): layer thickness. (c) Map of the teleseismic stations (black triangles) between 30° and 90° from the epicentre (black star).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** 3-D view of the model adopted for fault geometry. The four segments (grey shaded rectangles) best fit the distribution of aftershocks (black dots) of the Chengkung earthquake at depth, as well as the fault traces at the surface (grey lines). Black stars: Chengkung earthquake epicentre and hypocentre (linked by dashed line). Open triangles with references: accelerometric stations.</figcaption></figure>
<table><caption>**Table 2.** Parameters for the four segments of the fault model.</caption><thead><tr><th>Segment</th><th>Strike</th><th>Dip</th><th>Rake</th><th>l<sub>tot</sub></th><th>w<sub>tot</sub></th><th>dl</th><th>dw</th><th>Centre depth (km)</th></tr></thead><tbody><tr><td>1</td><td>24</td><td>70</td><td>65 ± 35</td><td>48</td><td>5.1</td><td>4</td><td>5.1</td><td>2.5</td></tr><tr><td>2</td><td>24</td><td>70</td><td>65 ± 35</td><td>48</td><td>7.5</td><td>4</td><td>3.7</td><td>8.5</td></tr><tr><td>3</td><td>24</td><td>45</td><td>65 ± 35</td><td>48</td><td>11.3</td><td>4</td><td>3.8</td><td>16</td></tr><tr><td>4</td><td>24</td><td>30</td><td>65 ± 35</td><td>48</td><td>28.0</td><td>4</td><td>7.0</td><td>27</td></tr></tbody></table>

*Notes*: Strike/Dip/Rake: focal mechanism parameters. *l*<sub>tot</sub>, *w*<sub>tot</sub>: length and width in km. Each segment is subdivided into subfaults of length (*dl*) and width (*dw*) in km. Centre depth: depth (km) of the centre of each segment.

confirmed by our final inversion and also by the field data (Lee *et al.* 2006) and therefore it was unnecessary to discretize it finely. A rupture initiation with a hypocentre at 18 km depth showed the best agreement with all the data.

## 4 RESOLUTION TESTS

In order to validate our results, we carried out two resolution tests with increasing complexity.

In a first resolution test (Fig. 6), the input model is a simple rectangular slip patch located at the centre of the fault model

(Fig. 6a). The rupture initiates and propagates from the hypocentre with a constant velocity rupture (*V*<sub>r</sub> = 2.5 km s<sup>-1</sup>). From this model, seismological and geodetic synthetic data are generated, and then taken as input data for separate and joint inversions. The results indicate that none of the separate inversions (Figs 6c–e) are able to fully recover the input model. We note, however, that the best result issued from the separate inversions was obtained with the strong motion data using only five stations. The joint inversion (Fig. 6b) displays the best-reconstructed slip map, indicating that a complete modelling should involve the three sources of data (GPS, accelerometric and teleseismic). This requirement partly results from the distribution of stations (e.g. the absence of GPS and accelerometric stations east off the Coastal Range). The teleseismic stations help filling the gap left by the onland GPS and accelerometric stations, but tend to spread out the slip map. This first test showed that the central part of the fault model was well resolved by the joint.

In a second resolution test, the input model is composed of three slip patches, each one characterized by a specific rupture velocity (*V*<sub>r1</sub>, *V*<sub>r2</sub>, *V*<sub>r2</sub>, Figs 7a left-hand panel and c left-hand panel). Only the result of the joint inversion is discussed herein. The three patches are recovered at the proper depth (Fig. 7a, on right-hand panel).

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Fault slip model, first resolution test corresponding to a rectangular slip zone (outlined by a thick dashed line). White arrows: slip vectors; white triangle: hypocentre. (a) Synthetic model with a constant rupture velocity for all calculation points. (b) Final joint inversion. (c)–(e) Separate data set inversions.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Fault slip model, second resolution test with three patches of constant displacement at all points and dashed lines outlining the patch contour. (a) Synthetic slip model (left-hand panel) and inversion from the joint inversion of teleseismic, strong motion and GPS data on the right. The open arrows show the slip vectors. (b) Rupture timing: hypocentral distance of rupture onset as a function of time. Open and grey symbols for synthetic model and inversion, respectively (one symbol for each subfault). Subfaults from each inversion define three subgroups (circles, square and triangle) as in the synthetic model. A linear regression was performed for each subgroup (straight lines) to determine the rupture velocity ($V_{r1}$, $V_{r2}$ and $V_{r3}$). The regression coefficient $R^2$ is indicated for each case. (c) Model (left-hand panel) and inversion (right-hand panel) slip map with corresponding $V_r$ ($V_{r1}$, $V_{r2}$ and $V_{r3}$) by subfaults. Each patch is characterized by a constant rupture velocity in the model. The resulting velocities from the joint inversion are shown with the same symbols, illustrating the good fit between velocity and location. (d) Model source time function (black line) compared to the computed one (dashed grey line). Insert: source time function obtained with the real data, for comparison.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Slip inversion results with real data (compare with synthetic experiments of Figs 6 and 7). Open triangle: hypocentre. Black dots: point sources on the fault model. (a)–(c): results for separate teleseismic, strong motion and GPS data inversions, respectively. (d): result for joint inversion using the three data sets. (e): overall source time function from the joint inversion (seismic moment rate as a function of time).</figcaption></figure>

The uppermost slipping zone, closest to the GPS and strong-motion stations, is particularly well reconstructed. The slip zone comprising the hypocentre is retrieved, but a certain degree of spreading remains in its deep part. The third patch to the SSW is found more widespread and its maximum slip is mislocated near the model edge. The portion of our fault model above 26 km depth is therefore well resolved, except for the SSW edge, with varying velocity rupture. The source time function (STF hereafter) corresponding to the release of the seismic moment rate with time, comprises two pulses, as does the real data joint inversion (Fig. 7d and insert). We investigated how rupture velocity ($V_r$) could be resolved. Rupture velocities are determined by linear regression on the hypocentral distance of the rupture onset versus time (Fig. 7b). In the synthetic model, the rupture extends from the hypocentre with $V_{r1}$ (3.2 km s⁻¹) in the NNE slip patch, corresponding to the first pulse of the STF (Fig. 7d). After a time gap of 5 s, the two other patches slip at the same time with different rupture velocities ($V_{r2} = 2.2$ km s⁻¹ for the SSW zone, $V_{r3} = 1.5$ km s⁻¹ for the upper patch) leading to the second pulse of moment release rate. This rupture scenario is well recovered in space (Fig. 7c) and time (Figs 7b and d) by the joint inversion. The error on $V_r$, as defined as the difference between the value in the synthetic model and the inverted one, does not exceed 0.1 km s⁻¹, even within the SSW patch where we

previously noticed a decrease in resolution. This test thus shows that we can resolve variations of the rupture velocity occurring during the rupture process with a rather good accuracy.

Since the resolution tests were performed without noise, they represent the optimal resolution, which could be achieved. We also performed tests with up to 15 per cent of additional random noise incorporated in the datasets, and the results remained stable.

## 5 COSEISMIC SLIP DISTRIBUTION: THE RESULTS

From the three real data sets, we conducted separate and joint inversions. For both the separate seismological inversions (teleseismic and strong-motion, Figs 8a and b), we obtained a relatively linear and narrow slip zone propagating from the hypocentral area to the SSW, composed of two main patches more or less connected. In the strong motion case the slip pattern is more oblique and reaches larger depths in the SSW. This tendency is still clearer in the GPS case where the two main patches became completely connected (Fig. 8c). The rupture area (herein defined as the slipping zone that underwent at least 20 cm of slip) reaches 650 km² according to the strong-motion and 705 km² from teleseismic data (Table 3), with an average slip of 70 cm in both cases. The rupture area from GPS is
<table><caption>**Table 3.** Main results of the separate and joint inversions for the Chengkung earthquake.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Separate inversion</th><th rowspan="2">Joint inversion</th></tr><tr><th>GPS</th><th>Strong-motion</th><th>Teleseismic</th></tr></thead><tbody><tr><td>Rupture area (km²)</td><td>597</td><td>650</td><td>705</td><td>812</td></tr><tr><td>Average slip (cm)</td><td>93</td><td>66</td><td>69</td><td>71</td></tr><tr><td>Maximal slip (cm)</td><td>193</td><td>184</td><td>191</td><td>195</td></tr><tr><td>Seismic moment (dyne.m)</td><td>1.98 × 10<sup>26</sup></td><td>1.68 × 10<sup>26</sup></td><td>1.93 × 10<sup>26</sup></td><td>2.07 × 10<sup>26</sup></td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Comparison between observed and computed data, for strong motion and teleseismic waveforms. Observed seismograms in black, computed in grey. Teleseismic stations are plotted over the focal mechanism (*P* and *SH*) given by our slip model at the hypocentre (strike/dip/rake = 24/45/74).</figcaption></figure>

597 km² (Table 3), slightly smaller than with seismological data but the average slip value is higher, 93 cm. The total seismic moment (*M*₀) was computed as the sum of the seismic moments of all the individual subfaults, using a unique value of rigidity (μ) of 3.2 × 10¹¹ N m⁻²:

$$M_0 = \mu \sum_i S_i \Delta u_i \quad (1)$$

with *S*ᵢ and Δ*u*ᵢ the area and the slip of the ith subfault, respectively.

The joint inversion of the three data sets (Fig. 8d) displays the two main slip patches observed with seismological data more or less connected but with a global shape similar to the GPS results. The average slip is 71 cm for a seismic moment of 2.07 × 10²⁶ N.m on a total rupture area of 812 km² (Table 3). The main first slip patch starts at the hypocentre and propagates to the South between

12 and 23 km. The second slip patch is located more to the South and at greater depth (20–27 km). A third slip zone with smaller slip values is observed in the far South at 17 km depth, with a branch rising towards the surface (up to 5 km depth) in the central part of the fault model. The possibility to model properly the three data sets simultaneously highlights the efficiency of the joint inversion (Figs 9 and 10).

Although filtered below 0.8 Hz the teleseismic *P* waves, as well as their modelling are dominated by frequencies below 0.3 Hz. We observe that the *SH* component of the teleseismic waves is less well resolved, especially for stations located South of the epicentre, whereas all the corresponding *P* waves are correctly modellized. This particular problem has also been described by Marson-Pidgeon *et al.* (2000) for very similar source-stations path. Their study was realized in the same conditions as our, with stations located between
Joint inversion data fit

<figure><img src="image_11.png" /><figcaption>Figure 10. Comparison between observed GPS data and computed horizontal and vertical displacements of the Chengkung earthquake. See also Fig. 3.</figcaption></figure>

30° and 90° from the epicentre for an earthquake located in Japan near to Taiwan and a filtering similar to ours. The authors considered that the ray path to those stations crossing several subduction zones is at the origin of the perturbations of the signal. This waveform distortion due to a plunging structure such as a subduction zone, was demonstrated for the SH component of the teleseismic waves by Igel & Ita (1997), and is likely to be the source of our modelling difficulties.

We investigated the range of variability of each parameter of the fault model. The simulated annealing algorithm requires a starting model to initiate the cost function. It also requires a seed number to initiate the random generator. The convergence path is expected to be influenced by both the starting model and the seed number. In a first step we performed 16 inversions with the same data set and algorithm but with different initial conditions as was done in Liu et al. (2006). For 10 models, the starting model is changed by varying the initial value of each parameter within its allowed boundary range. For the next six models, we used a different seed for the random generation. These 16 inversions provided a total of 1105 728 models, which allowed us to compute the standard deviation (σᵢ) for each parameter with respect to our best model in

a similar way as Piatanesi et al. (2007):

$$ \sigma_i = \sqrt{\frac{\sum_{j=1}^{j=n} \frac{(m_{ij} - m_i)^2}{E_j}}{\sum_{j=1}^{j=n} \frac{1}{E_j}}}, \quad (2) $$

where mᵢⱼ is the ith parameter of the jth model, mᵢ the ith parameter of the best model and Eⱼ the cost function of the jth model. The results are shown in Fig. 11. The standard deviation on slip amplitude (Fig. 11b) correlates with the slip values (Fig. 11a). For all subfaults where slip is larger than 60 cm, the standard deviation represents less than 30 per cent of the slip value. The rake deviation (Fig. 11c) is larger in the shallow part of the fault model.

We consider afterwards the evolution of the rupture process throughout time deduced from the joint inversion. The STF (Fig. 8e) displays two main pulses of seismic moment rate that may correspond to at least two stages in the rupture time history. More precisely, the first pulse is about two times larger than the second one and indicates a first stage with larger slip values. In order to better establish the rupture timing, we determined how rupture velocity evolves as a function of time (Fig. 12) by plotting
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Standard deviation on slip and rake. (a) Slip map from the joint inversion, the 60 cm slip boundary is underlined (thick line). Open triangle: hypocentre. (b) The slip deviation (cm) correlates the slip values [by comparison with (a)]. (c) Note that the rake deviation (°) is larger for the shallower part of the fault plane.</figcaption></figure>

the hypocentral distance of rupture initiation versus time for sub-faults that slipped more than 20 cm (Fig. 12a). This graph displays a clear deceleration of the rupture in three stages. For the three time intervals, $V_r$ is determined by linear regression (Fig. 12a) and the corresponding slip maps are presented in Fig. 12(b). To investigate the range of variability for $V_r$, we used the same 16 inversions with varying initial conditions as described above. We extracted the best solution of each inversion and built a composite graph of hypocentral distance versus time in Fig. 12(c). For each of those 16 models we divided the subfaults into four groups within the three same time intervals as before. For each group a linear regression was performed. From these regressions we deduced the different $V_r$ values and their corresponding regression coefficient $R^2$ to finally compute the standard deviation on $V_r$ as:

$$ \sigma_i = \sqrt{\frac{\sum_{j=1}^{j=16} R_j^2 [V_{r_{ij}} - V_{r_i}]^2}{\sum_{j=1}^{j=16} R_j^2}} \quad (3) $$

with $V_{r_{ij}}$ the ith $V_r$ of the jth model and $V_{r_i}$ the ith $V_r$ of our best-fitting model. Fig. 12(c) displays the regression lines for the 16 graphs with the associated standard deviation on $V_r$ that never exceeds 10 per cent of the $V_r$ value.

From the resolution tests, we know that it is possible to discriminate between different rupture velocities occurring simultaneously, and we can infer a rupture scenario. From 0 to 5.5 s after the rupture initiation, the rupture propagates fast, at $3.4 \pm 0.3$ km s$^{-1}$ ($V_{r1}$) along the two main slip patches. Then, the rupture velocity decreases to $2.4 \pm 0.1$ km s$^{-1}$ ($V_{r2}$) and the rupture propagates upwards. Finally, after 15 s, the rupture propagation slows down to $1.6 \pm 0.1$ km s$^{-1}$ ($V_{r3}$ up) for a group or subfaults located in the upper part of the rupture zone, and down to $1.1 \pm 0.1$ km s$^{-1}$ ($V_{r3}$ down) for the subfaults located in the lower southern part of the slipping zone.

In all the inversions previously shown, the three data sets were equally weighted. In order to assess the impact of the data weight on the inversion results we have tested six different combinations of weights individually varying between 0.6 and 1.4, with a fixed sum, 3, for the three data weights (Fig. 13). Larger variation of weights approximates the individual data set inversions. In the six tested cases presented in Fig. 13, the main results display two main

slip patches and a secondary one, corresponding to a rupture scenario that always includes three time intervals with four distinct velocity rupture. As Fig. 13 shows, the velocity rupture decrease as time increases. Incidentally, the inversion behaviour throughout weight variation remains constant and supports our choice to present equally weighted inversions only.

## 6 DISCUSSION AND CONCLUSION

In this study we obtained the coseismic slip distribution of the 2003 Chengkung earthquake by combining dense GPS measurements with near-source strong motion and far-distance broad-band records. By modelling the waveforms, we could determine the rupture time history, in particular variations of the rupture velocity, an aspect not addressed in previous studies. It is worth noting that the use of three major sources of data (teleseismic, accelerometric and geodetic, the latter involving both the horizontal and vertical displacements) provides better constraint on the fault slip pattern of the Chengkung earthquake than the use of separate sources. The joint inversion, taking these three complementary data sources, is therefore crucial.

From the hypocentre at 18 km depth, the rupture propagated laterally towards the SSW for about 30 km. The rupture grew also in the down- and up-dip directions, with most of the slip located between 12 and 26 km depth. Two slip patches dominate, with slip values exceeding 1.5 m (Fig. 8d). The absence of large slip near the surface, as was effectively observed in the field (Lee et al. 2006), is attested by the resolutions tests, which indicate that slip, or the absence of slip, is well constrained in the upper part of the fault model.

The rupture process is characterized by a clear deceleration with time (Fig. 12a). The main two slip patches ruptured during the first 5s, with a fast rupture velocity $3.4$ km s$^{-1}$. During the next 10 s, rupture velocity lowered to $2.4$ km s$^{-1}$. Finally, for the next 10s, it slowed down to $1.6$ and $1.1$ km s$^{-1}$ (Fig. 12b). The whole rupture lasted 30s and the seismic moment was $2.07 \times 10^{26}$ dyne cm ($M_w$ 6.8).

In comparison to previous studies based on static displacements from GPS or strong motion records (Ching et al. 2007; Hu et al. 2007; Cheng et al. 2009), our fault model is slightly more
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Rupture timing of the Chengkung earthquake. (a) Hypocentral distance of subfaults as a function of rupture onset time. In three time steps, four groups of subfaults are differentiated, depending on rupture velocity, each group represented by a specific symbol (black circles, grey squares, triangles and open triangles). Black lines: linear regression for each group. Corresponding rupture velocities ($V_{r1}$, $V_{r2}$, $V_{r3up}$ and $V_{r3down}$) and regression coefficients ($R^2$) are also shown. (b) Contribution of time interval to the final slip distribution. Open triangle: hypocentre. Dashed vertical line: separation between time periods. (c) Hypocentral distance of subfaults as a function of rupture onset time for the 16 best models obtained in varying the initial conditions (see text) with corresponding regression lines. For each $V_t$, the standard deviation $\sigma$ is shown.</figcaption></figure>

concentrated in space, with higher peak slip values. It is usual that inversions based on geodetic data alone result in smoother slip distributions. Higher frequency seismological data generally provide more details in the rupture model. Moreover, previous studies included a certain degree of smoothing in their results, which tends to expand the slip area, whereas our inversion did not incorporate a specific smoothing operator.

In the 2003 Chengkung Earthquake crisis the coseismic slip distribution is thus well determined, so that one can now focus on the initial aftershock sequence. During the first 2 d after the Chengkung earthquake, a strong decrease in the number of aftershocks per hour

is observed (Fig. 14a). To define the most active area we selected events with magnitude greater than 3.0 and located less than 3 km away from the fault plane (as defined earlier in Fig. 5). We have orthogonally projected the aftershocks foci onto the multi surface fault plane to easily compare their location with the coseismic slipping zone (Figs 14b and c). During the main shock day, most aftershocks occurred within the main shock slip zone, including the four events of $M_w$ greater than 5. The uncertainties on the aftershocks location, estimated to be a few kilometres, prevents us from determining if aftershocks occurred exactly on the coseismic slip maxima or in their immediate vicinity. However, a minor proportion of the
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Inversion results with varying data weights. The main characteristics of slip map and velocities rupture remain constant. In slip map the hypocentre is represented by the open triangle.</figcaption></figure>

aftershocks occurred outside the coseismic zone, to the South and in the deepest part of the fault. Regarding the day following the main shock, a similar seismic repartition was observed, with less events but with a large aftershock ($M_w = 5.4$) occurring in the southern part of the fault (Fig. 14c). Note that our model does not display any rupture at the location of this aftershock. In this subarea, the eastern dip of the fault is farther from the onshore GPS and strong motion stations and is thus less resolved. As a consequence a portion of fault slip may be missing in the deeper part of our coseismic model. As another hypothesis, the post-seismic afterslip below the coseismic rupture zone may have triggered these deep aftershocks. In any case, we have verified with a specific resolution test that slip is not well constrained in this area.

Ching *et al.* (2007) and Cheng *et al.* (2009) found post-seismic slip to occur at shallow depth in the 3 months following the main shock, immediately above the coseismic patch that we identified

in the hypocentral area. This was confirmed by rapid post-seismic creep at the surface, exceeding the small coseismic surface slip (Lee *et al.* 2006). If post-seismic afterslip occurred in the deep part of the fault, at greater distance from the GPS stations, it could have been undetected, or hardly differentiated from shallower slip. Despite these side uncertainties, the 2003 Chengkung earthquake demonstrates that the Chihshang Fault behaves in a way which is not commonly documented for other reverse faults, with a coseismic rupture essentially confined at depth (below 10–12 km), and a shallow part of the fault accommodating deformation by rapid post-seismic slip and slow interseismic creep.

Our joint inversion modelling provides a detailed coseismic slip map for the 2003 Chengkung earthquake using realistic listric geometry for the causative fault. It also reveals how the rupture has propagated along the fault surface. Understanding the relation between the coseismic behaviour and the post-seismic slip of the
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Aftershock distribution and relation with the coseismic slip. (a) Number of aftershocks per interval of one hour for the first 2 d after the main shock. (b) Fault model and aftershocks (dots) in cross-section. Grey dots: aftershocks located less than 3 km from the fault model. (c) Surface projection of the fault model, slip map and aftershocks (white squares) of magnitude ($M_w$) larger than 3.0 located near the fault (distance <3 km) for the same 2 d as in (a) and orthogonally projected on the fault model. The rupture area of the Chengkung Earthquake is outlined by a dashed line (slip larger than 20 cm). Black stars: aftershocks of $M_w > 5$. White star: Chengkung Earthquake epicentre.</figcaption></figure>

southern segment of the LVF, one of the most active faults of the island of Taiwan, can help for seismic hazard assessment.

## ACKNOWLEDGMENTS

This work was done in the frame of the cooperation in Earth Sciences between France and Taiwan, encouraged by the French CNRS

and the National Science Council of Taiwan (International laboratory ADEPT), and also the Nice University, the Institut Français à Taipei (IFT) and the Bureau de Représentation de Taipei (BRT) in Paris. The first author received a PhD grant from the French Ministry. The constructive comments of two anonymous referees and the Editor allowed us to significantly clarify and improve the manuscript.
REFERENCES

Angelier, J., Chu, H.T. & Lee, J.C., 1997. Shear concentration in a colli-
sion zone: kinematics of the active Chihshang fault, Longitudinal Valley,
eastern Taiwan, *Tectonophysics*, **274**, 117–144.

Angelier, J., Chu, H.T., Lee, J.C. & Hu, J.C., 2000. Active faulting and earthquake hazard: the case study of continuous monitoring of the of the Chihshang Fault, Taiwan, *J. Geodyn.*, **29**, 151–185.

Bouchon, M., 1981. A simple method to calculate Green's functions for elastic layered Media, *Bull. seis. Soc. Am.*, **71**, 959-971.

Chang, C.P., Angelier, J. & Huang, C.Y., 2000. Origin and evolution of a
melange: the active plate boundary and suture zone of the Longitudinal
Valley, Taiwan, *Tectonophysics*, **325**, 43–62.

Chen, H.H. & Rau, R.J., 2002. Earthquake locations and style of faulting in a active arc-continent plate boundary: the Chihshang fault of eastern Taiwan, in AGU Fall Meeting, San Francisco.

Chen, H.Y., Yu, S.B., Kuo, L.C. & Liu, C.C., 2006. Coseismic and postseismic surface displacements of the 10 December 2003 (Mw 6.5) Chengkung, eastern Taiwan, earthquake, *Earth Planets Space*, **58**, 5–21.

Cheng, L.W., Lee, J.C., Hu, J.C. & Chen, H.Y., 2009. Coseismic and Post-seismic Slip Distribution of the 2003 Mw=6.5 Chengkung Earthquake in Eastern Taiwan: elastic modeling from inversion of GPS data, *Tectono-physics*, doi:10.1016/j.tecto.2007.11.021.

Ching, K.-E., Rau, R.-J. & Zeng, Y., 2007. Coseismic source model of the 2003 Mw 6.8 Chengkung earthquake, Taiwan, determined from GPS mea-surements, *J. geophys. Res.*, **112**, B06422, doi:10.1029/2006JB004439.

Christensen, N.I., 1996. Poisson's ratio and crustal seismology, *J. geophys. Res.*, **101**, 3139–3159.

Delouis, B., Giardini, D., Lundgren, P. & Salichon, J., 2002. Joint inversion of InSAR, GPS, teleseismic, and strong-motion data for the spatial and temporal distribution of earthquake slip: application to the 1999 Izmit mainshock, *Bull. seism. Soc. Am.*, **92**, 278–299.

Fowler, C.M.R., 1990. *The Solid Earth: An Introduction to Global Geophysics*, p. 472, Cambridge University Press, New York.

Ho, C.S., 1986. A synthesis of the geologic evolution of Taiwan, *Tectono-physics*, **125**, 1–16.

Hsu, T.L., 1962. Recent faulting in the Longitudinal Valley of eastern Taiwan, *Mem. Geol. Soc. China*, **1**, 95–102.

Hu, J.C. et al., 2007. Coseismic deformation revealed by inversion of strong motion and GPS data: the 2003 Chengkung earthquake in eastern Taiwan, *Geophys. J. Int.*, **169**, 667–674.

Igel, H. & Ita, J., 1997. The effects of subduction zones on teleseismic SH waves: a numerical study, in *Upper Mantle Heterogeneities from Active and Passive Seismology*, pp. 333–341, ed. Fuchs, K., Kluwer Academic Publishers, Moscow.

Kim, K.H., Chiu, J.M., Pujol, J. & Chen, K.C., 2006. Polarity reversal of active plate boundary and elevated Oceanic Upper Mantle beneath the collision suture in Central Eastern Taiwan, *Bull. seism. Soc. Am.*, **96**, 796–806, doi:710.1785/0120050106.

Kuochen, H., Wu, Y.M., Chang, C.H., Hu, J.C. & Chen, W.S., 2004. Relo-
cation of eastern Taiwan earthquakes and tectonic implications, *Terres.
Atmos. Ocean. Sci.*, **15**, 647–666.

Kuochen, H., Wu, Y.M., Chen, Y.G. & Chen, R.Y., 2007. Mw 6.8 Chengkung earthquake and its associated seismogenic structures, *J. Asian Earth Sci.*, **31**, 332–339.

Lee, J.C., Angelier, J., Chu, H.T., Yu, S.B. & Hu, J.C., 1998. Plate-boundary strain partitioning along the sinistral collision suture of the Philippine and Eurasian plates: analysis of geodetic data and geological observation in southeastern Taiwan, *Tectonics*, **17**, 859–871.

Lee, J.C. et al., 2004. The 2003 Mw 6.5 Chengkung Earthquake and seis-
motectonics: earthquakes along a Creeping Fault at the Plate Suture, in
*International Conference in Commemoration of 5th Anniversary of the
1999 Chi-Chi Earthquake*, Taiwan, Abstract volume.

Lee, J.C., Chu, H.T., Angelier, J., Hu, J.C., Chen, H.Y. & Yu, S.B., 2006. Quantitative analysis of surface coseismic faulting and postseismic creep accompanying the 2003, Mw = 6.5, Chengkung earthquake in eastern Taiwan, *J. geophys. Res.*, **111**, doi:10.1029/2005JB003612.

Liu, P., Custódio, S. & Archuleta, R., 2006. Kinematic inversion of the 2004 M 6.0 Parkfield Earthquake including an approximation to site effects, *Bull. seism. Soc. Am.*, **96**, S143–S158, doi:110.1785/0120050826.

Marson-Pidgeon, K., Kennet, B.L.N. & Cambridge, M., 2000. Source depth and mechanism inversion at teleseismic distances using a neighborhood algorithm, *Bull. seism. Soc. Am.*, **90**, 1369–1383.

Nabelek, J.L., 1984. Determination of earthquake source parameters from inversion of body waves, *PhD thesis*. Mass. Inst. of Technol., Cambridge.

Örgülü, G., Delouis, B., Huang, B.S. & Legrand, D., 2005. Discrimination of the fault plane by waveform modeling: a case study for moderate-sized earthquakes in Taiwan, *Bull. seism. Soc. Am.*, **95**, 1825–1840.

Piatanesi, A., Cirella, A., Spudich, P. & Cocco, M., 2007. A global search inversion for earthquake kinematic rupture history: application to the 2000 western Tottori, Japan earthquake, *J. geophys. Res.*, **112**, B07314, doi:07310.01029/02006JB004821.

Savage, J.C., 1980. Dislocations in seismology, in *Dislocations in Solids*, ed. Navarro, F. R. N., Amsterdam, North-Holland.

Shyu, J.B.H., Chung, L.H., Chen, Y.G., Lee, J.C. & Sieh, K., 2007. Re-evaluation the surface ruptures of the November 1951 earthquake series in eastern Taiwan, and its neotectonic implications, *J. Asian Earth Sci.*, **31**, 317–331.

Teng, L.S., 1990. Tectonic evolution of late Cenozoic arc-continent collision in Taiwan, *AAPG Bull.*, **74**, 1004–1005.

Tsai, Y.B., 1986. Seismotectonics of Taiwan, *Tectonophysics*, **125**, 17–38.

Wu, Y.M. et al., 2006. Coseismic versus interseismic ground deforma-
tions, fault rupture inversion and segmentation revealed by 2003 Mw
6.8 Chengkung earthquake in eastern Taiwan, *Geophys. Res. Lett.*, **33**,
doi:10.1029/2005GL024711, L02312.

Yu, S.B., Chen, H.Y. & Kuo, L.C., 1997. Velocity field of GPS stations in the Taiwan area, *Tectonophysics*, **274**, 41–59.

Yu, S.B. & Kuo, L.C., 2001. Present-day crustal motion along the Longitu-
dinal Valley Fault, eastern Taiwan, *Tectonophysics*, **333**, 199–217.