G J I S e i s m o l og y A synthetic seismicity model for the Xianshuihe fault, southwestern China: simulation using a rate- and state-dependent friction law

We performed a numerical simulation of seismic activity along the Xianshuihe fault, a highly active strike-slip fault approximately 350 km long, located along the eastern margin area of the Tibetan plateau, southwestern China. Historical earthquake data over the last 300 yr indicate repeated periods of seismic activity, and migration of large earthquakes along the fault during active seismic periods. To understand the characteristics of historical seismicity and to obtain insight into seismic potential, we performed a numerical simulation of slip behaviour along the fault using a 2-D elastic plate model. The friction on the model fault obeys a laboratory-derived rate- and state-dependent friction law, while the long-term slip rates along the fault are assumed to be consistent with geologically and geodetically estimated slip rates. To simulate segmented rupture behaviour we introduced non-uniformity into friction parameters along the fault. The simulation results indicate that seismic rupture is arrested in regions with frictional properties that are highly velocity-strengthening or large values of characteristic slip distance over which frictional stress evolves. In the regions where seismic rupture is arrested, post-seismic sliding occurs, causing time-dependent stress transfer along the fault. Post-seismic slip histories in highly velocity-strengthening regions are well approximated by logarithmic time functions, while those in regions with large characteristic slip distances increase more rapidly in the initial stages of post-seismic sliding. This difference in post-seismic sliding produces the difference in time-dependent stress transfer and, in consequence, the statistical characteristics of earthquake recurrence. Fault interactions near a fault branch point produce complexity in the simulated earthquake sequence. Simulated rupture commonly triggers additional rupture in a neighbouring segment with a time delay of several years. The simulated seismic activity shows good agreement with observational data with respect to the following statistical features: (1) successive seismically active periods are separated by periods of relative quiescence of 100– 200 yr duration (2) and the duration of active period is a few tens of years. Strictly speaking, neither the time- nor slip-predictable models are able to describe the simulated local slip history for each fault segment. The ratio of the standard deviation to the average value of the recurrence intervals of simulated earthquakes ranges from 0.05 to 0.22, depending on the fault segment. The cumulative distribution function of the recurrence intervals of simulated earthquakes at each segment approximately obeys a Brownian passage time distribution or a lognormal distribution, which is commonly used in the statistical evaluation of earthquake occurrence.

. Map showing the Xianshuihe fault and related faults in the Sichuan-Yunnan block, southwest China. Legend key: 1, active strike-slip fault; 2, active normal fault; 3, active reverse fault; 4, major (thick) and secondary (thin) faults and 5, national border. Inset map shows the location of Fig. 1 in terms of continental China and its main active faults. Modified from Wen et al. (2002). related to collision of the Indian and Eurasian plates (e.g. Molnar & Tapponnier 1975). Previous studies have described up to 15 mm yr −1 of left-lateral slip on the Xianshuihe fault, as determined from geological criteria (Wen et al. 1988;Allen et al. 1991) and GPS observations (King et al. 1997). Approximately 20 earthquakes of M ≥ 6 have occurred along the fault since 1700 (Department of Earthquake Disaster Prevention, Sate Seismological Bureau 1995; Department of Earthquake Disaster Prevention, China Seismological Bureau 1999). The rupture areas of past large earthquakes along the Xianshuihe fault indicate that the fault is segmented and that large earthquakes have occurred repeatedly on each segment (e.g. Allen et al. 1991;Wen et al. 2002). The fault records alternating periods of seismic activity and quiescence (e.g. Huang & Yang 1988;Ma et al. 1989;Wen et al. 2002), and seismicity appears to migrate along the fault during active periods (e.g. Yi & Wen 2000;Papadimitrou et al. 2004). In the seismically active period from 1725 to 1816, earthquakes occurred in most of the fault segments (e.g. Yi & Wen 2000). The most recent earthquake sequence began in 1893 and ruptured many segments. The sequence includes the 1955 Kanding (M = 7 1/2), 1973 Luhuo (M = 7.6) and 1981 Daofu earthquakes (M = 6.9). The source processes of the 1973 and 1981 earthquakes were investigated by modelling body and surface waves to show that the focal mechanisms are consistent with the surface trend of the fault (Zhou et al. 1983a,b). As some fault segments remain unbroken during the recent earthquake sequence, several studies have highlighted the existence of seismic gaps along the Xianshuihe fault (e.g. Allen et al. 1991). Papadimitrou et al. (2004) calculated Coulomb stress changes due to large earthquakes along the Xianshuihe fault, and found that the observed migration of seismicity can be elucidated by coseismic static stress changes and steady tectonic loading. Furthermore, Papadimitrou et al. (2004) calculated stresses within the current seismic gaps to evaluate seismic potential in these areas. Wen et al. (2002) used historical earthquake recurrence data and the empirically determined distribution of earthquake recurrence time to calculate the probability of earthquake occurrence along individual fault segments.
Migration of seismicity during a seismically active period has been observed also for some large active faults (e.g. Ambraseys 1970;Wood & Allen 1973;Rydelek & Sacks 2001). Stein et al. (1997) showed that the migration of seismicity along the North Anatolian fault, Turkey, can be elucidated by the positive changes in Coulomb failure stress. Since the migration of earthquakes along long geological faults may be a common phenomenon, understanding its mechanism is important for evaluating the potential of earthquake occurrence.
Synthetic models of seismicity are useful for understanding the mechanics of fault sliding behaviour, such as the complexity of seismic cycles and fault interaction (e.g. Ward 1991Ward , 1996Cornell et al. 1993). Such models are also useful for investigating the statistical characteristics of earthquake recurrence, because the recorded earthquake history is insufficient for statistical analysis (e.g. Robinson 2004). Many previous studies of synthetic seismicity have assumed simple frictional properties for model faults, such as constant static/dynamic friction and velocity-weakening friction. On the basis of several kinds of experimental results, Dieterich (1979) and Ruina (1983) derived rate-and state-dependent friction laws. As the friction laws properly take into consideration the rate-and slipdependence of frictional strength and time-dependent restrengthening of faults, they can explain various experimental observations, such as the dynamic behaviour of unstable slip, the occurrence of aseismic sliding and the slip nucleation process (e.g. Tullis 1988;Scholz 1990;Marone 1998). Rate-and state-dependent friction laws can simulate much more realistic and complicated sliding behaviour than simpler friction laws; they are therefore, commonly used to model seismic cycles at relatively simple plate boundaries (e.g. Tse & Rice 1986;Rice 1993;Stuart & Tullis 1995;Kato & Hirasawa 1999). Kato (2001) presented a rate-and state-dependent friction model for simulating seismic cycles at intersecting reverse faults, and examined the mechanical interaction of the faults and the statistical characteristics of simulated earthquake recurrence. For numerical simulations of seismic cycles along strike-slip faults using rate-and state-dependent friction laws, Ben-Zion & Rice (1995) introduced depth-dependent pore pressure distribution to generate complex earthquake cycles and Hillers et al. (2006) examined the effect of spatial variations of characteristic slip distance on complex slip pattern and earthquake statistics. However, no previous study has attempted a numerical simulation with rate-and state-dependent friction laws for seismic cycles on long segmented strike-slip faults while assuming realistic geometries of existing geological faults, despite the fact that such a model may provide an important insight into the mechanics of complex earthquake sequences.
It should be noted that the rate-and state-dependent friction laws were derived from experiments at sliding velocities much lower than seismic slip velocities. Frictional behaviour during seismic slip might be significantly different from that expected from the rate-and state-dependent laws, due to several possible factors, such as thermal pressurization (e.g. Mase & Smith 1987) and frictional melting (e.g. Tsutsumi & Shimamoto 1997). At present, it is dif-ficult to conduct numerical simulations of seismic cycles for long segmented faults that take into consideration such complicated processes during seismic slip.
In the study reported herein, we employed a rate-and statedependent friction law within a numerical simulation of seismic cycles along a bent, branched and segmented strike-slip model fault that is similar in geometry to the Xianshuihe fault. Model parameters were chosen such that the simulated seismic cycles mimicked historical seismic activity along the Xianshuihe fault. On the basis of the numerical simulation, we discuss the mechanics of the segmentation of seismic activity, migration of seismic activity and statistical characteristics of the recurrence of the simulated earthquakes.

