EVIDENCE FOR TEMPORAL CLUSTERING OF LARGE EARTHQUAKES IN THE WELLINGTON REGION FROM COMPUTER MODELS OF SEISMICITY

Temporal clustering of large earthquakes in the Wellington region. New Zealand. has been investigated with a computer model that generates long synthetic seismicity catalogues. The model includes the elastic interactions between faults. Faults included in the model. besides the subduction thrust between the Australian and Pacific plates. are segments of the four major strike-slip faults that overlie the plate interface (Wairarapa, Wellington. Ohariu, and Wairau faults). Parameters of the model are adjusted to reproduce the geologically ohserved slip rates of the strike-slip faults. The seismic slip rate of the subduction thrust, which is unknown, is taken as 25% of the maximum predicted by the plate tectonic convergence rate, and its position fixed according to recent geodetic results. For comparison. the model was rerun with the elastic interactions suppressed, corresponding to the usual approach in the calculation of seismic hazard where each fault is considered in isolation. Considering earthquakes of magnitude 7.2 or more ("characteristic" events in the sense that they rupture most of a fault plane). the number of short (0-3 years) inter-event times is much higher with interactions than for the corresponding case without interactions (46% vs. 2% or all inter-event times). This reduces to 9% vs. 2% if the subduction thrust is removed from the models. Paleoseismic studies of the past seismic behaviour of the subduction thrust are clearly needed if the degree of clustering is to be tightly constrained. Although some other aspects of our model can he improved in future. we think that the probability of significant short-term clustering of large events. normally neglected in hazard studies. is very high. This has important implications for the engineering. insurance and emergency response communities.


INTRODUCTION
The distribution of large earthquakes in time, within some region with many faults, is often assumed to be random.Statistical studies of world-wide earthquake catalogues have shown, however.that some degree of temporal clustering is common [4].But the historical record of cmthquakes is usually too short and inhomogeneous to make quantitative estimates of clustering in smaller regions.Paleoseismic studies can extend the record for some individual faults.but with relatively large uncertainties.Another approach to this problem is the development of synthetic (computer) seismicity models that can generate long, homogeneous catalogues of earthquakes [1,15,16).We have previously developed such a synthetic seismicity model [8] and applied it to the Wellington region I 9].Here we report on further results using the latest available information on fault geometry and slip rates of the major strike-slip faults in the region 112].and on the position of the seismically active section of the underlying subduction thrust [2], which in previous models was only a guess.In particular, we now include a recently recognised extension of the Ohariu fault further to the northeast than previously thought.
The assumption of randomness of large events arises from the usual practice in estimating regional hazard of considering that large "characteristic" earthquakes occur on independent faults more-or-less periodically.at a rate corresponding to their different long-term slip-rates.For a large number of faults this usually results in a distribution of regional inter-event times close to Poissonian, with no clustering.So it is important to realise that our synthetic scismicity model does not just generate earthquakes on independent faults, but includes the elastic interactions between the faults.When a large event occurs on a particular fault it not only produces a stress drop on that fault but it changes the stresses on all other faults in the region.The changes can either retard or advance the occurrence of events on those other faults.The sign and magnitude of the effect depends on the fault geometries, senses of slip, and various mechanical parameters.Only in very simple cases is it possible to predict the sign of the effect.and never its magnitude, without detailed calculations.For example (Figure l ), a large earthquake on one segment of a strike-slip fault would be expected to enhance the probability of events on adjoining segments along strike.Closely spaced, parallel strike-slip faults tend to be mutually inhibitory.However, if quantitative information is needed, or in more complex cases, detailed calculations are required.We think that physically realistic models, made possible by advances in computer power and knowledge of seismogenesis, can provide this information.In this paper the specific question we wish to answer is: Do fault interactions in the Wellington region increase the likelihood of multiple large earthquakes ( other than aftershocks) within a short time, compared to what would be calculated assuming independent faults, and, if so, by how much~ Or is the naive notion that a large earthquake has "released the stress", and that no nearby large events are likely to follow soon, a good guess~ We think that this question has important implications for the engineering and insurance industries.The short historical record for the W cllington region suggests that clustering of large events might well be expected: the four largest events occurred in 1848 (M 7.1), I 855 (M 8.1 ). and two in 1942 (M 7 .2,7 .0).

