# Kinematic and spontaneous rupture models of the 2005 Tarapacá intermediate depth earthquake

S. Peyrat, P. Favreau

► To cite this version:

S. Peyrat, P. Favreau. Kinematic and spontaneous rupture models of the 2005 Tarapacá intermediate depth earthquake. Geophysical Journal International, 2010, 181, pp.369-381. 10.1111/j.1365-246X.2009.04493.x. insu-03605332

HAL Id: insu-03605332

https://insu.hal.science/insu-03605332

Submitted on 11 Mar 2022

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 Attribution 4.0 International License
# Kinematic and spontaneous rupture models of the 2005 Tarapacá intermediate depth earthquake

S. Peyrat and P. Favreau

Seismology Laboratory, Institut de Physique du Globe de Paris (CNRS-UMR 7154), 4 Place Jussieu, 75252-Paris, France. E-mail: peyrat@ipgp.jussieu.fr

Accepted 2009 December 21. Received 2009 December 20; in original form 2008 April 1

## SUMMARY

The 2005 June 13, an intraplate, intermediate depth earthquake ($M_w$ = 7.8, 98 km depth) occurred in northern Chile. Previous studies show that this earthquake is a slab-pull event with down dip extensional source mechanism. However, the physical origin and the role in seismic cycle of this kind of event are still under debate. We present a seismological source study from strong motion data, based on the simplified slip patch kinematic parametrization. We find two distinct sources. The main source is almost centred around the hypocentre and triggers the secondary one, situated deeper in the slab. The weaknesses and uncertainties of the kinematic models are partly resolved and discussed by constructing dynamic spontaneous rupture models with two radiative asperities and one embedding zone of non-radiative slip process, a required feature to fulfil the seismic moment constraint. Two kinds of models are compatible with observations: a crack mode model (combined with fast rupture velocity and small zone of non-radiative slip) and a pulse mode model (combined with slower rupture velocity and large zone of non-radiative slip). The fast rupture model is preferred for its physical simplicity. Replacing our result in the geometry and thermal contexts of the subduction, we conclude that the rupture has probably broken the whole lithosphere by encountering little resistance and without being very sensitive to the presence of the double seismic zone. This suggests that the rupture happens along a weak plane, on which large slab instabilities can take place.

**Key words:** Earthquake dynamics; Earthquake ground motions; Earthquake source observations; Seismicity and tectonics; Subduction zone processes; Intraplate processes.

## 1 INTRODUCTION

Chile is a high seismic area and, in this subduction zone, several significant gaps are identified. On Fig. 1, we report the last large events and the NEIC seismicity. Between Arica and Antofagasta, the last megathrust earthquake occurred in 1877 around Iquique and produced a destructive tsunami. The northern and southern parts of this segment have been recently ruptured during the following earthquakes: 1995 Antofagasta ($M_w$ = 8.0), 2001 Arequipa ($M_w$ = 8.4) and 2007 Tocopilla ($M_w$ = 7.7). Therefore, no major shallow earthquake occurred in the northern region of Chile for more than 100 yr. Thus, assuming a strong coupling, this area is probably under a considerable tectonic loading.

The 2005 June 13, 22:44 GMT, an intraplate intermediate depth earthquake ($M_w$ = 7.8) occurred just east of this seismic gap. Hypocentre is located at -20.0486° (lat), -68.28° (lon) and at 98 km depth. This earthquake was followed by a large number of aftershocks, which is not so common for this kind of earthquakes (Astiz et al. 1988, show that most of the intermediate depth earthquakes have no aftershocks). Aftershocks distribution recorded by a temporary network was used to discriminate the real fault plane. Their distribution is confined within a subhorizontal area of about

60 km long, elongated in north-south direction and slightly dipping to the west. Hence the fault that ruptured is a nearly horizontal plane and the study of Peyrat et al. (2006) confirms that the event is a horizontal intraplate slab-pull rupture.

Intraplate intermediate depth earthquakes in subduction zones are not rare, but, although they are sometimes very damaging, detailed studies (Ide & Takeo 1996; Cocco et al. 1997; Mikumo et al. 2002) are not common. The 1939 Chilean earthquake ($M_w$ = 8.3) was a complex intermediate depth earthquake (Campos & Kausel 1990), presumably occurring on a vertical plane, and it was one of the most destructive events to ever occur in Chile (clearly more destructive than interplate events of similar or even larger size). This may be due partly to the intraplate nature of the event and its location below the central valley of south central Chile.

Most intermediate depth events in northern Chile are located in the depth range between 100 and 270 km, with two distinct groups separated by an aseismic gap. The focal depths for most of the intermediate depth earthquakes occur in the relatively narrow depth range of 100–150 km, and some events are located at depth between 200 and 270 km (Stauder 1973; Hanuš et al. 1999). The origin of these earthquakes at intermediate depth is still controversial. The causes and the mechanisms by
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Tectonic setting: location of rupture areas of the last large earthquakes in northern Chile, location of the main shock of the 2005 June 13 Tarapacá earthquake, and seismicity (from NEIC) with event of magnitude greater than 4.5 during the period 1973–2005. Accelerometric data: station locations and corresponding vertical accelerations.</figcaption></figure>

which these earthquakes nucleate and propagate remain poorly understood.

Several complementary explanations are proposed. First, the britt-
tle/unstable behaviour at this depth can be due to thermomechani-
cal properties, in presence of fluid, causing dehydration inside the
downgoing slab and/or reactivation of pre-existing faults. In gener-
ating earthquakes, dehydration embrittlement may also create new
faults or reactivate pre-existing weak zones (Kirby 1995; Peacock
2001; Jung et al. 2004). Second, one invokes the state of stress and
the geometry of the slab (pull, bending/unbending, lateral segmen-
tation and undulations) (e.g. Cocco et al. 1997). By studying the
1997 Punitaqui earthquake in central Chile, Lemoine et al. (2002)
conclude that intermediate depth earthquakes can be related to the
heterogeneity of the stress field inside bended downgoing slabs. Ide
& Takeo (1996) illustrated the role of the stratified state of stress in
the slab: they showed that the rupture of the Japanese 1993 Kushiro-
Oki earthquake (M<sub>w</sub> = 7.6), very similar to Tarapacá, was channeled

in the slab between the two layers of the double Wadati–Benioff
seismicity zone (DSZ). Third, the seismic cycle of these events,
situated near the bottom of the seismogenic zone, can be influenced
by the seismic cycle of the interplate thrust events and reversely.
Astiz et al. (1988) suggests that before a major thrust event, the
interplate boundary is strongly coupled, and the subducted slab is
under tension at intermediate depth. For the 1997 Punitaqui event,
Gardi et al. (2006) quantified the effect of static stress perturbations
due to an intermediate depth intraslab earthquake on seismic cycles
of large interplate earthquakes. These cycles may be totally linked
through stress transfers at the scale of the entire subduction zone.
The Tarapacá earthquake might be a clear evidence of this cycle
and could be compared to the 1950 intermediate depth earthquake
(M<sub>w</sub> ≈ 8.0) (Kausel & Campos 1992), that happened before the
1995 Antofagasta earthquake.

To test these models, or at least to determine some regional vari-
ations, one needs to know not only the global characteristics of the
event, but also enough details of the rupture to access slip heterogeneities and eventually local stress drop and strength excess, even frictional specificities. During the last decade, several attempts were undertaken to determine the spatial and temporal distribution of slip over the fault area of large earthquakes. These source tomographies have further been developed and applied mainly to large earthquakes in California and Japan, where near-source waveforms were collected during earthquakes. Past experiences show that to obtain a reliable image of the details of the source process, many conditions are required: a good density and coverage of the seismic array, a reasonable distance from the source, a good quality of records, a good knowledge of the velocity model and a good parametrization of the source model, that is, consistent with the intrinsic possibilities of the seismic array and meaningful for the source process itself (Hartzell et al. 2007). The benchmark of Mai et al. (2007) is very illustrative. For the Tarapacá 2005 event, the limited number of strong motion data prevails on us to use a simplified slip patch parametrization (Vallée & Bouchon 2004).

In Section 2, we review the results of previous studies and point out the distinctive features of the strong motion waveforms that motivate our study approach. In Section 3, we present the results of the simplified kinematic inversion. However different results with the same level of fit are obtained depending of the inversion assumptions. The weaknesses and indeterminations of the kinematic models are partly resolved and discussed by constructing dynamic spontaneous rupture models in Section 4. Finally in Section 5, we discuss how this rupture process takes place in the geodynamic context.

## 2 MAIN SOURCE FEATURES