N U M E R I C A L M O D E L
We considered a 2-D elastic plane-stress model for the Xianshuihe fault (Fig. 2), assuming for the sake of simplicity that the shallower seismogenic elastic layer is decoupled from the deeper ductile substratum. The hypocentre locations of moderate and minor earthquakes along the Xianshuihe fault indicate that the seismogenic layer extends to a depth of approximately 30 km (Ma et al. 1989). This seismogenic zone thickness of 30 km is much less than the fault length of approximately 350 km; plane-stress modelling is, therefore, justified. Since we used a thin elastic plate model, the longrange elastic interaction was taken into consideration. However, we did not consider the effects on stress variation of the variation of aseismic sliding and viscoelastic deformation at deeper aseismic regions, though these effects may be important for time-dependent stress transfer (e.g. Rydelek & Sacks 2001). These effects should be considered in future studies.
The rigidity G and Poisson's ratio ν of the model elastic plate were 30 GPa and 0.25, respectively. The model fault geometry for the Xianshuihe fault shown in Fig. 2 was composed of 1752 straight line cells, each 250 m long. Static shear stress acting on the ith cell at time t due to tectonic loading and fault slip is described as where K ij is static shear stress at the centre of the ith cell due to uniform unit slip at the jth cell, as given by elastic dislocation theory (e.g. Hirth & Lothe 1982), u j is left-lateral slip at the jth cell, and V S j is the long-term tectonic slip rate at the jth cell. In calculating K ij , the model fault geometry (Fig. 2) is taken into consideration. In eq. (1), we assume that the shear stress at each cell is constant if the slip rates at all cells are equal to the long-term tectonic slip rates. We did not consider the perturbation of normal stress applied to the fault plane due to slip. High concentrations of normal stress will actually develop near fault bends, due to cumulative slip in pure elastic models; these stress concentration are not relaxed by fault slip and lead to the breakdown of the numerical simulation. In the real crust, local inelastic deformation, such as secondary faulting and microfractures, relax stress concentrations. Hence, neglecting the concentration of normal stress around fault bends is considered to be a good approximation in numerical simulation. As we simulated the entire seismic cycle, including coseismic slip, we used the radiation damping approximation (Rice 1993) to evaluate shear stress reduction due to high-speed seismic rupture. The shear-stress at the ith cell is then given by where V j (= du j /dt) is the left-lateral slip rate at the jth cell and β is S-wave velocity (3.27 km s −1 ). The frictional stress on the fault was assumed to follow a composite rate-and state-dependent friction law (Kato & Tullis 2001. The slip law and the slowness law are the most popular versions of the rate-and state-dependent friction law (e.g. Nakatani 2001). However, it is known that the slip law cannot quantitatively explain experiments of strengthening at very low sliding velocities and the slowness law fails to fit the experimental frictional behaviour around steady-state stable sliding. The composite law explains these experimental observations and is probably more appropriate for modelling seismic cycles. According to the composite law, τ is expressed by where V is the sliding speed, θ is a state variable, and A, B, L, τ 0 , V 0 and V c are constants. τ 0 and V 0 are reference stress and velocity, respectively; they do not affect simulated sliding behaviour. We used V c = 0.01 μm s −1 , following Kato & Tullis (2003). For A − B < 0, the steady-state frictional stress decreases with an increase in velocity, possibly leading to unstable (seismic) slip. Conversely, the steady-state frictional stress increases with velocity for A − B > 0. This acts to stabilize sliding behaviour and inhibits the generation of unstable slip. Ruina (1983) showed theoretically that frictional sliding becomes increasingly unstable for larger positive values of (B − A)/L. Fig. 3 shows the long-term slip rate V S and the friction parameters A, B, A − B and L assumed along the model fault in the present study, where Fault 1 denotes the main fault and Faults 2 and 3 are branches as shown in Fig. 2. The frictional parameters shown in Fig. 3 correspond to the case that shows the best agreement between simulation results and observation data (see the next section for details). The long-term slip rates used in the present model are approximately the same as those of Wen et al. (2002), who compiled slip rates determined from geological data and GPS observations. These assumed slip rates are also consistent with values recently estimated by Iwakuni (2004) mainly on the basis of GPS data. We divided the model fault into segments S1-S9, as shown in Figs 2 and 3, following Wen et al. (2002). The mechanical origin of fault segmentation is probably related in part to geometric structures such as fault bends and branches. However, such geometric features are not apparent at all segment boundaries, which indicates spatial non-uniformity of stress or frictional properties along the fault. Yi et al. (2004) examined the activity of small earthquakes along the Anninghe-Zemnhe fault, which is located at the southern extension of the Xianshuihe fault ( Fig. 1), and found non-uniformity in the activity of small earthquakes and in the b-value of the Gutenberg-Richter relation, which supports the contention that non-uniformity of stress or frictional properties are present along the fault. In addition, significant post-seismic sliding was observed at certain locations along the Xianshuihe fault following the 1973 M = 7.6 Luhuo earthquake, while little post-seismic sliding was found in regions of large coseismic slip (Allen et al. 1991;Bilham 1992;Ge 1992). The occurrence of post-seismic sliding usually indicates velocity-strengthening (A − B > 0) frictional properties (Marone et al. 1991), while seismic slip indicates velocity-weakening (A − B < 0) frictional properties. Aseismic sliding also occurs in regions that have very large L values, as indicated by Hillers et al. (2006). This fact also supports the existence of non-uniform frictional properties along the Xianshuihe fault.
We performed several tens of simulation runs by varying the spatial distribution of friction parameters along the model fault to examine the effect of non-uniform frictional properties on the complexity of earthquake cycles. Furthermore, we tried to find a distribution of frictional parameters that can express the characteristics of earthquake cycles along the Xianshuihe fault through trial and error. The characteristics for which we tried to fit the simulation to observations (e.g. Allen et al. 1991;Wen et al. 2002) are as follows.
(1) Seismically active periods occurred repeatedly at intervals of ∼200 yr, and M7 class earthquakes occurred at some segments during each active period.
(2) Each active period lasted a few tens of years. Since the recurrence intervals and the magnitudes of simulated earthquakes depend on A − B and L values in each segment (e.g. Stuart 1988), we varied these values to satisfy (1). As will be discussed later, A − B and L values in segment boundary regions affect the durations of seismically active periods. We therefore, varied these values to satisfy (2).
Of the many simulation runs, we will discuss six cases (see below). For the six cases, we assumed large positive A − B value of (A − B) 1 and/or large L value of L 1 around the segment boundaries between S1 and S2, between S2 and S3 and between S3 and S4 (Figs 3b and c), because a region with A − B > 0 or large L acts as a barrier to rupture propagation (Kato & Hirasawa 1999). The A value was fixed at 0.4 MPa except for the segment boundary regions. In the central part of each segment where seismic slip is nucleated, A − B was assumed to be a negative value (A − B) 2 = −0.08 MPa. This value of −0.08 MPa at seismogenic regions is consistent with laboratory data of granite friction at hydrothermal conditions (Blanpied et al. 1995). In the remaining regions, A − B was assigned a small positive value (A − B) 0 = 5 × 10 −5 MPa. Note that A − B is proportional to the effective normal stress, which is given by the normal stress applied to the fault minus the pore fluid pressure. Since the pore pressure at hypocentral depths is not known, large uncertainty exists concerning A − B values for real faults. The L value was assumed to be constant at L 0 = 0.03 m, except for the segment boundary regions where L = L 1 (Fig. 3c). Laboratory-derived values of L are typically 1-100 μm, and experiments show that L increases with the characteristic scale of roughness of sliding surfaces or the thickness of the fault gouge layer. L values of natural faults, are therefore, thought to be much larger than laboratory values (e.g. Marone & Kilgore 1993). The L 0 values of 0.03 m assumed in the present model are similar to those of previous simulation studies of earthquake cycles that use rate-and state-dependent friction laws (e.g. Tse & Rice 1986;Stuart & Tullis 1995;Kato 2001). For the six cases discussed in the present paper, only (A − B) 1 and L 1 in the segment boundaries (Figs 3b and c) Table 1. (A − B) 1 and L 1 values in the regions between segments S1 and S2, between S2 and S3 and between S3 and S4 for cases 1-6. See were variable. The values of (A − B) 1 and L 1 are shown in Table 1.
In the segment boundary regions between S1 and S2, between S2 and S3 and between S3 and S4, B was assumed to be 0 and (A − B) 1 = A for cases 1-3. Positive A − B values with B = 0 were assumed, to represent creeping sections (e.g. Stuart & Tullis 1995) possibly due to the existence of serpentine (Reinen et al. 1994). A = 0.4 MPa all over the fault for cases 4-6. In case 1, we assumed We consider that the simulation result for case 1 approximates relatively well the characteristics of observed seismic cycles along the Xianshuihe fault. Accordingly, we discuss mainly the simulation results for case 1. It should be noted that the simulation result for case 1 does not fit the observations exactly and that similar or better pattern of seismic cycles may be obtained for other sets of spatial distributions of frictional parameters. It cannot be concluded that the spatial distribution of frictional parameters along the actual Xianshuihe fault was determined by the study reported herein.
For friction that follows the rate-and state-dependent law, a critical slip nucleation zone size h * exists (Dieterich 1992). For numerical simulations, the computation cell size h must be much smaller than h * to avoid artificial slip instability (Rice 1993). In the case of plane stress shear crack, h * is given by In the simulation reported herein, h/h * < 0.02, which guarantees appropriate numerical results.
An initial condition of the simulation was that V init = 0.01 m yr −1 and θ = L/V init for the entire fault. Eqs (2)-(4) were computed numerically using a Runge-Kutta method with an adaptive time step control to maintain numerical accuracy (Press et al. 1992) to obtain a record of slip and shear stress over a period of 24 678 yr for case 1. This time interval is long enough to provide sufficient numbers of simulated earthquakes for valid statistical analyses of recurrence intervals. The numbers of the computation time steps were the same for the six cases. This results in some difference in the period of simulated histories of slip and stress, because the computation time step is variable.

