A MODEL FOR MM INTENSITIES NEAR LARGE EARTHQUAKES

The attenuation model for Modified Mercalli intensities that is currently in use in New Zealand (Dowrick & Rhoades, 1999) was developed from the available intensity data from large local earthquakes in New Zealand, but it does not represent well the intensity patterns that are expected when large earthquakes occur on long faults (length 20 km or more). This is because very few such events have occurred in New Zealand in historical times. An attempt to account for elongated source geometries has resulted in a model which provides a plausible extension to the Dowrick & Rhoades model. It also addresses detail in the intensity data from New Zealand's four largest historical earthquakes, that has not previously been accounted for. In development of the new model, stochastic terms have been added to represent the effects of asperities or areas of large slip on the rupture surface and to account for uncertainty in the fitting of the original data.


INTRODUCTION
have prepared a detailed mathematical model of the Modified Mercalli intensities likely to be experienced in earthquakes in New Zealand.Their work supersedes that done by Smith (1978Smith ( , 1990)), in particular because it uses a revised set of intensity data and moment magnitudes.The work is particularly important for hazard modelling, because it provides a tool for estimating the likely ground motion at any place, resulting from an earthquake elsewhere, given the parameters of its location, magnitude and focal mechanism, and because MM intensity tends to correlate well with damage.With further information about the vulnerability of assets, it allows risk assessment to be done, and this is useful for a number of applications in planning and insurance.
The Dowrick & Rhoades (1999) formula (hereinafter referred to as the DR formula) specifies coefficients for an intensity function of the form (1) where I is intensity, Mw is the moment magnitude, r a measure of distance, h a measure of depth, and A1 to A4 are coefficients.The integer MM intensity is derived from the value of I by truncating the decimal, not by rounding.The formula provides for the variety of earthquake mechanisms and tectonic regimes in its definition of coefficients A I to A4.Further, it defines elliptical isoseismals by specifying intensity as a function of distance in two normal directions: those of the major and minor axes.So in order to use the formula, it is necessary to know not only the location and moment magnitude of the earthquake but also the mechanism (reverse, normal, strike-slip), the tectonic regime (crustal, volcanic, slab event, interface, etc) and the orientation of the major axis of the elliptical isoseismals.For an earthquake on a known fault, this last can be assumed to be the strike of the fault.If the orientation is not known, the formula has coefficients that can be used to define circular isoseismals.
Where a source is an extended one, such as on a fault where the surface trace is well defined, it is common in modelling other measures of strong ground motion, such as peak ground acceleration and spectral measures, to account for the geometry by using the closest distance to the rupture surface as the distance measure in the attenuation formula.This may be the perpendicular distance to the surface trace or, for dipping faults, it may be to a point at depth.Dowrick and Rhoades have not done this, and it is useful to consider why that may not be appropriate for intensity modelling.At a site near a rupturing fault, the particular pulse that has the highest peak ground acceleration is likely to have come from the nearest part of the rupture surface.So the shortest distance to the rupture surface is a plausible distance measure to use in an attenuation formula.But intensity is inherently different from peak ground acceleration.It probably represents an integral of the entire ground motion history, certainly incorporating the duration of the signal.So each location where the earthquake is experienced must be considered to be irradiated by the entire rupture surface.
In using an elliptical formula, Dowrick & Rhoades (1999) addressed the fact that isoseismals are likely to be elongated near the source of a large earthquake, but in fact they had very few clear examples of large earthquakes with long surface ruptures in their database, so their formula cannot be expected to model such source geometry well.
This study takes the approach that, in the absence of adequate data from large earthquakes with long source geometries, an acceptable model for estimating ground motion near such sources must instead draw on geophysical understanding of the earthquake source.It therefore proposes a plausible intensity model, which is derived from the DR model but modifies it to take into account the characteristics of an elongated source and the few available data from large earthquakes.It also incorporates a stochastic component in order to address other characteristics of large sources.