Teleseismic body waves and geodetic data (GPS and INSAR) have been already used by Peyrat *et al.* (2006) to invert mechanism, fault plane, rupture extension and source time function of the main shock. Azimuthal distribution of stations well constrained the dip of the vertical nodal plane but left a large trade-off between strike and rake due to uncertainties on the low dip nodal plane. A sensitivity analysis was performed to find a set of admissible models and the preferred solution was chosen to be compatible with geodetic analysis. Then the preferred fault mechanism (strike/dip/rake) is (189°/24°/–74°), very close to the global centroid moment tensor (GCMT) mechanism (182°/23°/–81°). The source time function was a simple pulse of 17 s of duration, with a large moment release during the first 12 s. Delouis & Legrand (2007) also used teleseismic body waves and strong motion data in a combined inversion. Both studies showed average compatible results, except the following differences: they found, respectively, slow and fast rupture velocity (76 per cent and 98 per cent of S-wave speed at the focal depth); consequently, the slipping domain is larger and the source time function exhibits two pulses in Delouis & Legrand (2007) (a very narrow pulse followed by a large one); finally, both studies used slightly different mechanisms (respectively a unique 24W dipping plane versus two 15W/35W dipping planes). Therefore, it can be said that teleseismic body waves, geodetic data and long period surface wave (GCMT) have already determined, independently, the mechanism, the source duration and the seismic moment with an overall good agreement. What remains unclear for this event is the dimensions of the rupture and consequently the rupture velocity. Furthermore, the level of complexity of the rupture is still unknown. In the study of Delouis & Legrand (2007), a slip distribution has been inverted but part of the complexity of their model is introduced by need of

the method, as in all methods that aim to invert general maps of slip distribution. This study tends to constrain these elements, with the near field strong motion data only, since the other data have not been able to show a clear evidence of source complexity.

The main shock was recorded by a permanent strong motion network operated by the University of Chile and we then investigate source kinematic from these near strong motion records (see Fig. 1). Due to the depth of the source (around 100 km) none of the stations can be considered very close to the source. Only two of them, PICA and IQUI are near source and exhibit near field signal.

To infer main source features, we analyse different representations of *P* and *S* waveforms at the most shaken stations PICA, IQUI and ARIE. In particular, we analyse amplitude of acceleration (up to 3 Hz), and amplitude and polarization of velocity for different frequency bands (up to 1 Hz). Through these various treatments, we find a robust and persistent feature on all the stations: a significant minimum of seismic radiation at around 5.6 s after the *P* (or *S*) arrival. To emphasise this behaviour, we plot a stack of normalized *P* and *S* velocity waveforms on Fig. 2 with stations PICA, IQUI and ARIE, in the frequency band 0.01–0.25 Hz. On the *S* wave, we note that the minimum apparent source duration is about 11 s (the signal after 11 s can be attributed to source, structure and free surface effects). Finally, we remark that both episodes of strong radiation have different polarization. A possible hypothesis is that a second source (or a distinct asperity) has started to radiate at a different place at 5.6 s, while the first source may have already released most of its energy.

The particularities of the waveforms at PICA, IQUI and ARIE suggest that the source has radiated from at least two main asperities, the first one triggering the other one. Previous teleseismic approach by Peyrat *et al.* (2006) has not been able to separate these two asperities, probably because their signals interfere. With the combined teleseismic/strong motion inversion, Delouis & Legrand (2007) found a source time function with two episodes: a short one (from 0 to 3 s) and a long one (from 3 to 12 s). However the first

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Main features of strong motion records. Dotted and plain curves represent $L_2$ stacks of the normalized *P* and *S* amplitude waveforms at stations PICA, IQUI and ARIE. Before stacking the stations with $L_2$ norm, each station is normalized by its own $L_2$ norm computed on the observation duration (12 s for *P* wave and 15 s for *S* wave). The resulting curves are dimensionless velocity waveforms. The data are filtered in the frequency band 0.01–0.25 Hz. The signal of the *P* wave is not represented after 12 s because of the *S*-wave arrival. On both *P* and *S* waves, a clear minimum of amplitude at 5.6 s separates two episodes of strong radiation. At 11 s, most of the main rupture process is done (on the *S* wave, the slower decreasing signal after 11 s can be partly attributed to structure and free surface effects).</figcaption></figure>
<table><caption>**Table 1.** 1-D velocity model for the Tarapacá region.</caption><thead><tr><th><i>P</i> velocity (km s<sup>-1</sup>)</th><th><i>S</i> velocity (km s<sup>-1</sup>)</th><th>Density (10<sup>3</sup> kg m<sup>-3</sup>)</th><th>Depth (km)</th></tr></thead><tbody><tr><td>4.8</td><td>2.7</td><td>2.5</td><td>0</td></tr><tr><td>6.2</td><td>3.5</td><td>2.7</td><td>10</td></tr><tr><td>6.7</td><td>3.8</td><td>2.8</td><td>20</td></tr><tr><td>7.3</td><td>4.2</td><td>3.1</td><td>50</td></tr><tr><td>7.7</td><td>4.4</td><td>3.3</td><td>60</td></tr><tr><td>8.1</td><td>4.6</td><td>3.4</td><td>80</td></tr></tbody></table>

episode releases few moment and we prefer to associate it to the very beginning of the rupture process, the nucleation, rather than to a significant source expansion, that we seem to identify on PICA, IQUI and ARIE.

These elements motivate the development of the following alternative approach in order to study the source process with only the strong motion data set: we invert the source process as a set of slip patches rupturing from the hypocentre at constant velocity, with or without a delay of triggering between the patches. This approximation needs very few parameters and it is clearly justified here because the number of waveforms and their azimuthal distribution are far from optimal (see Fig. 1).

# 3 SLIP PATCH BASED KINEMATIC INVERSION WITH STRONG MOTION DATA

## 3.1 Kinematic inversion methodology

For the 2005 Tarapacá event, up to the frequency of 0.5 Hz, the velocity waveforms at stations PICA, IQUI and ARIE are relatively clear and simple due to the depth of the source (around 100 km): the angles of incidence are small so that the waveforms are relatively unaffected by the lateral heterogeneities. Therefore, concerning wave propagation, we use a simple layered model (Table 1) and we compute the Green's functions with the spectral discrete wavenumber method of Bouchon (1981), implemented in the 'axitra' code of Coutant (1990). Note that the gain due to this simple wave propagation context is partly cancelled by the poor intrinsic resolution of the accelerometric array. We choose to perform a low frequency inversion, up to 0.25 Hz, which is enough to image the main complexities of the source that radiates during 12 s. However this low frequency approach, with few stations, cannot uniquely determine a complex source model and for this reason, we parametrize the source process as a succession of slipping patches rupturing from the hypocentre. This parametrization was introduced by Vallée & Bouchon (2004) to model far field waveforms for teleseismic events. Di Carli et al. (2009) recently introduce this method to invert kinematic and spontaneous rupture models with strong motion data.

Each patch is representative of a crack of elliptical contour, with a centred elliptical slip distribution. A slip patch on the fault is determined by six parameters: two radii, one orientation, two centre coordinates and one slip amplitude. The mechanism remains constant and the rupture develops circularly from the hypocentre, at constant rupture velocity. Slip patches can intersect. A time offset can be added for each patch to simulate a delay of triggering. This parametrization allows to construct simple and complex rupture scenarios. The last parameter to discuss is the local rupture duration: if it is large (i.e. a non-negligible fraction of the time needed for the rupture to break one patch) it will have a significant role in the synthetic waveforms. Usually, mechanical models tend to show

that the source time functions are very complex and variable. To avoid this complexity, and possible trade-off with other parameters, we decide to define the local source time function for the sliding velocity as a triangular function with only 1 s of duration. The highest frequencies generated by this short source time function are not considered due to the low frequency filtering. Technically, to produce seismograms, the spatio-temporal slip history on the whole fault is discretized into a smooth source point distribution, following the kinematics explained just above.

Since the number of parameters to invert is small, we use a systematic, fully non-linear inversion method based on the neighbourhood algorithm (NA) developed by Sambridge (1999a). This technique is now widely used in source studies (François-Holden et al. 2007; Monelli & Mai 2008; Di Carli et al. 2009). Denoting a model *m*, the corresponding synthetic data *d*(*m*) and *o* the observations, the cost functions is defined as $|d(m) - o|_{L_2}$, that is, the $L_2$ norm of the differences between synthetics and data velocities, at stations IQUI, PICA, ARIE, POCO, TCP and CALA. We do not normalize the amplitudes between the stations to emphasis the role of the near-field observations, which are almost not subject to wave propagation complexity and noise. The non-linear inversion algorithm is a direct search method. NA algorithm uses simple geometrical concepts to go through the parameter space, with an optimized and variable refinement. At each iteration, models are calculated randomly. From one iteration to another, the new sample is recalculated in only some area of the parameter space corresponding to the Voronoi cells (nearest neighbour regions) of the best previous models (models with lowest misfits). NA algorithm differs from other methods such as genetic algorithm or simulated annealing in the fact that it requires only two control parameters: the number of sample $n_s$ at each iteration and the number of cells $n_r$ in which new samples are searched.

## 3.2 Preliminary strong motion data treatment