Slip histories
A portion of the simulated slip histories reduced by constant-rate long-term slip, u − V s t, for case 1 is shown in Fig. 4; the locations of the observation points on the fault are points 1-19 in Fig. 3. When a fault is locked (V = du/dt = 0), reduced slip u − V s t decreases with time at a constant rate V s . A sudden increase in u − V s t indicates an earthquake. No jerky slip motion is observed at points 3, 5 and 7 in Fig. 4; aseismic sliding is expected at these points because A − B > 0. Apparent stick-slip motion is observed at some points (e.g. points 9 and 10) with A − B > 0 in Fig. 4. The apparent sudden slip is rapid post-seismic slip rather than seismic, as will be discussed below. Fig. 4    commonly occur in clusters. Seismically active periods, where slip occurs at many points within a few tens of years, occur repeatedly, and successive active periods are separated by quiescent periods lasting 100-200 yr. We measured the recurrence interval of simulated earthquakes in each segment to obtain the average value, standard deviation and the coefficient of variation (COV), which represents the standard deviation divided by the average value and will be discussed in 3.2, as tabulated in Table 2. The recurrence intervals of earthquakes in segment S6 were not measured, because slip behaviour in the northern part of S6 is different from that in the southern part, as seen in the simulated reduced slip histories at points 10 and 11 in Fig. 4. As slip occurs simultaneously in segments S7 and S9 (points 12-15 in Fig. 4), the recurrence intervals of earthquakes in those segments were combined. The recurrence intervals of simulated earthquakes in the northern part of Fault 1 (segments S1-S4) are relatively short (∼200 yr), while intervals for the southern part of Fault 1 (S7 and S9) are relatively long (∼300 yr), and those for Faults 2 and 3 (S5 and S8) are approximately 500 yr (Table 2). These simulation results are roughly consistent with the recurrence intervals of past earthquakes on the Xianshuihe fault, as estimated from historical and geological data (e.g. Wen et al. 2002).
To investigate detailed sliding behaviour during periods of active seismicity, we examined simulated slip histories on expanded timescales (Figs 5 and 6), which correspond to the active periods A and C in Fig. 4. In active period A, rupture first occurs in segment S1, at the northern end of Fault 1. Approximately 30 yr later, rupture occurs in segment S4, which is approximately 100 km from segment S1; within 3 yr, the rupture migrates bilaterally to S3 and S2 in the north and S6, S7 and S9 in the south. The total duration of the active period is about 30 yr. Fig. 6 shows simulated slip histories during the active period C in Fig. 4. Rupture starts in segment S3 and then jumps to S6, S7 and S9. Over the total duration of the active period of about 40 yr, seismic activity migrates to segments S8, S4, S1 and S2.
The above examples of simulated slip histories during active periods indicate significant variation in the duration of an active period, the number of ruptured segments, the order of segment rupture and the time interval between successive ruptures. This variation probably occurs because the stress at each fault segment at any time depends on the complicated slip history along the fault. Fig. 7 shows the locations of the initiation of earthquake rupture of 11 simulated earthquakes in segment S4 for the time interval shown in Fig. 4  velocity first exceeds 10 mm s −1 during each simulated earthquake. The point of initiation of earthquake rupture is always located in the region of velocity-weakening friction (A − B < 0), consistent with the slip instability condition (Fig. 7; e.g. Ruina 1983). The location of earthquake rupture initiation is variable, probably dependent on complicated past slip histories and rupture tends to initiate close to neighbouring segments that housed preceding ruptures. Many of the ruptures shown in Fig. 7 were preceded by ruptures in segment S6. This indicates that repeated ruptures within individual segments are not necessarily identical with respect to the rupture process. This observation is supported by the variation in coseismic slip of simulated earthquakes in each segment (Fig. 4). It should be noted that frictional properties and normal stress are uniform in the velocity-weakening region of each segment. Results of the present simulation do not preclude the existence of characteristic locations where earthquake rupture always initiates given non-uniformity in frictional properties or stress.
In general, the duration of active periods in our model of seismicity is a few tens of years, and rupture occurs within segments that experience shear stress elevated by the effect of preceding ruptures, consistent with previous observations of the Xianshuihe fault (Papadimitrou et al. 2004). Delayed propagation of rupture is caused by time-dependent stress transfer due to post-seismic sliding in the velocity-strengthening region (A − B > 0), as demonstrated by Kato (2004). Fig. 8(a) shows post-seismic sliding, following the earthquake in segment S2 during the time interval shown in Fig. 6, at locations 5, 10, 15, 20 and 25 km from the southern end of the velocity-weakening region in segment S2. The maximum amplitude of post-seismic sliding is recorded near the end of the seismogenic zone with A − B < 0, and post-seismic slip amplitude decreases with distance from the seismogenic zone. At all locations, the rate of post-seismic slip decreases over time. Although the apparent time functions of simulated post-seismic sliding are different at various locations, they are well expressed by where t is time measured from earthquake occurrence, and u 0 , t 0 , p and q are constants. q ≈ 0 for all areas except the location 5 km from the velocity-weakening region in segment S2, where a rapid increase in slip was recorded immediately following the earthquake. Eq. (5) with q = 0 is the theoretical time function of post-seismic sliding for a simple spring-block model with rate-and state-dependent friction (Scholz 1990;Marone et al. 1991). Post-seismic sliding similar to that in Fig. 8(a) was commonly observed in velocity-strengthening regions with large positive A − B values in the present simulation. Significant post-seismic sliding was also observed in the simulation for large L in segment boundary regions for cases 4 and 5. Aseismic sliding is expected to occur for large L (e.g. Hillers et al. 2006). Fig. 8(b) shows simulated slip histories, following a simulated earthquake in segment S2, at the same locations on the model fault as those in Fig. 8(a), together with fitted logarithmic functions eq. (5). Post-seismic sliding caused by large L approximately obeys the logarithmic time function as well as that caused by large positive A − B. Comparing Figs 8(a) and (b), we find more rapid slip increase just after an earthquake for case 5. The p-value in eq. (5) should be positive in order to represent plausible slip histories over a long period; otherwise, the slip rate would be negative for large t. The determined p-values for five fitted functions in Fig. 8(a) for case 1 ranged from 2 mm yr −1 to 6 mm yr −1 , which are a little smaller than the assumed long-term slip rate of 10 mm yr −1 . However, negative p-values were obtained for the locations at 15, 20 and 25 km from the A − B < 0 region in Fig. 8(b) for case 5. This suggests that the logarithmic function eq. (5) may not be a very good approximation for post-seismic sliding, due to large L.
Since the 1973 Luhuo earthquake (M = 7.6) in segment S1 and part of S2, significant post-seismic sliding has been observed at Xialatuo, about 10 km SE of Luhuo (Fig. 1)   the source area of the earthquake (Allen et al. 1991;Ge 1992;Bilham 1992). Fig. 9 shows post-seismic sliding at Xialatuo (data from the Earthquake Administration of Sichuan Province, China) estimated from the time variation of a 133 m long baseline A − B (Fig. 9). This site was established in the mid-1970s and designed to cross the surface rupture of the 1973 M7.6 earthquake. Accordingly, postseismic data for about 5 yr immediately following the earthquake cannot be used. Cumulative sliding at this site following the 1973 earthquake has been observed as a steady and gradual elongation of the baseline A − B, which indicates that fault movement is leftlateral and aseismic. The general rate of post-seismic sliding has been slowing gradually, and the observed curve is well fitted by eq. (5), though the data just after the earthquake were not obtained. The time intervals between successive simulated earthquakes for case 1 may be shorter than those of historical earthquakes along the Xianshuihe fault; it is, therefore, also possible that the duration of simulated active periods, of a few tens of years, is shorter than observed values. For example, an active period occurred from 1725 to 1816 (e.g. Yi & Wen 2000). This discrepancy in the duration of seismically active periods might be improved by changing friction parameters along the model fault. The migration speed would then be lower for large positive A − B values in the velocity-strengthening region between seismogenic zones (Kato 2004).