THE NEED FOR AN ADJUSTMENT
The DR formula is a point-source formula.It is defined in that way, in terms of the basis functions that comprise it.When used to compute ground motion from large events, it is to some extent an extrapolation of data from smaller events, so it cannot be expected to represent elongated sources well, at short distances.This is for a very good reason: there are very few examples of large earthquakes in New Zealand that were accompanied by long ground ruptures.Stirling et al. (2000) have developed a seismicity model for New Zealand, which they have used for Probabilistic Seismic Hazard Assessment.This assessment is in terms of peak ground acceleration (PGA) and spectral acceleration, and there is a need to perform such an assessment in terms of MM intensity.The seismicity model identifies 305 active faults and a grid of sites at which parameters are defined to represent distributed seismicity.For each fault a likely mechanism, characteristic magnitude and mean return period are estimated.So it will be necessary to model intensities from earthquakes on these faults.Figure 1 shows the Wellington-Hutt Valley segment of the Wellington fault, which Stirling et al (2000) represents as four linear sections.
It is a vertical fault, with strike-slip movement, and is considered to have a characteristic magnitude of 7.3 (Stirling et al, 2000).Also shown in Figure 1 are the isoseismals as calculated by the DR formula.A major axis for the isoseismals has been defined by fitting a straight line through the mid-points of all four sections.
Figure 1 illustrates the main issue this study addresses.The DR formula predicts that the intensity drops to MM VIII at the ends of the fault.As will be shown in Figure 3 below, available data from large New Zealand earthquakes suggest that MM IX is more likely at the ends of the fault rupture.
The issue is even clearer when a larger earthquake is considered.Figure 2 shows a representation of the Alpine fault, for which the DR formula predicts an intensity of MM VI at the ends of the rupture, and this seems too low for an earthquake of magnitude 8.07.
It is clear that an adjustment will only be necessary for the largest earthquakes (possibly 7 and greater).It is also clear that there are few data to constrain the adjustment.But if the goal is a realistic model of the ground motion that is expected in large events, an adjustment must be made.

97
At the same time, another characteristic of large earthquakes can be addressed.In very large earthquakes the highest intensities (i.e.MM X and maybe greater) are not generally observed to occur along the whole length of the rupture (e.g.Takemura & Ohno, 1996).Rather, they tend to occur at isolated locations, and this is understood in terms of the existence of asperities on the rupture surface or of variation of the amount of slip along its length.So in extending the high intensities along the fault trace to take account of its length, it may be possible to retain the localised occurrence of the highest intensities.This is addressed in Section 3.3 below.

DETAILS OF THE NEAR-FAULT INTENSITY MODEL
The proposed model has four aspects.
• Within an elliptical region for which the major axis is the fault trace, intensities are assumed to be equal everywhere except in a localised region where they are higher.

•
Outside that ellipse, the formula merges continuously into the DR formula at a distance of two fault lengths.

•
The localised high intensity region can occur at a random location along the fault trace, and intensities there are given by the DR formula with distances measured from the local centre.

•
No change has been made to the DR formula in the minor axis direction, i.e. normal to the fault trace.

A plateau ellipse of high intensity at the fault
It is proposed that the elliptical isoseismals of the DR formula be extended, in a direction parallel to the fault strike, in order to accommodate the fault geometry.
Figure 3 shows the isoseismal data as used by Dowrick & Rhoades (1999) for the four largest New Zealand earthquakes in historical times.These data are the distances from the centre of the strong motion pattern to the isoseismals, in the direction of the major axis.Details are given in Table 1.In each figure the solid line gives the intensity predicted by the DR formula.Note that this is in the direction of the major axis for the ellipses defined in the DR formula.For each earthquake, the observed isoseismal radius for intensity MM IX exceeds that given by the DR formula.For the 1929 Buller earthquake (Fig. 3b) the MM X isoseismal also exceeds the DR value.It is therefore proposed that intensity MM IX be regarded as the intensity that accompanies surface faulting, i.e. that whatever the magnitude, intensity MM IX is expected to occur out to the ends of the fault trace.
Intensity is of course defined as an integer parameter.It is for this reason that Roman numerals are commonly used.But formulae for modelling intensities necessarily produce decimal numbers.Dowrick and Rhoades (1999) followed the convention suggested by Smith (1978), that the conversion from decimal values to integers be by truncation of the decimal, rather than rounding.This is so that the calculated intensity takes integer values at the isoseismals.So a calculated intensity of 7.6, for instance, occurs somewhere between the MM VII and MM VIII isoseismals, where the intensity is regarded as MM VII.In the sense of a quantity that can be calculated, it is helpful to consider intensity as a continuous parameter, but for which experimental measurements are approximate.
What decimal value should be chosen for the generic "faultbreak" intensity?The data in Figure 3 suggest MM IX, which implies 9.0 as a lower limit.9.5 is a reasonable upper limit, as mid-range MM IX.I have selected 9.2 as a realistic yet conservative choice.The broken lines in Figure 3 show the current proposal for an intensity model, where it exceeds that of the DR formula: a plateau at intensity 9.2 which extends as far as the end of the fault trace.Figure 3 shows fault lengths as given in Table 1 (see sources listed there).The model merges with the DR model at some greater distance.
Figure 3 has a logarithmic distance scale, so the differences between the Dowrick-Rhoades formula and the proposed one do not appear to be large.Table 2 compares the half-lengths of these four faults with the distances at which the DR formula predicts intensity 9.2.
175° The procedure for modifying the DR formula can be summarised as follows: a.
Does the DR formula predict an intensity greater than 9.2 at the epicentre, and less than 9.2 at the end of the fault?Only if both these conditions hold will any changes to the DR formula be made.
b.An elliptical plateau of intensity 9.2 will extend to the end of the fault trace (major axis), with minor axis that given by the DR formula for intensity 9.2.See Figure 4, where x = p at the end of the fault trace.c.Determine at what distance along the fault trace the DR formula predicts that intensity 9.2 occurs (x = a in Figure 4).