The use of strong motion displacement waveforms to build the time misfit function is not always easy for two reasons. First, the displacement misfit function depends always of the length of the time window, because of the static final displacement. Second, the static and a wide range of low frequencies can be contaminated by several factors: instrument, acquisition, installation (quality of coupling to the ground) and rotational motion (which induces frame- relative acceleration and tilting of the ground). There is no universal correction scheme (Iwan et al. 1985; Boore 2001; Boore & Bommer 2005). In some cases, acceleration data contain step-like or transient offsets. The adjustment (the baseline model) may take the form of a polynomial (Graizer 1979) or multiple linear segments (Iwan et al. 1985), but this processing is very non-unique.

In the present strong motion data, we cannot recover any stable final displacement and consequently we build the misfit function with velocities; this provides flexibility to define the length of the time window, but static information vanishes. The velocity waveforms still exhibit strong low frequency bias, so that we perform a baseline correction. In Fig. 3 (top panel), we show the result of our treatment for the north component of PICA. The baseline model on the acceleration is continuous and piecewise linear: an initial constant level, a transient linear segment during the strong motion and a final constant level. The constant levels are easily identified before and after the signal. The transient regime is adjusted by minimizing the residual average velocity at the end of the signal. On Fig. 3 (bottom panel), we plot the corrected and uncorrected
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Top panel: piecewise linear baseline corrections performed on the PICA north acceleration component. Bottom panel: 0.01–0.25 Hz filtered PICA north velocity component obtained with and without the application of the baseline correction. The filter is designed as a 0.01 Hz fourth-order highpass Butterworth combined to a 0.25 Hz second-order twice-passed lowpass Butterworth. The difference between the corrected and uncorrected signals is large in this frequency range.</figcaption></figure>

velocity waveforms, in the frequency band 0.01–0.25 Hz [highpass fourth-order and lowpass (twice-passed) second-order Butterworth filters]. The uncorrected velocity waveform exhibits a spurious low frequency content. Instead, the corrected waveform is well flattened and most of the unphysical low frequency content is removed. With this manipulation, we recover part of the true low frequency content that is traditionally removed by filtering.

### 3.3 Results of the kinematic inversions

We conduct several sets of inversions, with one, two and three slip patches and various ranges of parameter search. We present only the main results, that we obtain with one and two patches. On each patch, slip is distributed with an ellipsoidal shape, as in the crack model.

We start with the single patch model. We invert 7 parameters: the radii $a_1$, $b_1$, the orientation $\alpha_1$ taken from strike, the position, the value of the maximal slip $slip_1$ and the rupture velocity $v_r$. Instead of arbitrary centre coordinates, we define the position of the patch under the constraint that it remains in the ellipse. This let us introduce two parameters $0 < h_r < 1$ (dimensionless) and one angle $\alpha_h$ (taken from strike) to locate the ellipse relatively to the hypocentre. With $h_r = 0$ and 1, hypocentre is, respectively, at the centre and at the border of the ellipse. The ranges of search are given on Table 2. We find one optimal model denoted $k_a$. On Fig. 4, we show the slip distribution of this best model: large slip is concentrated in the hypocentral region and is mostly confined around the central part of the fault. On Table 2, we collect the misfit and the main characteristics of the inverted models. In model $k_a$, inverted rupture velocity is $3.5 \text{ km s}^{-1}$ and it can be considered as slow. Indeed, the rupture is mainly propagating along the strike direction, which corresponds to mode III, therefore assuming elastodynamic rupture theory with homogeneous stress and frictional properties on the fault, the rupture velocity can theoretically raise the S-wave speed, that is $4.6 \text{ km s}^{-1}$ at this depth. On Fig. 5, we see that model $k_a$ does not fit very satisfactory the data, especially the second S-wave arrival at station PICA at around 40 s, which corresponds to the second rupture episode that we highlighted in the introduction. In addition, the slip patch of model $k_a$ does not match the most southern part of the aftershocks, plotted on Fig. 4.

Now we examine the double patch models. In addition to the parameters of the first patch (linked to the hypocentre), the second patch is defined with seven parameters: the radii $a_2$, $b_2$, the orientation $\alpha_2$, the centre coordinates $x_e$ (towards strike) and $y_e$ (towards

<table><caption>**Table 2.** Top table: intervals of parameter search for each inversion. Bottom table: misfits and main rupture characteristics of the four inverted kinematic models $k_a$, $k_b$, $k_c$ and $k_d$.</caption><thead><tr><th></th><th>a<sub>1</sub><br>(km)</th><th>b<sub>1</sub><br>(km)</th><th>α<sub>1</sub><br>(degree)</th><th>r<sub>h</sub><br>none</th><th>α<sub>h</sub><br>none</th><th>slip<sub>1</sub><br>(m)</th><th>v<sub>r</sub><br>(km s<sup>-1</sup>)</th><th>a<sub>2</sub><br>(km)</th><th>b<sub>2</sub><br>(km)</th><th>α<sub>2</sub><br>(degree)</th><th>x<sub>e</sub><br>(km)</th><th>y<sub>e</sub><br>(km)</th><th>slip<sub>2</sub><br>(m)</th><th>t<sub>e</sub><br>(s)</th></tr></thead><tbody><tr><td>k<sub>a</sub> max</td><td>40</td><td>30</td><td>360</td><td>1</td><td>360</td><td>12</td><td>5</td><td>none</td><td>none</td><td>none</td><td>none</td><td>none</td><td>none</td><td>none</td></tr><tr><td>k<sub>a</sub> best</td><td><strong>33.68</strong></td><td><strong>16.16</strong></td><td><strong>198.58</strong></td><td><strong>0.26</strong></td><td><strong>356.4</strong></td><td><strong>10.02</strong></td><td><strong>3.35</strong></td><td>none</td><td>none</td><td>none</td><td>none</td><td>none</td><td>none</td><td>none</td></tr><tr><td>k<sub>a</sub> min</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>2.5</td><td>none</td><td>none</td><td>none</td><td>none</td><td>none</td><td>none</td><td>none</td></tr><tr><td>k<sub>b</sub> max</td><td>40</td><td>30</td><td>360</td><td>1</td><td>360</td><td>12</td><td>5</td><td>40</td><td>30</td><td>360</td><td>80</td><td>60</td><td>12</td><td>2</td></tr><tr><td>k<sub>b</sub> best</td><td><strong>30.10</strong></td><td><strong>12.59</strong></td><td><strong>194.39</strong></td><td><strong>0.26</strong></td><td><strong>1.7</strong></td><td><strong>4.62</strong></td><td><strong>3.47</strong></td><td><strong>34.98</strong></td><td><strong>5.38</strong></td><td><strong>323.9</strong></td><td><strong>38.64</strong></td><td><strong>45.03</strong></td><td><strong>4.15</strong></td><td><strong>-1.91</strong></td></tr><tr><td>k<sub>b</sub> min</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>2.5</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>-2</td></tr><tr><td>k<sub>c</sub> max</td><td>40</td><td>20</td><td>360</td><td>1</td><td>360</td><td>12</td><td>5</td><td>20</td><td>10</td><td>360</td><td>80</td><td>60</td><td>12</td><td>none</td></tr><tr><td>k<sub>c</sub> best</td><td><strong>25.29</strong></td><td><strong>13.41</strong></td><td><strong>12.58</strong></td><td><strong>0.23</strong></td><td><strong>143.3</strong></td><td><strong>4.98</strong></td><td><strong>3.48</strong></td><td><strong>18.72</strong></td><td><strong>7.21</strong></td><td><strong>178.2</strong></td><td><strong>53.50</strong></td><td><strong>40.65</strong></td><td><strong>5.32</strong></td><td>none</td></tr><tr><td>k<sub>c</sub> min</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>2.5</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>none</td></tr><tr><td>k<sub>d</sub> max</td><td>40</td><td>30</td><td>360</td><td>1</td><td>360</td><td>12</td><td>5</td><td>40</td><td>30</td><td>360</td><td>80</td><td>60</td><td>12</td><td>none</td></tr><tr><td>k<sub>d</sub> best</td><td><strong>35.26</strong></td><td><strong>19.98</strong></td><td><strong>10.57</strong></td><td><strong>0.23</strong></td><td><strong>179.1</strong></td><td><strong>2.77</strong></td><td><strong>4.22</strong></td><td><strong>35.58</strong></td><td><strong>7.41</strong></td><td><strong>192.15</strong></td><td><strong>58.26</strong></td><td><strong>41.84</strong></td><td><strong>2.83</strong></td><td>none</td></tr><tr><td>k<sub>d</sub> min</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>2.5</td><td>0</td><td>0</td><td>0</td><td>40</td><td>20</td><td>0</td><td>none</td></tr></tbody><thead><tr><th></th><th>Model</th><th>Misfit</th><th>M<sub>o</sub> (Nm)</th><th>Max. slip (m)</th><th>v<sub>r</sub> (km s<sup>-1</sup>)</th><th>Delay t<sub>e</sub> (s)</th><th>Rise time (s)</th></tr></thead><tbody><tr><td>k<sub>a</sub></td><td>24.5 per cent</td><td>2.64 × 10<sup>20</sup></td><td>9.92</td><td>3.35</td><td>none</td><td>1 (fixed)</td><td></td></tr><tr><td>k<sub>b</sub></td><td>13.2 per cent</td><td>2.83 × 10<sup>20</sup></td><td>4.64</td><td>3.47</td><td>-1.91</td><td>1 (fixed)</td><td></td></tr><tr><td>k<sub>c</sub></td><td>13.8 per cent</td><td>2.72 × 10<sup>20</sup></td><td>5.26</td><td>3.48</td><td>0 (fixed)</td><td>1 (fixed)</td><td></td></tr><tr><td>k<sub>d</sub></td><td>13.7 per cent</td><td>2.97 × 10<sup>20</sup></td><td>2.82</td><td>4.22</td><td>0 (fixed)</td><td>1 (fixed)</td><td></td></tr></tbody></table>

**Notes:** The values of the best model are in bold. Model $k_a$ was obtained by searching one single patch and $k_b$, $k_c$ and $k_d$ by searching two patches. In $k_b$, a delay of triggering between the two patches has been searched. In $k_c$ and $k_d$, this delay is fixed to zero. For $k_c$, the range for the patch dimensions has been limited to smaller values to force the finding of a rupture model with lower rupture velocity. Model $k_d$ is the best model for the large range of parameter search and without the delay of triggering.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Kinematic inverted models. Slip colour scale ranges from blue (0 m) to red (5 m). The star is the epicentre and red dots are the aftershocks. Model $k_a$ is the single patch model. Model $k_b$, $k_c$ and $k_d$ are three inverted double patch models, with one main patch attached to the hypocentre and one secondary patch. These three models are obtained with different additional assumptions. For $k_b$: an additional delay of triggering for the secondary patch, for $k_c$: a limited range of exploration for the ellipses dimensions to force the finding of a slower rupture velocity and for $k_d$: a wide range of exploration for the ellipses dimensions. The main characteristics of these models are reported on Table 2</figcaption></figure>

dip) taken from the northeast corner of the prescribed 80 × 60 km² fault domain, the maximal slip *slip*₂ and a delay of triggering *t*ₑ. By inverting the 14 parameters, we find one model $k_b$ with a better fit than model $k_a$ (see Table 2 for details). However, in model $k_b$, negative delay, stretched shape, far position and oblique orientation of the second patch suggest that the resulting kinematics do not correspond to a reasonable rupture process but rather to an artefact.