Complex earthquake recurrence patterns
The data shown in Fig. 4 for case 1 indicate that relatively regular stick-slip behaviour has occurred in the northern (points 1, 2 and 4 in segments S1 and S2) and southern parts (points 12-15 in S7 and S9) of Fault 1 and in Fault 3 (points 18 and 19 in S8). The COV (= the standard deviation divided by the average value) of the recurrence intervals of simulated earthquakes are relatively small in these segments (Table 2). Irregular stick-slip cycles with significantly variable recurrence intervals (large COV values) and coseismic slip amounts are observed in the middle part of Fault 1 (points 6, 8 and 9 in segments S3 and S4) and in Fault 2 (points 16 and 17 in S5). This variability in the irregularity of stick-slip cycles arises from the difference in the degree of fault interaction. In the middle part of Fault 1 and Fault 2, stress perturbation, due to fault slip in the northern and the southern part of Fault 1 and Fault 3, affects slip behaviour significantly. In addition, fault interaction at the fault branch point near points 10 and 16 is likely to have produced complex stress histories. Slip behaviour is variable during seismically active periods along segments S4, S5 and S6, which connect to the fault branch point. During some seismically active periods, such as A and C in Fig. 4, seismic slip in segments S4 and S6 does not propagate into S5. During other active periods, such as B in Fig. 4, seismic slip occurs in segments S4 and S5, but is absent in S6. This behaviour, in which only one of the subparallel (Fig. 2) segments S5 and S6 ruptures during an active period, is easy to understand because slip in one of the segments reduces shear stress in the other. However, in some seismically active periods, such as D in Fig. 4, seismic slip occurs along segments S4, S5 and S6. In these cases, the decrease in shear stress due to slip on subparallel faults in segments S5 and S6 is overcome by the effect of shear-stress increase related to slip along other segments, such as S4 and S7. It should be noted that the Ganzi-Yushu and Anninghe faults are adjacent to the Xianshuihe fault to the north and south, respectively. The two adjacent faults are left-lateral strike-slip faults, and have recorded large earthquakes (e.g. Wen et al. 2002). Slip-induced stress on these adjacent faults is expected to influence the slip behaviour of the Xianshuihe fault, especially at the northern and the southern ends of the fault. It is important that the relationships between these faults be investigated in future studies. Fig. 10 shows schematic slip histories reduced by constant-rate long-term slip for (a) periodic recurrence of characteristic earthquakes, (b) time-predictable earthquake recurrence, in which the occurrence time of the next earthquake can be predicted from the total slip of the preceding earthquake and (c) slip-predictable earthquake recurrence, in which the slip amount of the next earthquake can be predicted from the time elapsed from the preceding earthquake (Shimazaki & Nakata 1980). The time-predictable model implies that the fracture strength is constant, while the slip-predictable model implies that the residual frictional stress is constant provided that loading rate is constant. The simulated reduced slip histories at some points in Fig. 4 (case 1) seem to be approximated by the periodic recurrence model (Fig. 10a). However, strictly speaking, the simulated slip behaviour in many segments shown in Fig. 4 does not seem to obey periodic recurrence, the time-predictable model, or the slip-predictable model. This indicates that neither the fracture strength nor the residual frictional stress is constant.
In order to examine quantitatively recurrence pattern of simulated earthquakes, we measured local maxima s max and local minima s min of simulated reduced slip history at the central point of the velocityweakening (A − B < 0) region (thick line parts in Fig. 2) of each segment and determined their average values s max and s min and standard deviations δs max and δs min , respectively, as illustrated in Fig. 10(d).
The determined values are shown in Table 3, together with the average value u of coseismic slip of simulated earthquakes at the central point of the velocity-weakening region for each segment. δs max / u and δs min / u may represent irregularity of earthquake recurrence patterns; (a) δs max / u = δs min / u = 0 for the model of periodic recurrence of characteristic earthquakes (Fig. 10a), (b) δs min / u = 0 for the time-predictable model (Fig. 10b) and (c) δs max / u = 0 for the slip-predictable model (Fig. 10c). The values of δs max / u and δs min / u in Table 3 show that (1) both the values are much larger than 0 for each segment, (2) these values vary significantly dependent on the segment and (3) δs max / u and δs min / u take similar values for each segment. These results indicate that neither the time-predictable model nor the slip-predictable model explains the simulated earthquake recurrence patterns in the strict sense and that the irregularity of earthquake recurrence patterns depends on the segment.
On the basis of estimated slip amounts of historical earthquakes, Shimazaki & Nakata (1980) and Sykes & Quittmeyer (1981) suggested that the earthquake sequence in some regions can be approximated by time-predictable models, which have been used in the  Table 3. The average value of coseismic slip u of simulated earthquakes, the average value of local maxima of reduced slip s max just before an earthquake and its standard deviation δs max , the average value of local minima of reduced slip s min just after an earthquake and its standard deviation δs min and the ratios δs max / u and δs min / u obtained for the central point of the velocity-weakening (A − B < 0) region of each segment for case 1.  probabilistic forecasting of large earthquakes (e.g. Working Group on California Earthquake Probabilities 2003; Wen et al. 2002). However, detailed analyses of strain accumulation and release history along the San Andreas fault over the past 1500 yr suggest that neither the time-nor slip-predictable models can account for the estimated earthquake history (Weldon et al. 2004). On the basis of strain accumulation inferred from geodetic data, Murray & Segall (2002) concluded that the time-predictable model fails to explain the recurrence of M ≈ 6 earthquakes at Parkfield, California. Results of the present simulation indicate that the time-predictable model cannot be applied if significant fault interaction has occurred. Fig. 11 shows the recurrence interval of simulated earthquakes at segment S2 in relation to the sequence of seismic events (event number) for cases 1-6 (Table 1). These data provide no indication of a significant increase or decrease in the recurrence interval over time. The recurrence interval varies randomly for case 1; we cannot find any limit cycle oscillations. For case 2, the scatter of recurrence intervals changes significantly at event number ∼60, which indicates that the earthquake recurrence pattern changes at this point. For case 3, a similar change in the scatter of recurrence interval is found at event number ∼70, though it is less distinct than in case 2. For cases 4 to 6, where there are no significant barriers due to large positive A − B in segment boundary regions (Table 1), the variations of recurrence intervals of simulated earthquakes show patterns seemingly different from those for cases 1 to 3. The recurrence interval does not drift randomly, but exhibits some regularity, especially for cases 5 and 6. The recurrence intervals are about 200 yr for most events in case 5, and change abruptly to be short and then long. This pattern of abrupt change in recurrence interval appears repeatedly. Similar patterns are observed for case 6. Fig. 12 shows the time interval from the earthquake occurrence in segment S1 to that in segment S2 during each seismically active period, plotted against the event number for cases 1-6. The time interval ranges from −18.8 to 32.1 yr for case 1, and oscillates randomly between positive values (where rupture propagates from S1 to S2) and negative values (where rupture propagates from S2 to S1). This indicates that the simulated earthquake sequence is complex in nature. For cases 2 and 3, abrupt transitions of the time interval from negative to positive values and from positive to negative values are observed, respectively. The diagrams of the time interval versus the event number for cases 4-6 are very different from those for cases 1-3. The range of fluctuations of the recurrence interval is small for most periods, and there are some impulsive changes in the time interval for cases 4-6. It should be remarked that the unit in the time interval for case 6 (Fig. 12f) is the day, while for the other cases it is the year. The relatively small values of the time intervals for case 6 come from a weak barrier to rupture propagation in the segment boundary region between S1 and S2 because of small positive A − B and small L (Table 1). Comparing Figs 11 and 12, we find correspondence between the recurrence interval patterns for simulated earthquakes in S2 and the patterns of the time interval from earthquakes in S1 to S2. For cases 1-3, the scatter of the recurrence interval tends to be large for negative values of the time interval. For cases 5 and 6, impulsive changes in the recurrence interval correspond to impulsive changes in the delay time interval. This suggests that the earthquake recurrence pattern in S2 is significantly influenced by the interaction between S1 and S2.
The differences in (1) the recurrence patterns of simulated earthquakes ( Fig. 11) and (2) the time-dependent migration of seismic activity along the fault (Fig. 12) between the cases with large positive A − B barriers in the segment boundary regions (cases 1-3) and the cases without large positive A − B barriers (cases 4-6) probably arise from the difference in the time-dependent stress transfer in the segment boundary regions. As shown in Fig. 8, post-seismic slip histories for cases 1 and 5 differ. For case 1 with large positive A − B values in the segment boundary region, post-seismic slip continues for a long time with gradually decreasing sliding velocity, resulting in shear stress transfer for a long time. By contrast, for case 5 without large positive A − B values in the segment boundary region, post-seismic sliding velocity diminishes rapidly in comparison with case 1. This results in a greater chance of an earthquake trigger shortly after an earthquake in a neighbouring segment.
To examine variation in the recurrence intervals of simulated earthquakes, we obtained the cumulative distribution functions of the recurrence intervals of simulated earthquakes in S2, S4 and S8 for case 1 (black curves in Figs 13a-c) and for case 5 (black curves in Figs 13d-f). Although the COV varies significantly (Table 2), the distribution functional forms are similar for different fault segments for case 1. The red line in each panel in Fig. 13 is the fitted distribution function of a Brownian passage time distribution, which is derived by assuming (1) that the critical stress for earthquake occurrence is constant and (2) that the loading process is expressed by Brownian perturbation in addition to constant rate loading (Matthews et al. 2002). The Brownian passage time distribution provides a sound explanation of the distributions of the recurrence intervals of past earthquakes, and has been used in probability models for earthquakes in the United States of America (Working Group on California Earthquake Probabilities 2003) and Japan (Earthquake Research Committee 2001). The distributions of the recurrence intervals of simulated earthquakes for case 1 are well fitted also by the lognormal distributions (blue curves in Fig. 13), which is commonly used in probability models of earthquake occurrence (e.g. Nishenko & Buland 1987;Working Group on California Earthquake Probabilities 1990;Wen et al. 2002). The functional form of the lognormal distribution is very similar to the Brownian passage time distribution, and the distributions can explain the recurrence intervals of past earthquakes almost equally well. In Fig. 13, blue curves often overlap red curves. It is noted that the lognormal distribution has a defect, such that the probability of earthquake recurrence decreases with time when the elapsed time from the preceding earthquake is much longer than the average recurrence interval. Thus, it is not desirable to use it for calculating probabilities of earthquake occurrence (Matthews et al. 2002). The Weibull distribution has also been used in probability models for earthquake (e.g. Rikitake 1976), and can explain fairly well the distribution of recurrence intervals of simulated earthquakes for case 1 (green curves in Fig. 13). As discussed in Nishenko & Buland (1987) and Matthews et al. (2002), the Weibull distribution has heavy tails at short time intervals and the fit to simulated earthquakes is a little worse than with the lognormal distributions and the BPT distributions. It is interesting to note that (1) the fracture strength at each segment is not constant, due partly to the time-dependent property of friction used in the present model and (2)  The cumulative distribution functions of the recurrence intervals of simulated earthquakes for case 5 (Figs 13d-f) have features very different from those for case 1. They cannot be approximated well by the Brownian passage time distribution, the lognormal distribution, or the Weibull distribution, because there are prominent peaks in their frequency distributions, as understood from Fig. 11(e). This suggests that the different characteristics in time-dependent stress transfer affect the statistics of earthquake recurrence.
COV, which is also called aperiodicity (Matthews et al. 2002), ranges from 0.05 to 0.22 for case 1 in the present simulation ( Table 2). The COV values for the other cases do not differ greatly from those for case 1. These COV values are similar to or smaller than the estimated values for large earthquakes at plate boundaries of about 0.2 (Nishenko & Buland 1987), estimates of 0.2 for the Sichuan-Yunnan region, China (Wen et al. 2002), and estimates of 0.3-0.7 determined by the San Francisco Bay Region (Working Group on California Earthquake Probabilities 2003). This indicates that the degree of fault interaction in the present model is not large enough to explain the complicated seismicity recorded along real faults. Possible origins of the relatively minor complexity recorded in the present model are as follows. (1) We considered only ruptures of segments of several tens of km in length, and did not consider the effects of small earthquakes on and around the Xianshuihe fault.
(2) It is not possible for large earthquakes to occur over several segments in the present model, due to the existence of a long region of velocity-strengthening frictional property between velocityweakening regions. If the length of the velocity-strengthening region is short, a compound earthquake may occur over several regions with velocity-weakening frictional properties (Kato 2004).
In many earthquake recurrence models used to provide statistical forecasts of earthquake occurrence, the COV value is assumed to be independent of fault segment because the sample size of recurrence intervals for a single fault segment is insufficient (e.g. Nishenko & Buland 1987). However, results of the present simulation indicate that the COV of recurrence intervals varies significantly, dependent on the degree of fault interaction. Kumamoto & Hamada (2005) examined the COV values of recurrence intervals of large earthquakes along active faults in Japan and found that the COV value tends to increase with the number of neighbouring faults. Their conclusion is consistent with the present simulation result. For more reliable statistical forecast models of earthquakes, more appropriate procedures for estimating the variation in recurrence intervals should be developed.