THE SYNTHETIC SEISMICITY MODEL
The theory behind our synthetic seismicity model has been described in detail elsewhere [8,9].Briefly, we consider each fault in the model as a single plane embedded in an elastic halfspace which is divided up into a large number of rectangular cells for which we specify several mechanical parameters ( coefficients of static and dynamic friction, stress-drop on failure, pore-pressure, etc).The faults can be of any orientation and sense of slip.Each fault is subjected to a "driving force" that slowly loads the cells towards failure.The stress, Tij, is monitored and when any cell fails (as determined by the Coulomb failure criterion, Tshear > Static strength) an "earthquake" begins.The amount of slip is determined by the specified stress drop.That initial failure induces stress changes on all other cells.These changes, calculated using dislocation theory [SJ, may immediately induce further cell failures, so leading to a larger earthquake; and so on, until eventually all cells are stable (Tshear < dynamic strength) and the earthquake is finished.Then the loading mechanism resumes.Any cell may fail more than once and each time the percentage stress drop is the same.This may lead to spatially variable stressdrop at the end of an earthquake.The number of ruptured cells and their displacements are used to calculate the earthquake's magnitude, leading in the encl to a synthetic catalogue of earthquakes.Note that induced stress changes arc calculated for cells in all faults, so it is possible for slip on one fault to initiate an earthquake on another fault.
For models with sufficiently small cell size our model generates a magnitude distribution similar to that for real scismicity (i.e., more small events than large ones, with a "bvalue" close to I).However, there is usually a distinct class of large "characteristic" earthquakes that rupture the whole fault and whose numbers are higher than expected from an extrapolation of the frequencies of smaller events, as is usually observed for well developed faults [ 11].The appearance of characteristic events does not come about due to the lack of smaller faults in the model, but is a natural result given reasonable mechanical properties.
In the model a cell's strength, S, is given by where, µ is the coefficient of friction (static or dynamic as appropriate), L is the normal stress due to lithostatic pressure, 8113 is the normal stress due to the loading mechanism and slip on other fault patches, Phase is some (unchanging) base pore pressure level (taken as a percentage of the lithostatic pressure).8P(t) is the change in pore pressure due to slip on other patches, t is time, and we have suppressed patch indices.Our sign convention is that tensile stresses are positive, but we retain the common usage of "pore pressure" (compression is positive).We assume that the medium is sufficiently permeable that the slow loading mechanism does not cause any changes in pore pressure.The induced changes in pore pressure are given by [10]: where f3 is Skempton's Coefficient which can range from Oto l depending on pore geometry relative to the stress changes.We assume that all induced changes in pore pressure decay exponentially with a time constant Tpore' and do not solve the equations for flow explicitly.
The specification of the driving mechanism and how it alters the stress on a fault with time can be clone in several ways, as discussed in our previous work.Various parameters are adjusted by trial-and-error until the long-term slip rates (and slip directions) predicted by the model approximate the geologically observed values.

APPLICATION TO THE WELLINGTON REGION
The Wellington region is one of east-west plate convergence and subduction of the Pacific plate beneath the Australian plate.The convergence, oblique to the strike of the major geologic structures, is thought to be partitioned onto arcnormal (i.e., NW-SE) and arc-parallel (NE-SW) components [13].The latter is taken up on a number of sub-paralleL dextraL strike-slip faults (Figure 2) that are known to have ruptured in large earthquakes in the past [ 12].The paleoseismic evidence also gives us values for the long-term slip rates and single event disphlcements for characteristic earthquakes on these faults [ 12].The arc-normal convergence is assumed to be taken up by slip on the plate interface and on imbricate thrust faults above it (eastern Wairarapa and offshore: also some on the Wairarapa fault).The position of the plate interface is well defined by the present day microscismicity [71 and recent geodetic surveys using GPS methods allow the identification of which part is currently "locked" and presumably slips in large earthquakes [2].However.no large subduction thrust earthquakes have occurred in historical times and the paleoseismic evidence for such events is unclear.Since such earthquakes :ire very common in other subduction zones it seems prudent to assume they also occur in the Wellington region.But it is uncertain how much of the arc-normal convergence we should assign to seismic slip on the plate interface.A value of 25%1 is our preferred guess (close to the median of values observed elsewhere), giving a recurrence time of large events of about 500 years.
We include in our model (Figure 3, Table I) the four major, upper-plate, strike-slip faults (Wairarapa, Wellington, Wairau.and taking the nearby Ohariu and Shepherd's Gully faults as one "Ohariu" fault), plus the subduction thrust.A cell size of approximately 5 x 5 km is used, a compromise between very small cells (most realistic but computationally expensive) and taking each fault as a single cell (not at all realistic but very fast computations).Such a cell size allows sufficient slip complexity to realistically model the interactions between closely spaced faults.We take each fault as a single plane that approximates the observed strike and dip.However, the

*
The two rates arc for along-strike and up-dip, respectively.in mm/yr: average over the fault plane for the subduction thrust.at the surface for the other faults.

** ***
A veragc recurrence time of earthquakes of magnitude 7 .2or more on that fault.Average slip for earthquakes of magnitude 7.2 or more.
because of relatively large side-steps near Te Marua and Waikanae, respectively.Because we cannot model the whole world, we limit the along-strike length or the region considered to 120 km, the approximate length of rupture on the Wairarapa fault during the great earthquake of 1855.
Mechanical fault properties have been set to averages of the many models (over 40) that we have previously investigated [9 J.These properties result in behaviour that includes some enhancement in the number of large "characteristic" earthquakes that rupture most of a fault plane over that expected from an extrapolation of the numbers of smaller events.We believe that the geological evidence and historical seismicity points to this sort of behaviour for the major Wellington region faults.
The models (with and without interactions) were run long enough, after an initial transient phase, to produce about 1500 "large" events (as defined below), a total model time of 200.000 years.The inter-event times for large events ranges from O to 1391 years, with an average of about 130 years.In reality the faults would undergo significant evolution over a period of 200,000 years, and this period is used only in the sense that it allows adequate sampling of all the interaction modes that might be pertinent now.We do not think that the stress state in the region, and the dates of recent large earthquakes, are anywhere near well enough known that we can initialise the model to the present conditions and expect to explicitly predict the near future.This is especially true since the synthetic catalogues are "chaotic" in that the details depend strongly on the initial conditions: models with only slightly different heginnings produce catalogues that over time become less and less similar.Only the long-term statistics are invariant.

RESULTS
Our results are presented in terms of the distribution of interevent times, dT, for "large" earthquakes.This requires that we define what constitutes a "large" earthquake.We have done this by picking a magnitude threshold, 7.2, above which an earthquake is "large".This cutoff was chosen on the basis of the frequency-magnitude distribution.taking it as that magnitude above which the number of earthquakes begins to increase over the number expected by extrapolation from lower magnitudes [9].Most earthquakes of magnitude 7.2 or more would produce an intensity of at least MM VIII in central Wellington City (using the attenuation relationship of [3]), t I 0 10 20 km 27 implying significant structural damage, and liquefaction at susceptible sites.On some faults earthquakes as small as magnitude 7 .2will not rupture the full extent of the fault plane: magnitudes of characteristic events arc about 8.1, 7.9, 7.5, 7.4, 7.3, 7.2, 7.6 on the subduction thrust, Wairarapa, Wellington South, Wellington North, Ohariu South, Ohariu North, and Wairau faults, respectively.

FIGURE 2. Surface traces of the major strike-slip faults in the Wellington region (modified from [12]).
The results (Figure 4) show a large increase in the number of short inter-event times (dT of 0-3 years) when fault interactions arc considered (45.7% of all dT's) as compared to the situation where interactions are ignored ( 1.8% ), a factor of 25.That is, the model predicts that nearly half of the large earthquakes, M 7.2 or more, occurring in the Wellington region will be followed within 3 years by

Surface projection of the region of seismic slip on the subduction thrust 175°E
another large earthquake.For intermediate inter-event times ( dT of 3-300 years) the relative numbers are in the opposite sense, as would be expected (because the total number of events is about the same).However, for long inter-event times (dT more than 300 years) the relative number for the situation with interactions is again higher.As discussed in our previous work [9] this pattern reflects the effect of both mutual enhancement and inhibition between faults.0 20 km

FIGURE 3. Faults included in the synthetic seismicity model. The surface projection of the seismogenic portion of the subduction thrust is marked by dashed lines.
These results can be expressed as ratios by dividing the numbers of c/T's for the model with interaction by Lhe corresponding numbers for the case without interactions (Figure 5).These ratios will be 1.0 for all dT's if interactions have no effect.If interactions lead to mostly inter-fault enhancement ( i.e., clustering of earthquakes) then the ratios would be > 1.0 for small c/T's and < l.O for large c/T's.If mutual inhibition predominated we would expect the opposite pattern.

DISCUSSION
Our results indicate that when fault interactions are considered there is a large enhancement in the likelihood of multiple large earthquakes within a short time ( < 3 years).as compared to the corresponding cases with no fault interaction.Although the real situation is more complex, if we take the simple view that, given a short inter-event Lime, the first event helped to "induce" the second, then we can examine the catalogue to see which faults tend to enhance each other.We find that large earthquakes on the subduction thrust are by far the mosl likely to induce other large earthquakes (47% of all short interevenl times), and that the "induced" earthquake is most often on the Wairarapa fault.
The opposite situation (earthquakes on the Wairarapa fault induce events on the subduction thrust) is also common.Given this predominant role of the subduction thrust, and that the amount of seismic slip on it is not well known, we have run another set of models identical to the first except that the subduction thrust is removed.We still find shortlerm (0-3 year) clustering of earthquakes, but the effect is much reduced: some 99; of all inter-event times fall into that range instead of Lhc previous 46%.This points out the importance of paleoscismic research into the past seismic behaviour of Lhe subduction Lhrust if tight constraints on the degree of clustering are desired.
Of course, if a large earthquake occurs we will know what fault was involved.and could refine the estimates of Lhe chances of another large event.The percentage of large earthquakes followed (dT< 3 years) by another large event are 86%, 68%, l 7C/c, 55%, 48%.10'7c, and 4<;,f, for the subduction thrust, Wairarapa.Wellington South, Wcllinglon Norlh, Ohariu South.Ohariu North, and Wairau faults respectively.
In our models (wilh or without the subduction Lhrust) a large earthquake on either of the segments of the Wellington fault has a 15% chance of "triggering" an earthquake on the other segment wilhin 3 years.Indeed, in the model without the subduction thrusl Lhis segment interaction is the main cause of the observed clustering (55% of all short inter-event times).Although this triggering rarely happens immediately (i.e., a big earthquake on one segment jumps Lhe gap and continues on the other segment with no delay) it is possible that some other sets of faull mechanical properlics, or dynamical effects we neglect, would make this more likely.This would reduce the clustering effect by coalescing Lwo moderately large earthquakes into one very large event.
Because we must limit the region we model, along strike triggering on olher faults is neglected, and clustering somewhat underestimated.The large earthquakes in 1848 (Awalere fault) and 1855 (Wairarapa fault) could he an example.

29
As mentioned above, the fact Lhal the number of very long inter-event limes is greater for the model with interactions than for the model without can be explained as the effect of mutual fault inhibition in some situations.A similar analysis to the above shows that a large event on one of the strikeslip faults Lends, on the average, to inhibit events of the other parallel faults.However, this is not always the case, depending on the faults involved and on details of the slip distribution.
We can understand why large subduction thrust earthquakes in our models might enhance the chance of a large earthquake on the overlying Wairarapa faull by examining the induced stress due to a "typical" characteristic earthquake on the subduction thrust.The induced stress can be expressed in terms of a change in the "Coulomb Failure Function" (OCFF) [6], defined in our models as oCFF = Otshear + µ(otn + oP) (3) where we take Otshear as positive if it favours slip in a dextral sense (as observed for large evenls on the Wairarapa fault).For a typical subduction thrust event Lhe average oCFF is +5.36 MPa, and is uniformly positive except for the extreme uppcr/southwesl corner of the fault.Since the partial stress drop following a past large Wairarapa evenl (aboul 20%) would have left lhe shear stress still fairly high (about 72 MPa), and the average strength is about 90 MPa, the average induced stress represents a "bringing Lhe clock forward" of about 30"/c (368 years) of the "characteristic recurrence" time (l 226 years).
In the results presented above we found that a time interval of O to 3 years after a large earthquake is the period in which there is an enhanced chance of another lan!e event.The "3" should not be taken as particularly well tonstrained: some different hut not unreasonable fault mechanical properties would lead to a different time period, from essentially O to about 10 years.But there is always some short time period, in comparison to the average inter-event time, in which enhancement occurs.
Another point we wish lo emphasise is that fault interactions destroy any simple periodicity of characteristic earthquakes on a given fault.For example, the mean recurrence time of events on the Wellington fault (south segment) is about the same in our model with or without interactions (553 and 490 years), but the ranges arc 86-1949 and 454-514 years respectively.This effect is mainly due to the "roughening" of the otherwise smooth stress buildup on a fault due to slips elsewhere.IL can also lead to the occurrence of doublets (two large events on different sections of the same fault within a short time).The situation is not so extreme if the subduction thrust is neglected hut it would be dangerous to extrapolate observations of past quasi-periodicity (as for the Wairarapa fault) in a region where interactions are likely to he strong.
There arc several aspects of our model that lead to some uncertainty as to the "exact" degree of clustering of large earthquakes to be expected in the Wellington region.First, as pointed out above, the position of the locked section of the subduction thrust is still not known precisely.We have used the position given in [2], but a more definitive result will emerge as recent geodetic results are analysed.In our previous work [9] a large number of possible positions were examined.almost all leading to some degree of clustering.Also.neither the percentage of the arc-normal convergence rate that is taken up by seismic slip on the interface, nor the amount of co-seismic slip on the interface arc well constrained.We hope that future paleoseismic studies will enable better estimates of these key parameters.A further complication could be time-variable coupling between the subducting and over-riding plates.as suggested by [14].Second, we have been forced by computational constraints to model only a limited region, and the effect of earthquakes on extensions of the major faults to the north-east and southwest is ignored.As time passes improved computer resources will help to address this problem.Third, some smaller faults that arc capable of generating earthquakes of about magnitude 7 .0-7.2 are not included in the model, again because of computational constraints.These faults include those involved in the two 1942 Wairarapa earthquakes.which are not well defined geologically.Fourth.the "driving mechanism" that loads the faults towards failure (described in [9]) is somewhat simplified.We eventually hope to improve this by using finite element modelling of the plate convergence process.

CONCLUSION
Although some aspects of our modelling remain uncertain we nevertheless consider that temporal clustering of large earthquakes in the Wellington region is very likely.It is not tenable to assume that if there is a large earthquake in the region then the "stress has been relieved" and that another similarly large event is unlikely in the near future.This certainly has implications for the earthquake engineering, insurance.and emergency management communities.For example, do the providers of earthquake insurance have the capacity to pay large claims twice in close succession~ Also, following a large earthquake there will be a number of badly damaged structures thal will still be standing and they will have low resistance to another strong event.While such structures may be regarded as unsafe due to aftershocks (for a month or two).the threat of another large earthquake within a few years is usually not taken into consideration.

FIGURE 4 .
FIGURE 4. Distribution of inter-event times, in years, for earthquakes of magnitude 7.2 or more.Top: with fault interactions; Centre: without fault interactions; Bottom: with interactions but the subduction thrust removed.

FIGURE 5 .
FIGURE 5. Ratios of the numbers of inter-event times (for events of magnitude 7.2 or more) for the case with fault interactions to the numbers for the case without interactions.Solid circles: model with seismic slip on the subduction thrust (preferred).Open circles: model without seismic slip on the plate interface.Note the logarithmic scales.

large earthquake with right-lateral slip on the vertical fault A (map view) would induce stress changes that would enhance the chances of a large event with similar slip and on similarly oriented faults in the shaded regions
FIGURE I.A ( e.g., fa ult B), and inhibit them in the other regions (e.g., most off ault C).