To propose another solution, we investigate models with $t_e = 0$ (no delay). With this reduced 13 parameters set, we find model $k_d$, with a second patch located southwest to the hypocentre (Fig. 4). Compared to model $k_b$, model $k_d$ has a more reasonable slip pattern not only because the second patch is closer to aftershocks but also because the first patch covers entirely the aftershock area. Furthermore, this model $k_d$ has a larger (better) seismic moment, a rupture velocity (4.22 km s⁻¹) just below the S-wave speed for an almost equal misfit of 13.7 per cent (see Table 2). Finally, we test the effect of a reduction of the parameter ranges: by reducing the maximal dimensions of the ellipses (see Table 2), we find model $k_c$. Model $k_c$ has a similar slip pattern than $k_d$ (see Fig. 4) but it is rescaled in space due to a lower rupture velocity (3.48 km s⁻¹).

To conclude here, the single patch model $k_a$ is not satisfactory. The qualitatively different double patch models $k_b$, $k_c$ and $k_d$ found with different assumptions produce almost the same waveforms (Fig. 5) and comparable fits (Table 2). However, by introducing the triggering delay, we have probably produced an artefact in model $k_b$. Instead, model $k_d$ (without the delay of triggering) seems a reasonable solution. By truncating the intervals of search of the ellipse

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Data and synthetic waveforms (velocities) for single patch model $k_a$ and for double patch models $k_b$, $k_c$ and $k_d$. All waveforms are filtered in the frequency band 0.01–0.25 Hz. Single patch model $k_a$ does not reproduce waveforms complexities on IQUI and PICA whereas double patch models fit better the data, and almost equivalently.</figcaption></figure>

dimensions, we find another model $k_c$, similar to $k_d$ with a smaller spatial extension and thus a slower rupture velocity. This shows that the spatio-temporal relocation of the slip patches is not obvious and that it is almost vain to search and invert for a third slip patch. We point out that all models have a lack of about 75 per cent of seismic moment compared to the GCMT value (5.32 × 10²⁰ N m). This systematic bias is mainly due to the short rise-time (1 s) and to the circular shape of the rupture front. These assumptions tend to produce larger velocity amplitudes on waveforms, for a given fault slip. This problem will be corrected in the spontaneous rupture models. Finally, model $k_d$, that was preferred by the inversion procedure has the better seismic moment, according to previous studies.

## 3.4 Uncertainty analysis of kinematic models

Model $k_a$ (one patch only) and model $k_c$ (reduced ranges of exploration) can be considered as intermediate models included in the set of models generated during the inversion of model $k_d$. Model $k_b$, determined with a more general inversion (with one more parameter) gives an unrealistic solution. Here, we discuss only the probability of models, with two patches, in the set of the models generated during its inversion by the NA algorithm.

The convergence of the parameters is a preliminary indicator of the solution reliability. On Fig. 6, we plot the inverted parameters versus the iterations of the NA algorithm (25 300 models were generated), with a grey colour scale indicating the level of misfit, for the model $k_d$. Models $k_b$ and $k_c$ shows similar convergence and are not shown in the figure. We observe that enough models are generated to obtain a stable unique branch of solution, with reasonable fluctuations compared to the intervals of search. In addition, parameters not concerned by periodicity (unlike angles) are not obstructed by the boundaries of search.
Models generated by neighbourhood algorithm

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Neighbourhood algorithm convergence to the kinematic model $k_d$. Each panel shows the value of an inverted parameter (vertical axis) plotted against the model number (horizontal axis). $a_1, b_1, a_2$ and $b_2$ are the lengths of the axes of the elliptic slip patches (km). $\alpha_1$ and $\alpha_2$ are their orientation (degree). $h_r$ and $\alpha_h$ control the position of the first ellipse relatively to the hypocentre. $slip_1$ and $slip_2$ are the maximal slips on each ellipse (m). $v_r$ is the rupture velocity (km s$^{-1}$). $x_c$ and $y_c$ are the position of the centre of the second ellipse (km).</figcaption></figure>

A supplementary uncertainty analysis is performed on the inversion using an appraising procedure: the Bayesian neighbourhood algorithm (NA-bayes or NAB) of Sambridge (1999b). This technique has already been tested by Monelli & Mai (2008) for kinematic rupture inversion with strong motion data. NAB algorithm facilitates error analysis on the set of models generated during the NA search stage (inversion). NAB makes quantitative inferences from the entire set produced by the direct search method, and allows measures of resolution and trade-offs, within a Bayesian framework. From the Bayesian viewpoint, the solution to the inverse problem is the marginal (posterior) probability density function (PDF) computed for each parameter. Using the whole set of models founded during the inversion, NAB computes a geometrical approximation of the marginal PDF by simply setting the known probability of each model to be constant inside Voronoi cell. New samples generated according to the approximate probability are produced using a Gibbs sampler. A Gibbs sampler generates samples by performing a random walk in the model space (see Sambridge 1999b). The main advantage of NAB is to approximate the marginal PDF without the need of running the forward simulation for each resampled models.

Following Sambridge (1999b), the marginal PDF can be written as

$$P(m) = A \exp \left\{ -\frac{1}{2} \left[ \frac{|d(m) - o|_{L_2}}{\epsilon|o|_{L_2}} \right]^2 \right\},$$

where A is a normalization factor and $\epsilon$ a parameter of regularization. Parameter $\epsilon$ is formally linked to the random (and supposed Gaussian) level of noise in the observed data. This parameter is fundamental for the appraisal stage of the NAB. As we mention in Section 3.2, our strong motion data exhibit several sources of uncertainty, including random and not random perturbations like the piecewise linear baseline seen on the accelerations. We have partly corrected this systematic bias. However, to stay conservative, we choose $\epsilon = 0.1$. This corresponds to a confidence of 90 per cent in the observed data.

Following this methodology we compute for each model parameter its marginal PDF. Applied to model $k_b, k_c$ and $k_d$, the NAB algorithm generates 100 000 resampled models. The resulting marginal PDF are shown on Fig. 7 for each inverted parameter and for several additional deduced parameters: the surfaces of the two patches ($s_1$ and $s_2$), the total seismic moment ($M_0$), the approximated rupture time of the first patch ($t_{r1}$) and the coordinates of the centroid on the

fault ($x_c$ towards strike and $y_c$ towards dip). Vertical bars indicate the position of the corresponding final models $k_b, k_c, k_d$ (with the best misfit).