D I S C U S S I O N A N D C O N C L U S I O N S
In the study reported herein, we performed numerical simulations of seismic activity along the Xianshuihe fault, southwestern China, assuming that (1) the long-term slip rates are consistent with geologically-and geodetically-derived values, (2) the friction on the fault obeys a rate-and state-dependent friction law and (3) nonuniformity of friction parameters along the fault is introduced to simulate segmentation of seismic activity. We performed a large number of simulation runs, and reported the simulation results for six of them. For the simulations, we found that seismic activity that is in good agreement with the historical earthquake sequence along the Xianshuihe fault over the past 300 yr can be obtained by selecting appropriate friction parameters (case 1). Important characteristics of simulated seismic activity that are consistent with observational data are as follows. (1) Successive seismically active periods are separated by periods of relative quiescence of 100 to 200 yr duration.
(2) The duration of an active period is a few tens of years. (3) The order of segment rupture along a fault is variable, and rupture occurs where shear stress has been increased by preceding earthquakes. We examined the distribution of the recurrence intervals of simulated earthquakes and found that the cumulative distribution function of the recurrence intervals in each segment approximately obeys a Brownian passage time distribution. This distribution is commonly used in probability models of earthquake occurrence, and the COV depends on the fault segment. The COV value tends to increase for fault segments with high fault interaction. This trend should be taken into account when constructing statistical models of earthquake recurrence, in order to forecast earthquakes more reliably.
We examined the effect of different frictional properties in segment boundary regions on post-seismic sliding, earthquake recurrence patterns, and the time interval between successive earthquakes in neighbouring segments. We found that post-seismic slip histories in regions with a highly velocity-strengthening region (large positive A − B values) are different from those in regions with a large characteristic slip distance (large L values). Post-seismic slip for large L increases rapidly in the initial stages, while slip for large positive A − B increases continuously for a long time. This causes a difference in time-dependent stress transfer due to propagating post-seismic sliding, which results in a difference in the patterns of recurrence intervals of simulated earthquakes and the delay time intervals between earthquakes in neighbouring segments. The simulated earthquake recurrence pattern sometimes changes abruptly in the order of segment rupture and the scatter of recurrence intervals. Similar observations of the abrupt change in the characteristics of earthquake pattern have been reported by Mitsui & Hirahara (2004) for their simulation with a model of multiple spring-blocks with rate-and state-dependent friction.
Complexity of earthquake recurrence can be simulated in numerical models in other studies assuming simple static/dynamic friction (e.g. Brown et al. 1991) or velocity-weakening friction (e.g. Carlson & Langer 1989). In these models, frictional properties on model fault planes are usually assumed to be uniform. However, the spatial variation of stress is produced spontaneously due to the dynamic effects of rupture (Cochard & Madariaga 1994), and a wide variety of earthquake sizes is simulated. These model studies suggest that stresses produced by many small earthquakes affect the occurrence of large earthquakes, and the models may be superior to the present model for representing earthquakes of different sizes. However, in our study, we intended to simulate sequences of large earthquakes along the Xianshuihe fault, where the segment structure of earthquake activity has been recognized. Accordingly, the introduction of non-uniform frictional properties can be justified for the present purpose. Zöller et al. (2005Zöller et al. ( , 2006 introduced creeping sections or non-uniform stress drop on the model fault planes with static/dynamic friction to investigate the frequency-size distribution of earthquakes. Such modelling may be applicable to simulations of seismicity along the Xianshuihe fault if the effect of smaller earthquakes is taken into consideration. In the present simulation, stress concentrations near a fault bend or branch point were not taken into consideration, as explained in Section 2. However, on real faults, a fault bend or branch sometimes controls the initiation or termination of earthquake rupture (e.g. King 1986). As velocity-strengthening frictional properties were assumed for all regions of fault bends and branch points except for the bend near point 11 (Fig. 2), earthquake rupture was arrested rather than initiated at these points. Stress concentrations near a fault bend or branch should relax over the long term due to inelastic deformation, as evidenced by the fact that horizontal deformation along the fault is transformed into continuous verti-cal deformation with increasing or decreasing volume. However, over a time interval shorter than the relaxation time of inelastic deformation, stress concentrations near a fault bend or branch may significantly affect the earthquake rupture process. These processes should be taken into consideration in future studies for more detailed models of earthquake rupture sequences on segmented strike-slip faults.
Finally, we present some suggestions concerning application to the mitigation of field hazards. The present simulation of a long time span provides statistical models of earthquake recurrence, which may be better than those obtained from records of past earthquakes over a very short-time period. The simulation result indicates that the Brownian passage time distribution can be used for the statistical model for earthquake occurrence, and the COV of the recurrence interval of large earthquake at a segment depends on the complexity of fault geometry. Different values of COV should be used for the statistical models. Measurement of post-seismic sliding may be useful for monitoring time-dependent stress transfer to unbroken fault segments.

A C K N O W L E D G M E N T S
We thank Shengli Ma and Changrong He for their support for the present cooperative research. We are grateful to Yehuda Ben-Zion (Editor) and two anonymous reviewers for valuable comments and suggestions, which greatly improved the manuscript. This research was funded by the Japan Society for the Promotion of Science (Project No. 14340128) and the National Basic Research Program of China (2004CB418405).