An algorithm for merging with the DR formula at greater distances.
In Fig 5,point P (x,y) is outside the plateau ellipse.Smith (1995) proposed an ad hoc procedure for modelling intensities near long fault sources, and found that the geometry of the source seems to have no effect on intensities at distances beyond about two fault lengths.Using this result, the intensity function for point P is defined as follows.
a. Determine x and y, the coordinates in directions parallel to and normal to the fault strike, respectively.b.Replace x by x -XcYc, where p and a are as in Figure 4 and if y > 4p Ye= 0 c.Compute the DR intensity at point P'.This procedure forces continuity with the DR formula at x or y equal to two fault lengths.The new formula is not smooth at these distances, but continuity is assured.This replacement of x by x -XcYc maps point P to P', closer to the origin.For negative values of x, symmetry about x = 0 is assumed.The elliptical isoseismalfor intensity 9.2 is expanded in the along-strike direction so that its major axis is the same length as the fault rupture.

The highest intensities are more concentrated
Intensities of about MM X and greater are not expected to be experienced along the whole length of the fault trace.Rather they are localised because the amount of slip varies along the rupture surface or because of the presence of asperities which generate stronger motions.Takemura & Ohno ( 1996) have contoured the strong motion field from the Hyogo-ken Nanbu (Kobe) earthquake, using a source model with maximum slip near one end of the fault source.But for the purposes of predicting ground motion in large earthquakes, it is not possible to determine where the high intensity region will be.
A prudent measure, therefore, is to allow it to be at a random location, somewhere along the fault trace.
For intensities higher than 9.2 the isoseismal pattern is therefore modelled as being centred not at the mid-point of the fault but at a random location somewhere along it, and the intensity is determined by the DR formula with its origin at that random location.A schematic is shown in For scenario studies, the high intensity region can be placed wherever desired.For stochastic studies it should be given a random location somewhere along the rupture surface.But in either case this location of the peak intensity may be displaced from the mid-point of the fault only as far as the point from which the DR intensity drops to 9.2 by the end of the fault.

UNCERTAINTY
Two sources of uncertainty have been built into the model.The first is the location of the high intensity region, as described in Section 3.3.This affects sites that are within the plateau ellipse (Fig 6).The modelled intensity at any such site may be 9.2, or it may be a higher value if the randomly assigned high intensity region happens to be close enough to the site.
A second source of uncertainly relates to the fact that the DR formula was derived by fitting intensity data.Dowrick & Rhoades (1999) evaluated uncertainties in three ways.Their term , is the between-earthquake error term, which relates to the degree of goodness of fit of the isoseismal radii., takes the value 0.26 intensity units for the major axis direction, 0.21 for the minor axis direction, So to account for this uncertainty I add to the computed intensity a normally distributed random error term with zero mean and standard deviation equal to the mean of these two values of ,, i.e. 0.235.Because this represents the between-event variation, the same error term should be added to determinations of intensity at all sites, for any one earthquake.Dowrick & Rhoades (1999)  These error terms mean that repeated computations of intensity at a particular site, from a given source, will not be identical.Instead they will be scattered, and this variation should be built into any hazard estimates that are made.