The marginals provide a useful way of testing how the sampled algorithm converges. The final models (vertical bars in Fig. 7) correspond for most of the parameters to the marginal PDF peaks, indicating an acceptable convergence of the previous NA algorithm. The width of the marginals is an indicator of the constraint that may be placed on each parameter: the most the distribution is narrow; the most the data provided information on parameter. Note that the marginal PDF for each parameter do not follow the same distribution for each model, and they do not show always a ‘Gaussian’ shape. The only marginals with ‘Gaussian’ shape correspond to slip, rupture velocities, moment, rupture time, centroid position and patches positions. For the other parameters, the skewed distributions are probably due to the limitations of the parameter space; another explanation is that some parameters are correlated, like the sizes (a and b) of the patches. Indeed, the surfaces of the patches are better constrained, except for model $k_b$. But the marginal distributions show that all these models share some large-scale features. By observing the spreading of the distributions and the position of the models, we make the following remarks. As expected, seismic moment and centroid are well constrained and consistent. Hopefully, surface ($s_1$), position ($r_h, \alpha_h$), slip ($slip_1$), approximate rupture time of the first slip patch ($t_{r1} = \sqrt{a_1b_1}/v_r$), are constrained and consistent, although a slip-surface trade-off is visible. We make the same remarks for the second patch.

## 4 DYNAMIC MODELLING

The preferred kinematic model $k_d$ shows that the source process is likely happening on two main asperities. The uncertainty analysis shows that this model is relatively well constrained, however we can force the kinematic inversion to find reasonable alternative models with slower rupture velocity, as in model $k_c$. Kinematic study (Section 3) is based on arbitrary assumptions, like the fixed short rise time (1 s). Hence one reasonable alternative is to build spontaneous rupture models, by trial and error, to discuss our results. To this end, we consider the results of five mechanical models $d_0, d_1, d_2, d_3$ and $d_4$ based on classical elastodynamic rupture theory and inspired by kinematic models $k_c$ and $k_d$. These models are constrained qualitatively by the strong motion data. In addition,
Marginal probability density functions

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Results of the resampling algorithm for the three kinematic models with two patches. For each model $k_b$, $k_c$ and $k_d$, each panel represents the marginal PDF for the inverted parameters and for some additional deduced parameters. The first row represents the PDF of the inverted parameters of the first slip patch, and the rupture velocity. The second row represents the PDF of the inverted parameters of the second slip patch. The third row represents the PDF of the following additional parameters computed during the resampling: the surfaces of the ellipses ($s_1$ and $s_2$), the total seismic moment, the approximated rupture time on the first ellipse ($t_{r1}$) and the coordinates of the rupture centroid on the fault ($x_c$ and $y_c$). Each x-axis is over the complete range of the parameters space. The vertical lines show the results of each model. Each curve is scaled to the same maximum height, not the same area.</figcaption></figure>
<table><caption>**Table 3.** Seismic moments of the dynamic models compared to previous studies.</caption><thead><tr><th>Moment M<sub>o</sub><br>unit</th><th>CMT<br>NEIC/USGS</th><th>CMT<br>GCMT</th><th>Teleseismic<br>Peyrat et al. (2006)</th><th>Geodetic<br>Peyrat et al. (2006)</th><th>Joined inversion<br>Delouis & Legrand (2007)</th><th>d<sub>0</sub></th><th>d<sub>1</sub></th><th>d<sub>2</sub></th><th>d<sub>3</sub></th><th>d<sub>4</sub></th></tr></thead><tbody><tr><td>10<sup>20</sup> N m</td><td>6.5</td><td>5.32</td><td>5.4</td><td>5.8</td><td>5.47</td><td>5.00</td><td>5.63</td><td>5.61</td><td>5.63</td><td>5.56</td></tr></tbody></table>