MORE COMPLICATED GEOMETRIES
Long faults are rarely mapped as straight lines.The Wellington fault, for instance, is modelled as four linear sections in Figure 1 but is actually much more complicated.A methodology is needed to represent such sources, for the purpose of calculating intensity.
The technique used in Figure 1 is proposed.That is, a straight line is fitted by linear regression to the mid-points of the linear sections.The projection of all the sections on to this straight line is used for the major axis of the plateau ellipse, with intensity 9.2 as described as in Section 3.1.
For dipping faults, Dowrick & Rhoades (1999) described a geometry in terms of the upper edge of the rupture surface.This seems to introduce lateral bias, in that intensities will be higher on the down-dip side of this line.Instead, it is proposed that the major axis of the isoseismal pattern be located over the centroid of the rupture surface, as shown in Figure 7.Some bias may still be present, due to the fact that surface locations on the down-dip side of this axis are further from the rupture surface than locations up-dip.But it is better than using the top of the rupture surface as an axis.
An attempt has been made to account for complicated geometries by a procedure such as that proposed by Smith (1995), in which an integration is performed over the rupture surface.These efforts have not proved successful, and the simpler approach of the preceding two paragraphs is preferred.
Figures 8 and 9 show the isoseismals for the fault sources in Figures l and 2, as predicted by the new formula.See also Table 2.

CONCLUSIONS
The DR formula for calculating likely intensities is limited in its application to earthquake sources associated with long fault ruptures, or indeed any extended source geometries.The new development provides a plausible means of accounting for elongated sources.It does explain characteristics of the isoseismals of the four largest shallow New Zealand earthquakes that are not predicted as well by the DR model.The development is also based on judgement about likely intensities near fault ruptures, and on considerations of the effects of rupture surface asperities on ground motion.
Modification of the DR formula paves the way for development of a probabilistic seismic hazard model that is expressed in terms of MM intensities, to complement the model that has already been developed (Stirling et al, 2000) for peak acceleration and spectral acceleration.
It is proposed that this new intensity model be used instead of the Dowrick-Rhoades model, unless and until a more rigorously obtained model for near-fault intensities is available.
Major axis DR axis .

. . . centroid
View along strike of a dipping fault.The DR fonnula detennines distances from the top surface of the rupture, but it is proposed to measure from the centroid so that the major axis of the isoseismal pattern is located above this point.

d.
For each point on this DR ellipse, determine x and y coordinates as parallel and normal to the fault trace, respectively.Scale the x-coordinate by pla to give the x-coordinate for the modified DR ellipse (y coordinate unchanged).
& Webb, In prep.)Data for events in Figure 7. Mw = moment magnitude, Mech= mechanism (reverse or strike-slip}, he= depth to centroid of rupture surface, ht = depth to top of rupture surface, L = source length.All data are as used by Dowrick & Rhoades (1999) except L, for which references are given.* Dowrick & Rhoades used a reverse mechanism for the 1931 event, but have since (pers.comm.)preferred a strike-slip mechanism.

FigFigure 5 .
Figure 5.Intensity at point P(x,y) is determined according to the DR formula as if it were at P'. Mapping from P to P' is described by equation 5.
Figure 9.Proposed isoseismal modelling for an Alpine fault earthquake (cf Fig.2) This particular random location of the high intensity region is part-way toward to southern end of the fault.

of fault half-length and the predictions of the DR formula for the four major earthquakes and for the Wellington and Alpine faults.
Dowrick & Rhoades (1999) error term a, and it is important to recognize what this refers to.It represents the goodness of fit of the individual isoseismal radii for an event.But the present exercise is to provide a tool for estimating spot intensities rather than (smoothed) isoseismal radii, so this term has not been included.Rather,I note the finding byDowrick & Rhoades (1999)that within a single isoseismal zone, the spot intensities are distributed with mean -0.18 and standard deviation 0.64.So at individual sites a further correction is added: a normal variate with these parameters.