*Notes:* Dynamic models *d*<sub>1</sub>, *d*<sub>2</sub>, *d*<sub>3</sub> and *d*<sub>4</sub> are constrained to fit the reference seismic moment *M*<sub>o</sub> = 5.6 × 10<sup>20</sup> Nm. NEIC/USGS is the National Earthquake Information Center solution (http://earthquake.usgs.gov/regional/neic/), GCMT is the Global CMT solution (http://www.globalcmt.org/).

models *d*<sub>1</sub>, *d*<sub>2</sub>, *d*<sub>3</sub> and *d*<sub>4</sub> are constrained by a reference value of seismic moment *M*<sub>o</sub> ≈ 5.6 × 10<sup>20</sup> N m, a compromise consistent with earlier works (see Table 3).

The dynamic rupture is simulated using a 3-D finite-difference method and strong-motion data are computed at stations with the spectral method (code axitra). At each point of the fault, we use a rupture criterion in which the yielding resistance τ<sub>y</sub> is defined as a first order ordinary differential equation in time

$$ \dot{\tau}_y(t) = \frac{1}{D_c} \{ V_c[\tau_s - \tau_y(t)] - \delta\dot{u}(t)[\tau_y(t) - \tau_d] \}, \quad \tau_y(0) = \tau_s, $$

where δu is the fault slip, τ<sub>s</sub> is the (static) maximal resistance and τ<sub>d</sub> is the (dynamic) residual resistance. D<sub>c</sub> is a characteristic slip weakening distance, and V<sub>c</sub> is a characteristic rate weakening.

This phenomenological friction law is a dimensional version of the one defined by Nielsen et al. (2000), without considering the viscosity part of their governing equation. This law also described in Nielsen & Carlson (2000) is similar to a rate and state friction law, and is based on the constitutive relation for friction between lubricated surfaces proposed by Carlson & Batista (1996). From many frictional forms that could be considered for our kind of problem, we choose this formulation because of the simplicity of the slip-weakening onset, the time-dependent healing and the steady-state rate weakening, the phenomena needed to describe a wide range of dynamical features. The main interest of this friction law is to deal with a smooth transition from cracklike to self-healing pulselike ruptures.

Rupture contains several regimes dominated by time dependent restoration at rest, slip weakening at rupture onset and rate weakening at steady state. Parameter D<sub>c</sub> and V<sub>c</sub> control respectively the sharpness of slip and rate weakening behaviours. In elastodynamic rupture experiments, D<sub>c</sub> regularizes and determines the rupture front thickness. A large enough V<sub>c</sub> (i.e. comparable to the sliding velocity behind the rupture front) selects the pulse rupture mode instead of the classical crack rupture mode. The last important parameter is the pre-stress τ<sub>0</sub>, which leads to the definition of the dynamic stress drop τ<sub>0</sub> − τ<sub>d</sub> (which gives energy to rupture) and to the strength excess τ<sub>s</sub> − τ<sub>0</sub> (which tends to stop it). Starting from the kinematic model k<sub>d</sub>, we build the dynamic models by defining four zones: the nucleation zone, the first and second elliptical asperities and the interstitial zone (leading between and around the asperities). The average values of pre-stress, strength excess, D<sub>c</sub> and V<sub>c</sub> are reported on Table 4. The detailed maps of the parameters are given only for model d<sub>1</sub> on Fig. 8.

To begin, we build two spontaneous rupture models d<sub>0</sub> and d<sub>1</sub> with purely slip dependent friction (V<sub>c</sub> = 0), that aim to mimic kinematic model k<sub>d</sub>, having a high rupture velocity. Model d<sub>0</sub> has an unbreakable interstitial zone (no slip) and it is constrained by strong motion only. Model d<sub>1</sub> has a breakable interstitial zone and it is constrained by strong motion and seismic moment. On Fig. 9, we plot the snapshots of their sliding velocity and elastic stress every second. Kinematic model k<sub>d</sub> is also represented for comparison. Both models d<sub>0</sub> and d<sub>1</sub> reproduce the features of k<sub>d</sub>, with a high rupture velocity. The average rupture velocity measured in the anti-plane

<table><caption>**Table 4.** Mechanical parameters of the five dynamic models d<sub>0</sub>, d<sub>1</sub>, d<sub>2</sub>, d<sub>3</sub> and d<sub>4</sub>.</caption><thead><tr><th>Model</th><th>Average<br>dynamic stress<br>drop (MPa)</th><th>Average<br>strength<br>excess (MPa)</th><th>Average<br>D<sub>c</sub><br>(m)</th><th>Average<br>V<sub>c</sub><br>(m s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>d0 / nucleation</td><td>15.0</td><td>0</td><td>0.5</td><td>0</td></tr><tr><td>d0 / ellipse 1</td><td>6.0</td><td>1.0</td><td>0.5</td><td>0</td></tr><tr><td>d0 / ellipse 2</td><td>11.5</td><td>1.05</td><td>0.5</td><td>0</td></tr><tr><td>d0 / interstitial</td><td>—</td><td>+∞</td><td>—</td><td>—</td></tr><tr><td>d1 / nucleation</td><td>15.0</td><td>0</td><td>0.5</td><td>0</td></tr><tr><td>d1 / ellipse 1</td><td>6.5</td><td>1.0</td><td>0.5</td><td>0</td></tr><tr><td>d1 / ellipse 2</td><td>11.5</td><td>1.5</td><td>0.5</td><td>0</td></tr><tr><td>d1 / interstitial</td><td>-3.5</td><td>0</td><td>0.5</td><td>0</td></tr><tr><td>d2 / nucleation</td><td>20</td><td>0</td><td>0.5</td><td>0</td></tr><tr><td>d2 / ellipse 1</td><td>8.5</td><td>6.5</td><td>0.66</td><td>0</td></tr><tr><td>d2 / ellipse 2</td><td>17</td><td>3.0</td><td>0.5</td><td>0</td></tr><tr><td>d2 / interstitial</td><td>-3.15</td><td>0</td><td>0.5</td><td>0</td></tr><tr><td>d3 / nucleation</td><td>20</td><td>0</td><td>0.5</td><td>0.3</td></tr><tr><td>d3 / ellipse 1</td><td>16</td><td>9</td><td>0.66</td><td>0.3</td></tr><tr><td>d3 / ellipse 2</td><td>17</td><td>6.5</td><td>0.5</td><td>0</td></tr><tr><td>d3 / interstitial</td><td>0</td><td>0</td><td>0.5</td><td>0</td></tr><tr><td>d4 / nucleation</td><td>20</td><td>0</td><td>0.5</td><td>0.3</td></tr><tr><td>d4 / ellipse 1</td><td>16</td><td>9</td><td>0.66</td><td>0.3</td></tr><tr><td>d4 / ellipse 2</td><td>17</td><td>3.5</td><td>0.5</td><td>0</td></tr><tr><td>d4 / interstitial</td><td>-2.65</td><td>0</td><td>0.5</td><td>0</td></tr></tbody></table>

*Notes:* The parameters are averaged on each zone of the model, that is, the nucleation zone at the hypocentre, the two elliptical asperities and the interstitial (and surrounding) zone lying around and between the two asperities. Models d<sub>0</sub> and d<sub>1</sub> have crack-like mode and high rupture velocity (close to the 4.6 km s<sup>-1</sup> S-wave velocity). d<sub>0</sub> has no interstitial rupture but lower seismic moment. Models d<sub>2</sub>, d<sub>3</sub> and d<sub>4</sub> differ from models d<sub>0</sub> and d<sub>1</sub> mainly on the large asperity (ellipse 1), where they have low rupture velocity (down to 3.5 km s<sup>-1</sup>), due to the combination of significant strength excess (6.5 MPa) and increasing D<sub>c</sub> from the hypocentre. On the large asperity, model d<sub>2</sub> is crack-like and model d<sub>3</sub> and d<sub>4</sub> are pulse-like, due to the introduction of rate dependence in the friction (V<sub>c</sub> = 0.3 m s<sup>-1</sup>). Model d<sub>2</sub>, d<sub>3</sub> and d<sub>4</sub> have smaller dimension and larger stress drops than d<sub>0</sub> and d<sub>1</sub>. Models d<sub>3</sub> and d<sub>4</sub> have even larger stress drop (16 MPa in average and locally 32 MPa) to sustain the pulse-like rupture mode. Model d<sub>3</sub> and d<sub>4</sub> have, respectively, small and large interstitial rupture zone. Model d<sub>4</sub> exhibits a very large interstitial (then very surrounding) area of low radiative slip. All interstitial zones have negative stress drop.

direction is close to the S-wave speed (4.6 km s<sup>-1</sup>). The waveforms computed at stations PICA, IQUI and ARIE reproduce qualitatively well the data (see Fig. 10). During the search, we notice that the position of the second asperity and its time of triggering are highly constrained by the shape of the large positive S wave pulse on the east component at IQUI and by the timing of the second S wave arrival at PICA. We remark that slip is larger in dynamic models d<sub>0</sub> and d<sub>1</sub> than in kinematic model k<sub>d</sub>. Indeed the dynamic rupture front does not radiate as circularly (efficiently) as in the kinematic model k<sub>d</sub> and, furthermore, the dynamic rise time is much larger than the kinematic one (1 s). Finally, the only difference between model d<sub>0</sub> and d<sub>1</sub> is the seismic moment, still too low of 11 per cent
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Fault maps of the dynamic parameters of the dynamic model $d_1$. Left-hand panel: distribution of shear pre-stress (or dynamic stress drop) $\tau_0 - \tau_d$. Right-hand panel: distribution of shear strength excess $\tau_s - \tau_0$. Around the hypocentre (black star), particular distributions of parameters (larger stress drop and no strength excess) are needed to nucleate the rupture. Here parameter $D_c = 0.5$ m is homogeneous. The stress drop is slightly negative in the interstitial zone between the two asperities, in order to obtain a non-radiative slip process.</figcaption></figure>

in model $d_0$. This lack of seismic moment is recovered in $d_1$ by letting the rupture propagate in the interstitial zone, in order to produce a halo of low (non-radiative) slip. This adjusted non-radiative slip is obtained with null strength excess and slightly negative stress drop (see Table 4 and Fig. 8). To conclude, model $d_1$ is a convenient rupture model, minimalist, purely slip dependent, with high rupture velocity.

We now propose alternative solutions to model $d_1$, having lower rupture velocity on the first main asperity. For this, we construct three different models $d_2$, $d_3$ and $d_4$, to be compared with model $d_1$. These models have almost the same seismic moment, very close to the reference value. Their comparison is shown through rupture histories and sensitive waveform comparisons on Fig. 11. In alternative models $d_2$, $d_3$ and $d_4$ a slow average rupture velocity (down to $3.5 \text{ km s}^{-1}$) is obtained on the first asperity by using higher strength excess and by increasing linearly the critical slip distance $D_c$ from the hypocentre to stabilize as possible the rupture velocity. Model $d_2$ is purely slip dependent. On Fig. 11, we see that waveforms

<figure><img src="image_10.png" /><figcaption>**Figure 10.** Comparison of the velocity waveforms of the dynamic models $d_0$ and $d_1$ with the strong motion data. Waveforms are dominated by S wave but P wave and intermediate field are visible on PICA and IQUI. All waveforms are filtered in the frequency band 0.01–0.25 Hz. A general good agreement is observed. The complexity of the waveforms is similar. Phases of arrivals and pulses are correct. Although local differences of amplitude and pulse durations may be observed (especially on the north components), amplitudes are globally satisfactory for such models built by trial and error. Models $d_0$ and $d_1$ produce waveforms with insignificant differences.</figcaption></figure>

produced by $d_2$ have a large lack of amplitude, expressing an insufficient moment release rate, may be combined to a lack of directivity at PICA. In models $d_3$ and $d_4$, we introduce the rate effect ($V_c = 0.3 \text{ m s}^{-1}$) to produce a pulse rupture mode. Model $d_3$ produces more amplitude, however too much at IQUI and not enough on the first pulse at PICA. Furthermore, to respect the moment constraint, model $d_3$ needs a higher slip in the interstitial zone and consequently a transient supershear rupture is provoked between the asperities. Finally, we propose to spread the halo of non-radiative slip and we end up with model $d_4$. Compared to $d_1$, model $d_4$ fits the data better on IQUI (the signal decreases more rapidly); however it fits the data less well on PICA (the first S arrival has not the good shape).

Hence models $d_2$ and $d_3$ can be rejected. Models $d_1$ and $d_4$ are convenient and difficult to discriminate with data. However, model $d_4$ suggests that a slow rupture on the first asperity may be associated to a pulse mode rupture and a very large embedding zone of

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Velocity and shear stress snapshots of the kinematic model $k_d$ and of the spontaneous rupture models $d_0$ and $d_1$. Final slip and stress distributions are represented on the last column at right. First row: velocity snapshots of the kinematic model $k_d$. Second and third row respectively: velocity and shear stress snapshots of the dynamic model $d_0$. Fourth and fifth row, respectively: velocity and shear stress snapshots of the dynamic model $d_1$. Source duration of all these kinematic and dynamic models is approximately 11 s. In model $d_0$, the interstitial area between the asperities is stuck and it is subject to a very high stress concentration. In model $d_1$, a halo of non-radiative diffuse slip is created in the interstitial area, where stress concentration have been strongly attenuated.</figcaption></figure>
<figure><img src="image_12.png" /></figure>

**Figure 11.** Left-hand side: Slips, rupture time and space-time characteristics of the four dynamic models *d*₁, *d*₂, *d*₃ and *d*₄ (with equivalent seismic moment). On the slip maps (first and third columns), black contour lines represent the arrival of the rupture front every second. The vertical dotted white lines are the axes along which we map the rupture in space and time. On the second and fourth columns, we map the sliding velocity as a function of space (along the dotted white lines) and time (horizontal axis). On the large asperity, models *d*₁ and *d*₂ have crack-like rupture, whereas models *d*₃ and *d*₄ have pulse-like rupture. Model *d*₁ has a rupture velocity close to the *S*-wave speed (4.6 km s⁻¹). Model *d*₂, *d*₃ and *d*₄ have lower rupture velocity (down to 3.5 km s⁻¹). Model *d*₂ has a slower initial moment release and model *d*₃ is radiative in the interstitial zone (local supershear transition). Right-hand side: Comparison of waveforms (dynamic models *d*₁, *d*₂, *d*₃, *d*₄ and data). Only the east components at IQUI and PICA are represented since they are very sensitive to the detail of the rupture process. To enhance the differences of the models, the frequency band is enlarged to 0.01–0.33 Hz, instead of the one used for inversion (0.01–0.25 Hz). All these models have been adjusted to fit qualitatively the data at best. For model *d*₂, the first *S*-wave pulse on PICA has too large duration and too small amplitude and the signal decreases too slowly at IQUI, hence corresponding to an insufficient moment release rate. For model *d*₃, the first positive *S*-wave pulse at PICA is too large and the amplitude at IQUI is too large, due to the radiation produced in the interstitial zone. Finally, models *d*₁ (fast crack-like rupture) and *d*₄ (slow pulse-like rupture associated to extended interstitial slip) fit better the strong motion data.

non-radiative slip. In that sense, model *d*₄ is physically more complex than model *d*₁.

## 5 DISCUSSION AND CONCLUSION

A determination of both the slip distribution and the highly radiative zones on the fault is crucial to interpret the importance of the Tarapacá earthquake in the seismic cycle of the subduction. The slip distributions of kinematic model *k*<sub>d</sub> and dynamic model *d*<sub>1</sub> are in good overall agreement with the final model found by Delouis & Legrand (2007). The common points with their study are: a maximum slip located in the region where the aftershocks are distributed, similar rupture velocity and directivity and extension of the slip domain far out of the aftershocks location to the west. The advantage of our simple solutions is to identify the main features of the source complexity. Our main significant results are: the identification of two very distinct sources and the confirmation that the rupture extends down in the slab. It is interesting to note that Delouis & Legrand (2007) had to increase the dip of the rupture plane to the west to fit teleseismic data. They also concluded that strong motion data require that slip overreach the western boundary of the aftershocks cluster. Also their resulting source image is slightly different, their conclusions are very consistent with our two sources model. And finally, whatever the dynamic model is considered (*d*₁ or *d*₄) the mechanical characteristics of the event (local stress drop and rupture velocity) can be considered as normal and do not break the usual scaling laws (dominated by shallower events).

The location of the Tarapacá event in the geodynamic context is difficult. Indeed seismic imaging gives good resolution only up to 100 km depth. Deeper, only the small Wadati–Benioff seismicity

remains the main tool of observation. Last decades, improvement of both seismic imaging and earthquake locations have shown that the main zone of intermediate-depth earthquakes lies at the top of the subducting plate but not within its cold interior (Abers 1992). Hacker *et al.* (2003) proposed that earthquakes occur in subducting slab where dehydration is expected and are absent from parts of slabs predicted to be anhydrous. They proposed four petrological and seismic distinct layers in subducting oceanic plate, with one seismic gap between double seismic zones corresponding spatially to the thermally cold core of the slab. Today, multiple layering of Wadati–Benioff zones small seismicity at intermediate depths is a much more common feature than previously thought and it seems to be identified almost everywhere (Brudzinski *et al.* 2007). Intermediate-depths seismicity in northern Chile delineates a double seismic zone (DSZ) (Rietbrock & Waldhauser 2004; Comte & Suárez 1994), although detailed resolution remains limited by the uncertainties of the hypocentre depth locations. Upper plane earthquakes are inferred to occur within the subducting oceanic crust, whereas lower plane earthquakes occur in the subducting oceanic mantle (Kirby 1995). In some cases, the DSZ have specific state of stress for both seismic zones, either strongly clear compressional or tensional, nevertheless, north Chile fault plane solutions for both layers show a wide variability (Rietbrock & Waldhauser 2004; Comte & Suárez 1994). Extensional faulting seems to predominate in the upper seismic zone with downdip extension. Although the lower seismic zone is poorly resolved, (Rietbrock & Waldhauser 2004) estimates that upper and lower band of seismicity are separated by only 10 km.

So, the true location of the slab at 100 km depth is not well known. Several sources of uncertainties comes from the seismic
<figure><img src="image_13.png" /><figcaption>**Figure 12.** West–east cross-section (20°S): Situation of the slab-pull rupture, in a geodynamic context extrapolated from the ANCORP'96 study (seismic profile, and temperature profile). Green to blue filled circle represent background seismicity ($M > 4.5$) from NEIC and Engdahl *et al.* (1998). Black circled circle represent the aftershocks. The blue to red area shows the slip distribution of the spontaneous rupture model $d_1$, projected on the cross-section with a colour scale ranging from 0 to 5 m.</figcaption></figure>

imaging, from the thermal structure of the subduction zones, from the location of the seismicity and from the underlying assumptions that link slab structure and seismicity. The uncertainty on the vertical slab location may be about 10 km. However, we propose a relocation of the event in the geodynamic context. Fig. 12 shows a west–east cross-section with the slip distribution of model $d_1$ replaced in the subduction of north Chile. We extrapolate it from the seismic reflection profile of geophysical experiment ANCORP'96 (ANCORP Working Group 2003) across the Andes at 21°S. Keeping in mind the large downdip extension of the source (50 km), we see that the Tarapacá earthquake ruptured through the entire lithosphere, as also suggested by Kausel & Campos (1992) for the 1950 intermediate depth earthquake that occurred landward of 1995 Antofagasta earthquake. For the Tarapacá earthquake, the rupture seems to start at the base of the crustal slab and then to propagate both in the oceanic crust and in the mantle lithosphere. Then, we point out the possibility that the second source was triggered at the base of the mantle lithosphere, close to the isotherm 500 °C, where a few aftershocks have been located.

Ide & Takeo (1996) constructed a dynamic rupture model of the intermediate depth 1993 Kushiro-oki earthquake ($M_w = 7.6$) that occurred also on a subhorizontal plane. The mechanism was downdip extension with high stress drop, and in this case, the rupture started near the lower plane of DSZ, propagated toward the region of downdip compression and was stopped presumably by this change of the stress state, and then did not reach the upper plane of the DSZ (Kosuga *et al.* 1996; Ide & Takeo 1996; Ozel & Moriya 1999). In Japan, the stresses in DSZ are in-plate compression for the upper zone and in-plate extension for the lower one (Kosuga *et al.* 1996). This earthquake had peculiar features in that the slip on the fault plane was restricted within an extremely small area between the upper and lower planes of the DSZ in the subducting slab beneath Hokkaido, showing a sort of decoupling of the two seismic zones.

Opposite to this 1993 Kushiro-oki earthquake, our study tends to confirm that Tarapacá earthquake have ruptured both seismic zones of the DSZ. It suggests that in north Chile the presence of the two seismic zones of the DSZ do not really influence such rupture. At that point, we suggest that the Tarapacá rupture is part of the process that could lead to a long-term slab separation. It is interesting to note that the separation nucleates in the actively seismic zone, close to the base of the crustal slab, and extents to the isotherm 500 °C, which corresponds to the vanishing of the seismicity and probably the lower boundary of the lithosphere. It is also interesting to note that the breaking occurred at the top of the 100–150 km depth seismicity cluster. Concerning the lateral (north–south) extension

of the rupture, from a detailed tomographic image Rietbroek *et al.* (2006) found evidence for a trench perpendicular tear in the slab of the subducting plate at 21°S. They suggested that the Tarapacá earthquake occurred at the tip of a scissor-like structure currently breaking up the subducting Nazca plate north of 21°S.

In general, the mechanism for any seismicity below 70 km is still an ongoing debate. Warren *et al.* (2008) identified fault planes between 100 and 300 km depth in Middle America, all subhorizontal. Similar results in Tonga-Kermadec and Central South America subduction zone suggests that the earthquake-generating mechanism is controlled by pressure rather than temperature or other tectonic parameters. Jiao *et al.* (2000) suggested that the intermediate-depth earthquake occur along faults formed prior to subduction in the outer rise. This hypothesis seems very reasonable for the Tarapacá rupture since it neither showed a large sensitivity to the layered slab structure nor encountered a large resistance (standard rupture velocity), as if rupture propagated easily on an inherited, well defined, plane of weakness. However, the presence of the second deep source, separated from the main one by an eventual low non-radiative slip moderates this conclusion: a complex intraslab seismic cycle may be possible.

## ACKNOWLEDGMENTS

This work was supported by the ANR French National Research Agency (ANR-05-CATT-014). We thank editor Massimo Cocco and two anonymous reviewers for their helpful comments on the manuscript. We are grateful to Malcolm Cambridge for providing his NA and NAB algorithms. We thank J.P. Vilotte, P. Bernard, J.B. de Chabalier and M. P. Bouin for extended discussions and R. Boroschek from the University of Chile for providing the strong motion data.

## REFERENCES

* Abers, G.A., 1992. Relationship between shallow- and intermediate-depth seismicity in the eastern Aleutian subduction zone, *Geophys. Res. Lett.*, **19**, 2019–2022.

* ANCORP Working Group, 2003. Seismic imaging of a convergent continental margin and plateau in the central Andes (Andean Continental Research Project 1996 (ANCORP'96)), *J. geophys. Res.*, **108**, doi:10.1029/2002JB001771.

* Astiz, L., Lay, T. & Kanamori, H., 1988. Large intermediate-depth earthquakes and the subduction zone process, *Phys. Earth planet. Inter.*, **55**, 80–166.

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

* Boore, D.M., 2001. Effect of baseline corrections on displacements and response spectra for several recordings of the 1999 Chi-Chi, Taiwan, earthquake, *Bull. seism. Soc. Am.*, **91**(5), 1199–1211.

* Boore, D.M. & Bommer, J.J., 2005. Processing of strong-motion accelerograms: needs, options and consequences, *Soil Dyn. Earthq. Eng.*, **25**(2), 93–115.

* Brudzinski, M.R., Thurber, C.H., Hacker, B.R. & Engdahl, E.R., 2007. Global prevalence of double benioff zones, *Science*, **316**, 1472–1474.

* Campos, J. & Kausel, E., 1990. The large 1939 Intraplate earthquake of Southern Chile, *Seism. Res. Lett.*, **61**(1), 43.

* Carlson, J.M. & Batista, A.A., 1996. Constitutive relation for the friction between lubricated surfaces, *Phys. Rev. E*, **53**, 4153–4165.

* Cocco, M., Pacheco, J., Singh, S.K. & Courboulex, F., 1997. The Zihuatanejo, Mexico, earthquake of 1994 December 10 ($M = 6.6$): source characteristics and tectonic implications, *Geophys. J. Int.*, **131**, 135–145.

* Comte, D. & Suárez, G., 1994. An inverted double-seismic zone in Chile: evidence of phase transformation in the subducted slab, *Science*, **263**, 212–215.
Coutant, O., 1990. Programme de Simulation Numérique AXITRA, Rapport LGIT, Université Joseph Fourier, Grenoble, France.

Delouis, B. & Legrand, D., 2007. Mw 7.8 Tarapacá intermediate depth earthquake of 13 June 2005 (northern Chile): fault plan identification and slip distribution by waveform inversion, *Geophys. Res. Lett.*, doi:10.1029/2006GL028193.

Di Carli, S., Francois-Holden, C., Peyrat, S. & Madariaga, R., 2008. Kinematic and dynamic inversion of the 2000 Tottori earthquake based on elliptical subfault approximations, *Geophys. Res. Abstr.*, **11**, Abstract EGU2009-11330.

Engdahl, E., van der Hilst, R. & Buland, R., 1998. Global teleseismic earthquake relocation with improved travel times and procedures for depth determination, *Bull. seism. Soc. Am.*, **88**(3), 722–743.

François-Holden, C., Di Carli, S., Sladen, A. & Madariaga, R., 2007. Non linear kinematic source inversion of the 2000 Tottori, Japan earthquake with variable rupture velocity, *Geophys. Res. Abstr.*, **9**, Abstract EGU2007-A07712.

Gardi, A., Lemoine, A., Madariaga, R. & Campos, J., 2006. Modeling of stress transfer in the Coquimbo region of central Chile, *J. geophys. Res.*, **111**, B04307, doi:10.1029/2004JB003440.

Graizer, V.M., 1979. Determination of the true ground displacement by using strong motion records, *Izvestiya, Earth Phys.*, **15**(12), 875–885.

Hacker, B.R., Peacock, S.M., Abers, G.A. & Holloway, S.D., 2003. Subduction factory 2. Are intermediate-depth earthquakes in subducting slabs linked to metamorphic dehydration reactions? *J. geophys. Res.*, **108**, doi:10.1029/2001JB001129.

Hanuš., V., Slancová., A., Špičák, A. & Vaněk, J., 1999. Discontinuous Nature of the Lower Part of the South American Wadati-Benioff Zone in The Arica Elbow Region, *Stud. Geophys. Geod.*, **43**, 163–184.

Hartzell, S., Liu, P., Mendoza, C., Ji, C. & Larson, K.M., 2007. Stability and uncertainty of finite-fault slip inversions; application to the 2004 Parkfield, California, *Bull. seism. Soc. Am.*, **97**(6), 1911–1934.

Ide, S. & Takeo M., 1996. The dynamic rupture process of the 1993 Kushiro-Oki earthquake, *J. geophys. Res.*, **101**, 5661–5675.

Iwan, W.D., Moser, M.A. & Peng, C.Y., 1984. Strong-motion earthquake measurement using a digital accelerograph, *Bull. seism. Soc. Am.*, **75**, 1225–1246.

Jiao, W., Silver, P.G., Fei, Y. & Prewitt, C.T., 2000. Do intermediate- and deep-focus earthquakes occur on preexisting weak zones? An examination of the Tonga subduction zone, *J. geophys. Res.*, **105**(B12), 28 125–28 138.

Jung, H., Green II, H. & Dobrzhinetskaya, L., 2004. Intermediate-depth earthquake faulting by dehydration embrittlement with negative volume change, *Nature*, **428**, 545–549.

Kausel, E. & Campos, J., 1992. The Ms = 8 tensional earthquake of 9 December 1950 northern Chile and its relation to the seismic potential of the region, *Phys. Earth planet. Inter.*, **72**, 220–235.

Kirby, S., 1995. Intraslab earthquakes and phase changes in subducting lithosphere, *Rev. Geophys.*, **33**(Suppl.), 287–297.

Kosuga, M., Sato, T., Hasegawa, A., Matsuzawa, T., Suzuki, S. & Motoya, Y., 1996. Spatial distribution of intermediate-depth earthquakes with hor-

izontal or vertical nodal planes beneath northeastern Japan, *Phys. Earth planet. Inter.*, **93**, 63–89.

Lemoine, A., Madariaga, R. & Campos, J., 2002. Slab-pull and slab-push earthquakes in the Mexican, Chilean and Peruvian subduction zones, *Phys. Earth planet. Inter.*, **132**, 157–175.

Mai, P., Burjanek, J., Delouis, B., Festa, G., François-Holden, C., Monelli, D., Uchide, T. & Zahradnik, J., 2007, *EOS, Trans. Am. geophys. Un.*, **88**(52), Fall Meet. Suppl., Abstract S53C-08.

Mikumo, T., Yagi, Y., Singh, S.K. & Santoyo, M.A., 2002. Coseismic and postseismic stress changes in a subducting plate: possible stress interactions between large interplate thrust and intraplate normal-faulting earthquakes, *J. geophys. Res.*, **107**, doi:10.1029/2001JB000446.

Monelli, D. & Mai, P.M., 2008. Bayesian inference of kinematic earthquake rupture parameters through fitting of strong motion data, *Geophys. J. Int.*, **173**, 220–232, doi:10.1111/j.1365-246X.2008.03733.x.

Nielsen, S.B., Carlson, J.M. & Olsen, K.B., 2000. Influence of friction and fault geometry on earthquake rupture, *J. geophys. Res.*, **105**, 6069–6088.

Nielsen, S.B. & Carlson, J.M., 2000. Rupture pulse characterization: self-healing, self-similar, expanding solutions in a continuum model of fault dynamics, *Bull. seism. Soc. Am.*, **90**, 1480–1497.

Ozel, N. & Moriya, T., 1999. Different stress directions in the aftershocks focal mechanisms of the Kushiro-Oki earthquake of Jan. 15, 1993, SE Hokkaido, Japan, and horizontal rupture in the double seismic zone, *Tectonophysics*, **313**, 307–327.

Peacock, S., 2001. Are the lower planes of double seismic zones caused by serpentine dehydration in subducting oceanic mantle? *Geology*, **29**, 299–302.

Peyrat, S. et al., Tarapacá intermediate-depth earthquake (Mw 7.7 2005, northern Chile): a slab-pull event with horizontal fault plane constrained from siesmologic and geodetic observations, *Geophys. Res. Lett.*, **33**, L22308, doi:10.1029/2006GL027710.

Rietbrock, A. & Waldhauser, F., 2004. A narrowly spaced double-seismic zone in the subducting Nazca plate, *Geophys. Res. Lett.*, **31**, doi:10.1029/2004GL019610.

Rietbrock, A., Haberland, C. & Nippress, S., 2006. A tear in the subducting Nazca slab at 21° S revealed from accurate locations of intermediate depth seismicity, *EOS, Trans. Am. geophys. Un.*, **87**(52), Fall Meet. Suppl., Abstract S43D-08.

Sambridge, M., 1999a. Geophysical inversion with neigbourhood algorythm—I: searching a parameter space, *Geophys. J. Int.*, **138**, 479–464.

Sambridge, M., 1999b. Geophysical inversion with neigbourhood algorythm—II: appraising the ensemble, *Geophys. J. Int.*, **138**, 727–746.

Stauder, W., 1973. Mechanism of spatial distribution of chilean earthquakes with relation to subduction of the oceanic plate, *J. geophys. Res.*, **78** (23), 5033–5061.

Vallée, M. & Bouchon M., 2004. Imaging coseismic rupture in far field by slip patches, *Geophys. J. Int.*, **156**, 615–630.

Warren, L.M., Langstaff, M.A. & Silver, P.G., 2008. Fault plane orientations of intermediate-depth earthquakes in the Middle America Trench, *J. geophys. Res.*, **113**, B01304, doi:10.1029/2007JB005028.