Peter J Stafford

John B Berrill
email: john.berrill@canterbury.ac.nz

Jarg R Pettinga
email: jarg.pettinga@canterbury.ac.nz

New predictive equations for Arias intensity from crustal earthquakes in New Zealand

Keywords: Arias intensity, New Zealand, Crustal, Earthquakes, Ground-motion prediction, Attenuation 0.1

Arias Intensity (Arias, MIT Press, Cambridge MA

Introduction

New Zealand is commonly recognised as a country with significant levels of seismic activity by worldwide standards, yet very few empirical ground-motion models have been derived for the country. This study takes a step towards remedying this situation by presenting new predictive equations for Arias Intensity from crustal earthquakes in New Zealand over a wide range of magnitude and distance values. While the utility of Arias Intensity is now accepted, it has not yet been incorporated directly into many engineering applications; for example, vulnerability functions for building damage are very seldom expressed in terms of Arias Intensity. As the specific applications of Arias Intensity have not yet been defined for many circumstances, this study presents a suite of equations that users may select from in order to best meet their needs. The provision of a suite of equations also allows for the epistemic uncertainty associated with the prediction of Arias Intensity to be accounted for in probabilistic seismic hazard analyses.

Arias Intensity [START_REF] Arias | A measure of earthquake intensity[END_REF]) is an important measure of the strength of a ground motion, as it is able to simultaneously reflect multiple characteristics of the motion in question. Whereas most common scalar ground-motion measures, such as peak ground acceleration, peak ground velocity, or individual ordinates of spectral acceleration, reflect a very specific aspect of the ground motion, Arias Intensity, whilst remaining a scalar measure, is able to capture and represent multiple attributes of the overall ground motion. In particular, the ability to reflect the energy of a groundmotion signal distributed across a broad frequency band or, equivalently, to reflect the influence of the entire duration of a ground motion is a property of Arias Intensity that lends itself to prudent application in many structural and geotechnical applications. [START_REF] Travasarou | Empirical attenuation relationship for Arias intensity[END_REF] have demonstrated the effectiveness of using Arias Intensity as a predictor of the likelihood of damage to short-period structures, showing that the correlation between Arias Intensity and structural damage is, in this case, stronger than that between peak ground acceleration and structural damage. Arias Intensity is also a very useful ground-motion measure that can be used in geotechnical applications such as in the determination of the likelihood of rock falls and landslides [START_REF] Harp | Shaking intensity thresholds for rock falls and slides: evidence from 1987 Whittier Narrows and Superstition Hills earthquake strongmotion records[END_REF][START_REF] Gaudio | An approach to time-probabilistic evaluation of seismically induced landslide hazard[END_REF] or for estimating the susceptibility of a site to liquefaction [START_REF] Egan | Assessment of earthquakeinduced liquefaction using ground-motion energy characteristics[END_REF][START_REF] Kayen | Assessment of liquefaction potential during earthquakes by Arias intensity[END_REF][START_REF] Kramer | Ground motion intensity measures for liquefaction hazard evaluation[END_REF]. In applications such as these, it is not necessarily the peaks of a strong-motion record that most strongly govern the response. Rather, it is a combination of the frequency content of the motion, the duration of, or number of cycles in, the motion, and the amplitude of the motion. The Arias Intensity is able to capture all of these characteristics for any given ground motion.

The recognition of this utility has led to probabilistic seismic hazard analysis (PSHA) being conducted in terms of Arias Intensity in addition to the analyses carried out for more traditional ground-motion measures. An example is the recent PSHA for southeastern Spain by [START_REF] Peláez | A preliminary probabilistic seismic hazard assessment in terms of Arias intensity in southeastern Spain[END_REF]. In order to conduct hazard analyses in terms of Arias Intensity, stable empirical ground-motion models must be available for this purpose. However, in comparison with other measures of strong ground motion, there are very few models for Arias Intensity that could be used within a PSHA [START_REF] Wilson | (ed) Evaluating earthquake hazards in the Los Angeles Region -an earth-science perspective[END_REF][START_REF] Keefer | Predicting earthquakeinduced landslides, with emphasis on arid and semiarid environments[END_REF][START_REF] Wilson | Relation of Arias intensity to magnitude and distance in California[END_REF][START_REF] Sabetta | Estimation of response spectra and simulation of nonstationary earthquake ground motions[END_REF][START_REF] Travasarou | Empirical attenuation relationship for Arias intensity[END_REF]. Of the few models that are available, most must now be considered obsolete, as they do not make use of the large number of high-quality strong-motion records that have been compiled over the past decade or so. Some of the most recent work considering Arias Intensity has focussed upon earthquake-specific scaling that cannot be used within a PSHA framework [START_REF] Hwang | Attenuation relations of Arias intensity based on the Chi-Chi Taiwan earthquake data[END_REF]). Of the potentially usable models, the model of [START_REF] Travasarou | Empirical attenuation relationship for Arias intensity[END_REF] is based upon, by far, the largest dataset and uses the most comprehensive functional form for the regression model. A big weakness of the remaining models is their very limited consideration of site response; most are either derived for a specific site class or use dummy variables to distinguish between site classes. In addition, none of the models (with the exception of [START_REF] Travasarou | Empirical attenuation relationship for Arias intensity[END_REF] consider differences between ground motions generated by earthquakes of different mechanism.

When deriving predictive equations for strong ground-motion measures, an analyst must always confront the issue of which functional form to use for the model. Usually, the final model that is presented is arrived upon following a trial and error procedure; trying different functional forms, assessing the significance of the regression coefficients thus obtained, checking for strong correlations amongst the coefficients, considering the standard deviation of the prediction, etc. However, it is conceivable that the analyst may discover multiple functional forms that are able to model the data with a similar degree of efficacy. One specific model will always score better than others following the application of some statistical test such as the Akaike Information Criterion (AIC; [START_REF] Akaike | A new look at the statistical model identification[END_REF] or the Bayesian Information Criterion (BIC; [START_REF] Schwarz | Estimating the dimension of a model[END_REF]), but one must recognise that this score judges the ability of the model to capture trends in the particular dataset that is under consideration and does not preclude the possibility of better performance from alternative models when modelling ground motions from future earthquake events.

In modern PSHA, multiple ground-motion models are almost always considered in order to capture a degree of the epistemic uncertainty associated with predicting ground motions given some limited set of predictor variables, such as magnitude, distance and site class. Each of these models is usually taken from a different study drawing upon different strong-motion datasets that include recordings from different earthquakes. However, given that the derivation of predictive equations necessarily requires the consideration of various functional forms, it makes sense to present a suite of predictive equations having different functional forms derived from a single strong-motion dataset. Therefore, in the present study, rather than presenting a single ground-motion model, a suite of models is presented so that these may be used to capture epistemic uncertainty within a PSHA framework or, alternatively, so that the end user may select a model that they believe to be the most appropriate rather than having this decision made on their behalf.

An additional problem that is associated with utilising suites of predictive models to capture epistemic uncertainty in PSHA is that predictive equations are sometimes derived for various different horizontal-component definitions (e.g. [START_REF] Beyer | Relationships between median values and between aleatory variabilities for different definitions of the horizontal component of motion[END_REF]. Commonly, some average of two horizontal components is used in the derivation of predictive equations, as the resulting standard deviation is smaller and the regression analysis more stable. However, there are situations in which hazard estimates for single component definitions may be more appropriate [START_REF] Baker | Which spectral acceleration are you using?[END_REF]. For this reason, as well as obtaining coefficients for a suite of functional forms, coefficients are also derived for various different horizontal-component definitions for each model.

It is shown that all of the models that constitute the suite of equations that are presented have functional forms that are broadly consistent with theoretical considerations and that the differences between them reflect commonly held beliefs regarding how to most appropriately model strong ground motions.

Strong ground-motion dataset

Although New Zealand straddles a major plate boundary, there are relatively few strong ground-motion records when compared to other developed countries in similar tectonic environments such as Japan, Taiwan and western North America. The New Zealand dataset, which is openly available via the GeoNet project (http:// www.geonet.org.nz/) consists of recordings from events occurring in three distinctly different tectonic regimes: crustal, interface and subduction, as well as including recordings of waves that have passed through the central volcanic plateau where wave amplitudes are observed to attenuate at a heightened rate [START_REF] Zhao | Attenuation of peak ground accelerations in New Zealand earthquakes[END_REF][START_REF] Cousins | A model for the attenuation of peak ground acceleration in New Zealand earthquakes based on seismograph and accelerograph data[END_REF][START_REF] Eberhart-Phillips | Estimating slab earthquake response spectra from a 3D Q model[END_REF][START_REF] Mcverry | New Zealand acceleration response spectrum attenuation relations for crustal and subduction zone earthquakes[END_REF]. A rigorous development of predictive equations for strong ground-motion modelling in New Zealand would include treatment of each tectonic regime and make allowance for scenarios involving wave propagation through the central volcanic plateau (i.e. [START_REF] Mcverry | New Zealand acceleration response spectrum attenuation relations for crustal and subduction zone earthquakes[END_REF]. However, in the present study, the Arias Intensity is modelled for crustal earthquakes only, and no account is made for volcanic path effects.

The relatively low number of suitable records that are available for the regression analyses prompted us to follow the lead of [START_REF] Zhao | Attenuation of peak ground accelerations in New Zealand earthquakes[END_REF], [START_REF] Cousins | A model for the attenuation of peak ground acceleration in New Zealand earthquakes based on seismograph and accelerograph data[END_REF] and [START_REF] Mcverry | New Zealand acceleration response spectrum attenuation relations for crustal and subduction zone earthquakes[END_REF] and supplement the New Zealand dataset with foreign strong-motion records. The selection of suitable foreign records has been greatly facilitated by the dataset for the Next Generation Attenuation (NGA) models being made available via the internet (http://peer.berkeley.edu/nga/). This large, well-constrained dataset was restricted to exclude recordings from the 1999 Chi-Chi earthquake and related aftershocks. Thereafter, only free-field records for which site classification according to NZS1170.5:2004(Standards New Zealand 2004) could be made were included. In most cases, classifications made on the basis of shear-wave velocity are possible.

The New Zealand earthquakes that have been considered as part of this study are listed in Table 1, whilst the foreign records may be obtained from the NGA_Flatfile (Version 6, 22/4/2005). The composite ground-motion dataset is shown in terms of the magnitude-distance distribution in Fig. 1, in which contributions from the NZ and foreign datasets are identified as are contributions from different fault mechanisms and site classes. From inspection of Table 1 and Fig. 1, it can be seen that the dataset includes records from earthquakes having moment magnitudes ranging between 5.08 and 7.51 (of which the New Zealand data spans the range 5.08 to 7.20) and includes ground motions recorded at distances of up to 300 km (with the foreign data being limited to 20 km in Fig. 1). These limits on the predictor variables directly limit the strict range of applicability of the empirical ground-motion models that are derived from this dataset [START_REF] Bommer | The influence of magnitude range on empirical groundmotion prediction[END_REF]). In Table 1, the number of records that each event contributes to the total dataset is given, while in Table 2, the numbers of records coming from both NZ and foreign earthquakes may be seen. In this latter table, different numbers of foreign records are considered, and this feature will receive further attention later in the manuscript. From inspection of Table 1, it is immediately obvious that the number of records coming from each earthquake is relatively low and that there are considerably more reverse faulting events than there are strike-slip, normal or oblique events. However, although the New Zealand records are relatively scarce, those that have been included in the analysis are typically of good quality, and a considerable amount of attention has been paid to ensuring that the meta-data related to the records is as accurate as possible.

An extensive literature review was conducted in order to ascertain the most appropriate faultrupture models corresponding to the given earthquakes. For the majority of events in the dataset, focal mechanism solutions exist, and for some, 2003), ( 25) [START_REF] Mcginty | The 2003, Mw 7.2 Fiordland earthquake, and its near-source aftershock strong motion data[END_REF] additional constraint is available in the form of spatial aftershock patterns, geodetic modelling, elastic dislocation modelling, Coulomb stress change modelling and considerations of structural geology. Information from the many references cited in Table 3 was extracted in order to determine the finite fault rupture parameters for each event. In the worst-case scenario, the only information available for a particular event was the Harvard centroid moment tensor solution and an estimate of the earthquake magnitude and hypocentral position. In such cases, the fault rupture model was determined using the strike and dip of the focal mechanism solution coupled with the source-scaling relationships developed for New Zealand by [START_REF] Dowrick | Relations between earthquake magnitude and fault rupture dimensions: how regionally variable are they?[END_REF].

Selection of functional forms

Arias Intensity [START_REF] Arias | A measure of earthquake intensity[END_REF]) describes the cumulative energy per unit weight absorbed by an infinite set of single-degree-of-freedom (SDOF) oscillators having fundamental frequencies uniformly distributed in (0,∞) [START_REF] Travasarou | Empirical attenuation relationship for Arias intensity[END_REF]). In the most general case, the expression for Arias Intensity may be written as [START_REF] Kayen | Assessment of liquefaction potential during earthquakes by Arias intensity[END_REF]:

I xx (ξ ) = arccos (ξ ) g 1 -ξ 2 ∞ 0 a 2 x (t) dt, ( 1
)
where the term I xx (ξ ) represents the Arias Intensity experienced by SDOF oscillators with a damping ratio of ξ aligned in the x-direction, responding to ground shaking in the x-direction, defined by the acceleration time history a x (t).

Naturally, there will be similar expressions representing the response of SDOF oscillators aligned in orthogonal directions responding to orthogonal components of the same ground motion. In this sense, the expression given in Eq. 1 is just one term of a second-order tensor. The trace of this tensor I xx + I yy + I zz is an invariant, and consequently, all pairs of orthogonal axes passing through a predefined origin have the same Arias Intensity (i.e. I xx + I yy = const; [START_REF] Travasarou | Empirical attenuation relationship for Arias intensity[END_REF]. This point is worth emphasizing, as it implies that the use of the arithmetic mean of the two recorded horizontal components results in an orientation-independent measure of the Arias Intensity [START_REF] Boore | Orientation-independent measures of ground motion[END_REF].

The most common representation of the Arias Intensity is recovered for the case of zero damping in the SDOF oscillators. Given this condition, the expression in Eq. 1 reduces to that below in Eq. 2.

I xx = π 2g ∞ 0 a 2 x (t) dt (2)
The expression given in Eq. 2 portrays the Arias Intensity in terms of the integral of what may be considered as a transient stochastic signal. However, such a representation offers very limited insight into potential functional forms for the regression analysis. Therefore, in order to develop suitable functional forms for the predictive equations, we make use of Parseval's theorem, which states that the total power of a signal in both the time and frequency domains is equivalent. This theorem is restated for completeness in terms of natural frequency below.

∞ -∞ |a (t)| 2 dt = ∞ -∞ |A ( f )| 2 df (3)
Here, a(t) and A( f ) are an acceleration time history and the corresponding Fourier amplitude spectrum (FAS) of acceleration, respectively. Note that, as we are considering real-valued acceleration signals in the positive time-domain, the expression in Eq. 3 can be simplified to consider only positive times and frequencies without any loss of generality.

Following [START_REF] Boore | Simulation of ground motion using the stochastic method[END_REF], an expression for the FAS of acceleration may be given as in Eq. 4, in which C is a coefficient representing various parameters related to the source, M 0 is the seismic moment, f c is the corner frequency of the source spectrum, R is a distance measure, Q( f ) represents an anelastic attenuation function, β is the shear wave velocity and S( f ) represents a site transfer function.

A(f ) = CM 0 (2π f ) 2 1 + ( f/ f c ) 2 1 R exp - π f R Q( f )β S( f ), (4)
Equation 4 can then be incorporated into the description of Arias Intensity in order to obtain an expression containing familiar seismological descriptors rather than the transient acceleration time-series representation given in Eq. 2. Such a procedure has previously been applied by [START_REF] Olafsson | A theoretical attenuation model for earthquake-induced ground motion[END_REF] (for the derivation of a relationship for the closely related RMS acceleration), [START_REF] Travasarou | Empirical attenuation relationship for Arias intensity[END_REF] and [START_REF] Stafford | New empirical predictive equations for the Fourier amplitude spectrum of acceleration and Arias intensity in New Zealand[END_REF]. The purpose of introducing the above expression into that for the Arias Intensity is not to obtain an exact closed theoretical expression for the Arias Intensity, but rather to guide us in the selection of an appropriate functional form for the predictive models. With this in mind, it is prudent to make some simplifying assumptions to the expression in Eq. 4 in order to relax some of the frequency dependence inherent in this model. Two such simplifications shall be made. Firstly, the assumption of a frequencydependent anelastic attenuation function, Q( f ), may be relaxed to the case of a constant attenuation rate, Q 0 . This relaxation is not regarded as being overly inhibitive to the accuracy of the resulting model, as the true frequency-dependence of anelastic attenuation is not accurately known over a broad range of frequencies and is known to be regionally variable [START_REF] Abercrombie | A summary of attenuation measurements from borehole recordings of earthquakes: the 10 Hz transition problem[END_REF]. Additionally, good models have been obtained for the FAS using the assumption of constant anelastic attenuation in the past [START_REF] Anderson | A model for the shape of the Fourier amplitude spectrum of acceleration at high-frequencies[END_REF][START_REF] Olafsson | Estimation of earthquake-induced response[END_REF]). The second relaxation is made with regard to the site response term. While it is well known that site response varies with frequency and approximate models of this variation exist [START_REF] Boore | Site amplifications for generic rock sites[END_REF], it is also appreciated that our understanding of site response is limited (e.g. [START_REF] Boore | Can site response be predicted?[END_REF]). Consequently, for the purpose of incorporating the effect of site response into a predictive equation for Arias Intensity, we ignore any frequency dependence of this function and simply introduce constant site class factors instead. This is in keeping with almost all existing predictive equations for frequency dependent strong-motion indices [START_REF] Douglas | Earthquake ground motion estimation using strong-motion records: a review of equations for the estimation of peak ground acceleration and response spectral ordinates[END_REF].

One further adjustment to Eq. 4 can be made to make the derivation slightly more general. In Eq. 4, the geometric spreading is assumed to be perfectly elastic spherical spreading at a rate of R -1 . However, in reality, this theoretical spreading rate is seldom observed. Frequently, the rate of spreading is observed to be hyper-spherical, but there are also many examples of attenuation at rates lower than this. Therefore, in the following derivation, the geometric spreading rate is assumed to be R -ζ .

The form of the model for the FAS that we will carry into the development of the functional form of the Arias Intensity relationship can therefore be given by the expression below, where now the coefficient C has been expanded to reveal the crustal density, ρ, a reference distance, R 0 , the radiation pattern, R φθ , a partitioning factor, V, and a free surface amplification, F.

A( f ) = π R φθ V F M 0 f 2 ρβ 3 R 0 1 + ( f/ f c ) 2 1 R ζ exp π f R Q 0 β S i (5)
In Eq. 5, the term S i represents the site response with the subscript, i, corresponding to the various site classes incorporated into the final model. Now, inserting Eq. 5 into the expression for Arias Intensity in terms of frequency, we obtain the following equation.

I a = π 2g π R φθ VFM 0 ρβ 3 R 0 R ζ 2 S 2 i ∞ 0 f 4 1 + ( f/ f c ) 2 2 × exp - 2π f R Q 0 β df (6)
In order to evaluate the integral in Eq. 6, we follow [START_REF] Olafsson | A theoretical attenuation model for earthquake-induced ground motion[END_REF] and make a change of variable. Letting f = f / f c , κ Q = 2π R / Q 0 β, and then further letting λ = κ Q f c , the above integral expression can be recast as follows.

= f 5 c ∞ 0 f 4 1 + f 2 2 exp -λ f d f = f 5 c λ (7)
In this equation, the capital psi term is a function of lambda and involves forms of both the sine and cosine integrals, Si(λ) and Ci(λ), respectively.

= 1 - λ 2 {Ci (λ) [λ cos (λ + 3 sin (λ))] +Si (λ) [λ sin (λ -3 cos (λ))]} (8)
Unfortunately, upon first inspection, the indefinite form of the above integral expression precludes the selection of a simple functional form for a regression model based upon theoretical considerations. Fortunately however, as highlighted by [START_REF] Olafsson | A theoretical attenuation model for earthquake-induced ground motion[END_REF], the solution to the integral expression given in Eq. 7 can be very well approximated over a wide range of practical values by a surprisingly simple expression, i.e. by simply letting the expression for equal one. The deviation from unity is small over the range of lambda values corresponding to common engineering scenarios, with the worst case deviation being less than 20%.

It can also be noted that while the anelastic attenuation factor, Q 0 , has been assumed to be independent of frequency for this derivation; it is typically found to be a linearly increasing function of distance [START_REF] Anderson | A model for the shape of the Fourier amplitude spectrum of acceleration at high-frequencies[END_REF][START_REF] Olafsson | Estimation of earthquake-induced response[END_REF]. This is due to the fact that waves travelling greater distances typically penetrate deeper into more competent regions of the lithosphere where they experience relatively less attenuation than their counterparts in the shallower crust. It can therefore be appreciated that rather than being an increasing function of distance, the term κ Q is approximately constant with respect to distance.

Therefore, the theoretical model for Arias Intensity can be simplified to that in Eq. 9.

I a = π 3 R 2 φθ f 4 c M 2 0 S 2 i gρ 2 β 6 R 0 R 2ζ κ Q (9)
In the above equation, the terms for both the corner frequency and the seismic moment are essentially functions of moment magnitude [START_REF] Brune | Tectonic stress and spectra of seismic shear waves from earthquakes[END_REF][START_REF] Hanks | Moment magnitude scale[END_REF]. Making the relevant substitution and consequent rearrangements results in Eq. 10 in which the new terms a and d are constants from the M W -f c and M W -M 0 relationships respectively [START_REF] Stafford | Engineering seismological studies and seismic design criteria for the Buller region[END_REF].

ln [I a ] = ln π 3 R 2 φθ gρ 2 β 6 R 0 ln 10 (4a + 2d) +ln 10 M W -2ζ ln [R] -ln κ Q + 2ln S i (10)
It should also be recognised that the primary factor that controls the differing strengths of ground motions from different source mechanisms using the above theoretical model is the radiation pattern term, R φθ , that reflects the azimuthal variation of the strength of seismic waves leaving the source of a shear dislocation. In the final regression equation, factors for various styles of faulting will be included that essentially reflect the variation in the radiation pattern.

The theoretically derived model in Eq. 10 therefore provides very strong constraints upon the likely values of the parameters that should be obtained during the regression procedure. However, while the model form of Eq. 10 acts as a good guide in the selection of functional forms, there are many subtle variations from this base form that may result in better predictions of Arias Intensity. It must be remembered that many assumptions were made in the derivation of Eq. 10 and that the various functional forms finally adopted for the regression analysis will account, to some extent, for the influences of these various assumptions.

The most common modifications to theoretical magnitude and distance scaling adopted for the development of predictive equations are to include higher order scaling with respect to magnitude and to modify the distance scaling to account for finite-fault effects and/or magnitudedependent geometric spreading. There are many possible forms that these modifying terms may take and consequently, a large number of potential functional forms were initially analysed. All of these models began with Eq. 10 as a foundation and then considered alternative ways of accounting for some of the simplifying assumptions that were made in order to arrive at this fundamental form. In particular, [START_REF] Travasarou | Empirical attenuation relationship for Arias intensity[END_REF] found some dependence upon a stress parameter that resulted in nonlinear magnitude scaling that is at odds with the form of Eq. 10. For this reason, while linear magnitude scaling was initially preferred, models with higher order terms were considered and included within the final suite. The other key effect that must be considered is that Eq. 10 is based upon a point-source model and some adjustment must therefore be made to account for finite-source effects. There are numerous ways of accounting for such effects in terms of different functional forms, but it is often difficult to identify the most suitable form, as the functional terms that are introduced are often correlated with other terms, and making a distinction between what is physical and what is an artefact of the dataset is not straightforward. For example, what may appear to be magnitude saturation may, in part, be due to nonlinear site response, and what may appear to be magnitudedependent geometric spreading may be the result of finite-fault effects. The development of Eq. 10 is therefore indispensible for constraining and guiding the functional forms considered during the model development process.

Once the functional forms were identified, the coefficients for each model were determined using the nonlinear random effects procedure of [START_REF] Lindstrom | Nonlinear mixed effects models for repeated measures data[END_REF], and the associated log-likelihood, AIC and BIC were calculated. Four distinct datasets were considered, with the difference between the four depending upon the amount of foreign data included in order to constrain the scaling of the models at short sourcesite distances and large magnitudes (Table 2). The statistical significance of the model coefficients and the correlations among these coefficients were all checked, and the models were then rated according to the Akaike and Bayesian Informa-tion Criteria. The final four models that were selected following these considerations are given in Eqs. 11, 12, 13 and 14 and are hereafter referred to as Models 1-4.

lnI a,ij = c 1 + c 2 M w,i + c 3 ln r x,ij + exp c 4 M w,i +c 5 Z hyp,i + c 6 S C.ij + c 7 S D,ij + c 8 F R,i +η i + ε ij (11) lnI a,ij = c 1 + c 2 M w,i + c 3 ln r x,ij + exp c 4 M w,i +c 5 Z hyp,i +c 6 S C,ij + c 7 +c 8 lnI a,rock,ij S D,ij +c 9 F R,i + η i + ε ij (12) lnI a,ij = c 1 + c 2 M w,i + c 3 M w,i -6.5 2 +c 4 ln r x,ij + c 5 + c 6 Z hyp,i + c 7 S C,ij +c 8 S D,ij + c 9 F R,i + η i + ε ij (13) lnI a,ij = c 1 + c 2 M w,i -6 2 + c 3 ln M w,i 6 +c 4 ln r 2 x,ij + c 2 5 + c 6 Z hyp,i + c 7 S C,ij +c 8 S D,ij + c 9 F R,i + η i + ε ij (14)
In these equations r x may represent either the rupture distance, r rup , or the Joyner-Boore distance, r jb , as coefficients are obtained for both metrics. The term Z hyp is the hypocentral depth, S c and S D are dummy variables for site classes C and D of NZS1170.5:2004(Standards New Zealand 2004) and take values of one when that condition applies and zero otherwise, F R is a similar dummy variable that corresponds to reverse and reverse-oblique fault mechanisms. The coefficients are given by the terms c 1 -c 9 , the eventspecific fixed effects are given by η, and the model error is given by ε. In all cases, the subscripts i and j are indices for event and record, respectively. The base model, i.e. the case where all dummy variables are zero, represents the case of ground motions from a strike-slip or normal faulting event being observed on site class A or B of NZS1170.5:2004. In the above equations, the Arias Intensity that is predicted may be associated with a random horizontal component of ground motion (RN), the larger of two orthogonal components of ground motion (MX), the geometric mean of two orthogonal components (GM) or the arithmetic mean of these components (AM).

The η i and ε ij are assumed to be independent and normally distributed with zero means. During the regression analysis, it was found that the use of site classes rather than any continuous predictor variable for site response resulted in heteroskedastic variance with respect to site class. For this reason, the variance structure of the models includes different intra-event variability for sites on rock (site classes A and B) and sites on soil (site classes C and D). The standard deviation of the models may therefore be written as:

σ T = τ 2 + (1-S C ) (1-S D ) σ 2 rock + (S C + S D ) σ 2 soil ( 15
)
where σ T represents the total standard deviation of the model, τ 2 is the inter-event variance corresponding to the η i values and σ 2 rock and σ 2 soil are the intra-event variances corresponding to the ε ij values for rock and soil, respectively.

Results

Model coefficients have been obtained for a large number of models: four different component definitions and two different distance measures for each of the four selected functional forms. For practical purposes, it is prudent to identify one model that may be used in situations where a single estimate of Arias Intensity is required. In this case, the recommended model is Model 2, Eq. 12, in which the Joyner-Boore distance measure is used along with the arithmetic mean component definition. The regression coefficients for all of the models are given in Tables 4,5, 6 and 7, and in these tables, the variance components for each model are provided. Sometimes, people judge or evaluate alternative models by seeking to find the model that has the smallest standard deviation, but such a method of judgement does not take into account how complicated the model has to be in order to achieve this standard deviation. A preferable approach is to use the statistics of the log-likelihood, AIC and BIC in order to discern between models. Generally speaking, the higher the log-likelihood and the lower the AIC and BIC values for a given model and dataset, the better the performance of the model. Furthermore, one may use a likelihood-ratio test in order to make comparisons among alternative models and to assess whether or not adding some complexity to a model is statistically justified. Table 8 presents the log-likelihood values, AIC and BIC values for the four alternative models presented in Eqs. 11, 12, 13 and 14 for the case where the arithmetic mean of the two horizontal components is predicted.

From this table, it is clear that, at least in a statistical sense, Model 2 should be preferred over the others. This model also makes use of a functional form that is very closely related to the theoretical form in Eq. 10, and for this reason, we have some confidence that it will perform well near the limits of the range of applicability of the model. An example of the fit of the models to the observed data is given in Fig. 2 in which Model 1 is plotted along with observed data covering a relatively broad magnitude range of 0.5 magnitude units and covering all site and mechanism classes. It can be appreciated that, despite the considerable variability associated with considering such a wide range of magnitude, mechanism and site classes, the fit to the observed data is good across the entire distance range considered. The goodness of fit that is observed in Fig. 2 is not just applicable to Model 1 as can be appreciated from a consideration of the magnitude and distance scaling of all models as shown in Fig. 3. For comparative purposes, the model of [START_REF] Travasarou | Empirical attenuation relationship for Arias intensity[END_REF] is also plotted on this figure, and con- The statistics for the Joyner-Boore distance metric relate to the dataset having 385 records from 64 earthquakes while those for the rupture distance metric relate to the dataset having 344 records from 64 earthquakes df Degrees of freedom, AIC Akaike Information Criterion, BIC Bayesian Information Criterion sequently, the predictions are for the arithmetic mean of two horizontal components and for the rupture distance in order to make this comparison on a consistent basis.

Given the large suite of models that have been developed as part of this study, it is not feasible to present residual plots for all of the models. For this reason, only the residuals for the arithmetic mean are presented here in Fig. 4. In this figure, both the inter-event and intra-event components of the total error are considered with the inter-event residuals (the fixed effects) being plotted against magnitude and hypocentral depth while the intra-event residuals are plotted against magnitude and distance.

Discussion

Two, perhaps, unusual characteristics of the derived models immediately present themselves and require some justification. The first of these relates Fig. 2 Comparison between the observed and predicted ground motion for the preferred model (Model 1). The panel on the left plots events having magnitudes between 6.0 and 6.5, while the panel on the right plots events with magnitudes between 6.5 and 7.0. In both cases, the modelled curves (median and plus/minus one standard deviation) are plotted for the intermediate site class, C, and for strike-slip/normal, SN, mechanism Fig. 3 Magnitude and distance scaling of the predictive equations for the arithmetic mean of Arias Intensity. In both panels, the legend shows the line types associated with the different models. In the left panel, blue curves correspond to a rupture distance of 1 km, red to 30 km and green to 80 km. In the right panel, blue curves correspond to a moment magnitude of 5.5, red to 6.5, and green to 7.5. Also plotted are the curves for the Travasarou et al. (2003) model (TBA03) to the grouping of both site classes A and B into a single site class. Site class A corresponds to strong, hard rock sites for which the average shear wave velocity over the upper 30 m is in excess of 1,500 m/s. However, given that much New Zealand is geologically rather young and is undergoing considerable active deformation, the occurrence of such sites is relatively rare. Such sites are also very rare among the foreign records included within the dataset.

Site class B covers a very large range of shearwave velocity values, all of which can generically be referred to as rock. Given that there is no statistical basis with which to discard the records from site class A and additionally considering the relatively broad classification of rock sites proposed in NZS1170.5:2004, it makes sense to group these two classes together. The consequence of the rather crude site classification scheme is the variance structure proposed in Eq. 15. From Tables 4, 5, 6 and 7, one may appreciate that the total standard deviations corresponding to rock site conditions (site classes A and B) are significantly larger than their soil counterparts. It appears as though a significant reduction in the overall standard deviation of the model may be achieved if a more precise site classification scheme were to be adopted such as one making direct use of a continuous variable like the average shear-wave velocity over the upper 30 m. However, as the derived equations are intended for use within New Zealand, it makes sense to ensure that the site classification is consistent with the New Zealand Standard NZS1170.5:2004.

The second justification that must be made is with respect to the omission of a dummy variable to discriminate between ground motions from Fig. 4 Residual plots for Models 1-4 for the arithmetic mean. The panels on the left show the inter-event residuals plotted against both magnitude and hypocentral depth while the panels on the right show the intra-event residuals plotted against both magnitude and the Joyner-Boore distance. Residuals, i strike-slip and normal mechanism events. All of the models that have been developed initially included such a discriminating variable, but in all cases, the coefficient that was obtained was found to be statistically insignificant at the 95% confidence level. [START_REF] Travasarou | Empirical attenuation relationship for Arias intensity[END_REF] is the only other predictive model for Arias Intensity to make any distinction between fault mechanisms, and in their study, a significant difference between ground motions from strike-slip and normal mechanism events was found. The lack of difference found in the present study may well be an artefact of the limited number of records in the dataset upon which the analysis has been performed. However, given that [START_REF] Travasarou | Empirical attenuation relationship for Arias intensity[END_REF] come to the anticipated conclusion that ground motions from normal mechanism events are lower than those from strike-slip events, the models presented in the present study will most likely provide conservative estimates of ground motions associated with normal mechanism earthquakes. It should also be noted that, when the functional form adopted by [START_REF] Travasarou | Empirical attenuation relationship for Arias intensity[END_REF] is used to derive model parameters from the datasets of this study, very few of the coefficients are found to be statistically significant.

Two additional points arise when comparing the models presented herein with that of [START_REF] Travasarou | Empirical attenuation relationship for Arias intensity[END_REF]. The first is that, for Models 1, 3 and 4, the site response is modelled using simple dummy variables without consideration of any nonlinearity that may exist. For Model 2, nonlinear site-response is included for predictions of Arias Intensity for site class D. The form of this component of the model is based upon that used by [START_REF] Abrahamson | Empirical response spectral attenuation relations for shallow crustal earthquakes[END_REF] but without the parameter representing the onset of nonlinearity, as this was found to be statistically insignificant. This same form was applied for site class C but could not be supported by the data and was therefore dropped. The result is that, when using Model 2, site response for class C is likely to be overpredicted for strong ground-shaking, while it is likely to be under-predicted for weak motions. Nonlinear site response is a widely accepted physical phenomenon that should be incorporated into empirical ground-motion models where possible. However, the reason that all but one of the models developed herein neglects nonlinear site response is that the datasets are insufficient in both quantity and quality to adequately capture this effect in a stable manner. Model 2 is proposed as first attempt to incorporate this effect, and although the model is statistically robust, the currently adopted functional form will need updating as more high-quality accelerograms are obtained. Another factor that counted against including terms for nonlinear site response in all of the models was that, as the majority of the large-magnitude shortdistance data comes from foreign earthquakes, the modelling of nonlinear site-response would be almost independent of New Zealand data.

The second additional point that must be made is with respect to the adoption of magnitude-and amplitude-independent standard deviation terms. While Travasarou et al. (2003) define magnitudedependent inter-event variability and amplitudedependent intra-event variability, the dataset used in this study was not deemed complete enough to consider the inclusion of such effects. In addition, the residual plots shown in Fig. 4 do not indicate the presence of any obvious hetereoskedasticity with respect to magnitude, and this trend was also observed with respect to predicted Arias Intensity.

Apparent rates of geometric spreading in New Zealand

A feature of the derived models of quite some significance is the apparently high rates of geometric spreading that they all predict (see Fig. 3). Again, comparing the results obtained in the present study with those of [START_REF] Travasarou | Empirical attenuation relationship for Arias intensity[END_REF] and, in particular, considering the coefficients of Model 4 for the r rup distance measure and the arithmetic component definition (Table 5), as the functional form of this model is based upon that of these authors, one finds that the rate of geometrical spreading is significantly higher, cf. -2.3683 and -1.703 (Table 5). Some of this difference can be attributed to the difference between the nearsource constraint terms being greater for the Model 4 (cf. 11.0675 and 8.78), as these coefficients are correlated, but this effect would only account for a small part of the difference observed in Fig. 3. An explanation for this difference is provided when one considers that the terms that account for geometric spreading also account for scattering and anelastic material losses along the propagation path. New Zealand's tectonic environment dictates that the crust is highly fractured, and consequently, large amounts of scattering are observed in recordings of ground motions, as manifested by long trains of coda waves [START_REF] Stafford | New empirical predictive equations for the Fourier amplitude spectrum of acceleration and Arias intensity in New Zealand[END_REF]. Further evidence to support the values of the distance scaling terms found in the present study come as a result of the theoretical considerations made earlier. Following the derivation of the base functional form in Eq. 10, it was found that Arias Intensity should decay with distance at a rate twice that predicted by near-spherical geometric spreading (2ζ ). Models of peak ground acceleration [START_REF] Zhao | Attenuation of peak ground accelerations in New Zealand earthquakes[END_REF][START_REF] Cousins | A model for the attenuation of peak ground acceleration in New Zealand earthquakes based on seismograph and accelerograph data[END_REF] and the FAS of acceleration [START_REF] Stafford | Engineering seismological studies and seismic design criteria for the Buller region[END_REF][START_REF] Stafford | New empirical predictive equations for the Fourier amplitude spectrum of acceleration and Arias intensity in New Zealand[END_REF] developed for New Zealand have also found that these ground motion indices decay with distance at rates significantly greater than that predicted by spherical spreading. The findings of the present study are therefore consistent with other predictive models of ground motions developed for use in New Zealand. It should also be noted that, although records that were influenced by volcanic path effects were not included in the present study, previous studies, such as those of [START_REF] Eberhart-Phillips | Estimating slab earthquake response spectra from a 3D Q model[END_REF] and [START_REF] Mcverry | New Zealand acceleration response spectrum attenuation relations for crustal and subduction zone earthquakes[END_REF], have demonstrated that even greater rates of decay with distance can be expected in such cases. The functional forms in Eqs. 11, 12, 13 and 14, the coefficients in Tables 4,5, 6 and 7 and the residual plots in Fig. 4 all suggest that the Arias Intensity increases as the focal depth of the earthquake increases. Such a term is not ordinarily modelled within ground-motion predictive equations. When using the Joyner-Boore distance metric, one may envisage a depth term being incorporated to account for the increased distance associated with small magnitude events whose rupture areas do not approach the surface. However, for large magnitude events, such an interpretation does not make sense if one assumes that essentially equally strong waves are radiated across the fault plane. Additionally, the increase in intensity with depth contradicts this explana-tion. The physical significance of this term can be explained by the relatively thick crust that exists in New Zealand; Table 3 indicates that rupture surfaces of the earthquakes considered in this study extend down to as deep as 25 km. For earthquakes occurring at this depth, confining pressures are very high, and given that asperities tend to occur in the vicinity of the point of rupture initiation [START_REF] Mai | Hypocenter locations in finite-source rupture models[END_REF], the increase in Arias Intensity with increasing focal depth may be interpreted as being the result of higher stress-drop events occurring at greater depths (e.g. [START_REF] Allen | Depth variation of seismic source scaling relations: implications for earthquake hazard in southeastern Australia[END_REF]).

Recently, it has been demonstrated that the depth to the top of the rupture is an important predictive parameter with motions associated with surface-rupturing events leading to weaker amplitudes than their buried counterparts [START_REF] Kagawa | Differences in ground motion and fault rupture process between the surface and buried fault rupture earthquakes[END_REF][START_REF] Somerville | Differences in earthquake source and ground motion characteristics between sur-face and buried faulting earthquakes[END_REF]. The dependence on focal depth that is found in this study may partly reflect this phenomenon also.

The influence of the foreign strong-motion records

While the foreign ground-motion records that are included in the analysis are important for providing constraint to the predictive models at short distances and large magnitudes (see Fig. 1), it is important to try to gauge the extent to which the incorporation of this data influences the coefficients of the predictive equations. To this end, the regression procedure was performed for composite datasets that included varying amounts of foreign data, as dictated by the maximum distance for which foreign records were considered. The number of foreign events included for four different cases were previously specified in Table 2.

The influence that the foreign data has upon the coefficients and, consequently, the form of the predictive models is demonstrated for Models 1 and 3 in Fig. 5. From inspection of this figure, it can be appreciated that Model 1 is remarkably stable to changes in the number of foreign records included in the dataset, whereas Model 3 shows greater sensitivity. The combinations of magnitude and distance for which the greatest sensitivity is observed correspond to those regions for which there is the least amount of data, in particular at large magnitudes. The attenuation rates that are observed for both models appear to be almost independent of the amount of foreign data that is used. The objective of incorporating foreign data is to constrain the scaling of the Arias Intensity from large magnitude earthquakes at short source-site distances but without influencing the rates of attenuation of Arias Intensity within New Zealand. Figure 5 suggests that this objective has largely been achieved as the attenuation rates appear very stable and the models derived for the cases where the foreign data is limited to less than 20 km only differ slightly from the models derived for the cases with greater amounts of foreign data.

Limiting the foreign data to 20 km therefore appears to provide a good balance between ensuring near-source constraint whilst maintaining rates of attenuation relevant to New Zealand.

Differences between component definitions of Arias Intensity

The significance of ensuring that appropriate component definitions are used has recently been reinforced in various studies [START_REF] Bommer | On the use of logic trees for ground-motion prediction equations in seismic-hazard analysis[END_REF][START_REF] Baker | Which spectral acceleration are you using?[END_REF][START_REF] Beyer | Relationships between median values and between aleatory variabilities for different definitions of the horizontal component of motion[END_REF][START_REF] Boore | Orientation-independent measures of ground motion[END_REF]. Each of these studies suggests methods for moving between different component definitions used in predictive equations. Rather than relying upon empirical conversion factors, it makes sense to derive different coefficients for different component definitions directly during the regression analysis. Consequently, coefficients have been derived for four horizontal-component definitions for each of the functional forms that have been considered. The component definitions that have been considered are a random component, the larger component and the geometric and arithmetic means of two components, all with respect to orthogonal horizontal components. It may well be argued in light of orientationindependent measures of ground motion [START_REF] Boore | Orientation-independent measures of ground motion[END_REF]) that the arithmetic mean is the optimal component definition to use for the prediction of Arias Intensity, as this definition is invariant with respect to instrument orientation [START_REF] Travasarou | Empirical attenuation relationship for Arias intensity[END_REF]. However, different component definitions may be more appropriate for a given situation, and it is therefore useful to derive equations for multiple definitions directly. The differences between Arias Intensity values predicted using random, geometric mean or arithmetic mean are very similar, which may be appreciated from the similarity in the coefficients presented in Tables 4,5, 6 and 7. However, the difference between the predictions based upon the larger component and the arithmetic mean (and consequently the other component definitions) are significant. The ratios between the predicted larger and the predicted arithmetic mean Arias Intensities are both magnitude-and distancedependent and vary over the approximate range of 1.15-1.3. The value of the ratio within this range tends to increase with both decreasing magnitude and distance. The differences in the variance components may be inferred directly from Tables 4, 5, 6 and 7. Given that a 30% difference between the amplitudes of predicted ground motions based upon different component definitions may exist and that, furthermore, the standard deviation of the prediction is defined in a logarithmic sense, it is important that due consideration is given to ensuring appropriate component definitions are adopted for use in probabilistic seismic hazard analyses.

The degree to which epistemic uncertainty is accounted for

Although the suite of equations that have been presented herein go some way to capturing the epistemic uncertainty related to the prediction of Arias Intensity, it must be appreciated that the epistemic uncertainty that is accounted for in this manner is conditional upon the use of a single ground-motion dataset. Although different amounts of foreign data were considered, the final models were all derived for the case where foreign data was limited to distances of less than 20 km. In addition to the uncertainty related to the selection of specific functional forms for use in predictive equations, there is also uncertainty associated with how to compile a strong-motion dataset for regression purposes in the first place.

While the sensitivity of the results to the inclusion of the foreign data has been investigated as part of this study, the sensitivity of the results to the inclusion of records of various natures has been neglected. Examples of this aspect of epistemic uncertainty include: whether an analyst includes only free-field records or free field records and records from ground or basement levels of lightweight structures; whether or not records that may potentially be influenced by topographic effects are included; whether one restricts records to just those recorded on digital recording instruments; whether restrictions are placed upon the total usable bandwidth of a record; whether limits on distance are enforced, among others.

In each of the cases mentioned above, an analyst would obtain slightly different regression coefficients, and this difference should be included within a rigorous evaluation of the overall epistemic uncertainty associated with a PSHA involving multiple empirical ground-motion models. While the actual degree to which the epistemic uncertainty is accounted for is unknown, it is reasonable to postulate that a significant amount of the true (but unknown) epistemic uncertainty may be accounted for through the use of the presented suite of equations. In the case where weights are sought for the various models, the statistical values provided in Table 8 may be used to guide this process.

Conclusions

A new suite of predictive equations for the estimation of Arias Intensity from crustal earthquakes in New Zealand has been presented for various horizontal component definitions. In a departure from common presentations of predictive equations, four models, using different functional forms, are presented so that they may be used in a logic tree formulation within a PSHA. Differences in the amplitudes of predicted ground motions associated with different component definitions are shown to be significant, and care must therefore be taken when applying the equations to ensure that the most appropriate definition for the application in hand is adopted.

The equations that have been presented are shown to predict significantly different Arias Intensity values than those predicted using the model of [START_REF] Travasarou | Empirical attenuation relationship for Arias intensity[END_REF], which was derived using a worldwide dataset of strong-motion records. In particular, the reduction in the predicted values of Arias Intensity with distance is significantly greater for the newly developed New Zealand models. This greater rate of decay with respect to distance is shown to be consistent both with theoretical considerations and with analogous coefficients obtained in predictive equations for other measures of ground-motion in New Zealand. This high rate of decay with distance is postulated to be the result of relatively large degrees of scattering due to waves passing through a highly fractured New Zealand crust.

Fig. 1

 1 Fig. 1 Magnitude-distance distribution for the composite dataset with the foreign data limited to be within a distance of 20 km. The left panel distinguishes events according to both site class (A, B, C and D) and origin (NZ New

  , and h b correspond to the fault rupture length, rupture width, and depths to the top and bottom of the rupture surface respectively. Numbers in the Refs column correspond to the following references: (1) Anderson et al. (1993), (2) Webb and Anderson (1998), (3) Dowrick and Rhoades (1998), (4) Doser et al. (1999), (5) Zhao et al. (1997), (6) Stafford (2006), (7) Dziewonski et al. (1985), (8) Dziewonski et al. (1993), (9) Dziewonski et al. (1995), (10) Dziewonski et al. (1996), (11) Ekstrom et al. (2005a), (12) Ekstrom et al. (2005b), (13) McGinty and Robinson (1999), (14) McGinty et al. (1997), (15) Robinson et al. (2003), (16) Hincapie et al. (2005), (17) Arnadottir et al. (1995), (18) Dowrick and Rhoades (2004), (19) Matcham et al. (2006), (20) Abercrombie et al. (2000), (21) Robinson and McGinty (2000), (22) Gledhill et al. (2000), (23) McGinty et al. (2005), (24) Reyners et al. (

  Depth, Z hyp[km]

Fig. 5

 5 Fig.5Influence of the foreign data on the empirical models. The top panels show the influence of the foreign data on the magnitude and distance scaling of model 1 while the bottom panels show this influence for model 3. In

Table 1

 1 New Zealand crustal earthquakes considered as part of this studyThe Mech column indicates the mechanism of the earthquake. The N recs column specifies the number of recordings with rupture and Joyner-Boore distances within 300 km. References for these events are given in Table 3 S Strike-slip, N normal, O Oblique, R Reverse

	Event no. Event name	YYMMDD HHMM Mech. Latitude ( • ) Longitude ( • ) Z hyp (km) M w N recs
	1	Seddon	660423	0649	O	-41.63	174.40	19	5.75 1
	2	Inangahua Aftershock 680524	2057	R	-41.99	171.82	11	5.78 1
	3	Maruia Springs	710813	1442	S	-42.13	172.10	9	5.70 3
	4	Opunake	741105	1038	N	-39.65	173.63	17	5.44 1
	5	Milford Sound	760504	1356	O	-44.67	167.45	10	6.51 5
	6	Unnamed	790324	2106	R	-41.94	171.63	10	5.08 1
	7	Godley River	840624	1329	S	-43.60	170.56	13	6.12 2
	8	Lake Tennyson	900210	0327	S	-42.25	172.74	9.4	5.93 3
	9	Weber 2	900513	0423	R	-40.31	176.33	13	6.37 19
	10	Hawk's Crag 1	910128	1258	R	-41.89	171.58	10	5.79 7
	11	Hawk's Crag 2	910128	1800	R	-41.90	171.67	11	5.93 10
	12	Hawk's Crag 3	910215	1048	R	-42.04	171.59	9	5.42 5
	13	Wilberforce	920330	0702	R	-43.05	171.23	5	5.50 3
	14	Arthur's Pass	940618	0325	R	-42.98	171.48	5	6.71 12
	15	60km N White Island 941215	1120	S	-37.27	177.53	12	6.31 6
	16	Offshore East Cape	950205	2251	N	-37.65	179.49	10	7.09 10
	17	East Cape Aftershock 950210	0145	N	-37.92	179.51	10	6.49 5
	18	Cass	951124	0619	O	-42.96	171.83	9	6.25 8
	19	Thompson Sound	001101	1036	R	-45.12	167.00	18	6.20 2
	20	Jackson Bay	011207	1927	R	-44.16	168.71	4	5.80 7
	21	Karamea	020504	1259	R	-41.39	172.34	5	5.77 4
	22	Fiordland	030821	1212	R	-45.13	166.93	18	7.20 21
	23	Fiordland Aftershock 030821	1412	R	-45.29	166.84	21	6.11 8

Table 2

 2 Numbers of events and records for the New Zealand and foreign crustal earthquakes

		Limiting distance for foreign records
		r rup,max	r rup,max	r rup,max	r rup,max	r jb,max	r jb,max	r jb,max	r jb,max
		10 km	20 km	30 km	40 km	10 km	20 km	30 km	40 km
	No. events	50	64	65	68	53	64	65	68
	No. NZ events	23	23	23	23	23	23	23	23
	No. foreign events	27	41	42	45	30	41	42	45
	No. records	220	344	474	591	266	385	501	604
	No. NZ records	144	144	144	144	144	144	144	144
	No. foreign records	76	200	330	447	122	241	357	460
	The shaded columns indicate the datasets used to derive the final predictive equations

Table 3

 3 Fault rupture models for the New Zealand events detailed in Table1, indexed by event number

	Event no.	M w	Strike ( • )	D i p( • )	L(km)	W(km)	h t (km)	h b (km)	References
	1	5.75	232	68	4.8	4.8	16.8	21.2	1,3,4
	2	5.78	179	49	5.0	5.0	9.1	12.9	2,3,4
	3	5.70	242	83	4.6	4.6	6.7	11.3	1,3,4
	4	5.44	195	41	3.4	3.4	15.9	18.1	2,3
	5	6.51	48	47	16.3	10.0	6.4	13.6	3,4
	6	5.08	30	60	2.2	2.2	9.0	11.0	3,5,6
	7	6.12	62	83	9.2	7.5	9.3	16.7	3,7
	8	5.93	55	89	8.0	5.0	6.0	11.0	1,3,5,13,14
	9	6.37	215	35	19.0	15.0	8.0	16.6	2,3,5,15
	10	5.79	42	30	8.0	5.0	5.0	7.5	1,3,16
	11	5.93	8	48	8.0	5.0	5.0	8.7	1,3,16
	12	5.42	213	37	3.3	3.3	8.0	10.0	3,4
	13	5.50	70	19	3.6	3.6	4.4	5.6	3,8
	14	6.71	221	47	16.0	12.0	1.0	9.8	4,5,17-21
	15	6.31	194	77	12.2	8.6	7.8	16.2	9
	16	7.09	215	58	37.7	15.1	3.6	16.4	10
	17	6.49	220	41	15.8	9.8	6.8	13.2	10
	18	6.25	176	45	9.0	8.5	4.0	10.0	5,18,19,22
	19	6.20	171	57	12.5	12.5	14.0	24.5	15
	20	5.80	48	45	5.0	3.0	2.9	5.1	23
	21	5.77	17	60	5.0	5.0	2.8	7.2	6,11,18
	22	7.20	30	30	35.0	20.0	13.0	23.0	12,24
	23	6.11	140	65	9.1	7.5	17.6	24.4	12,24,25
	The columns L, W, h t

Table 4

 4 Regression coefficients for the four considered component definitions and two distance metrics for Model 1, Eq. 11

	Component	Arithmetic mean	Geometric mean	Maximum comp.	Random comp.
		(AM)		(GM)		(MX)		(RN)
	Distance metric	r rup	r jb	r rup	r jb	r rup	r jb	r rup	r jb
	c 1	-5.6655	-6.6177	-5.7243	-6.6655	-5.4261	-6.4007	-5.7625	-6.7075
	c 2	2.3762	2.4626	2.3842	2.4699	2.3619	2.4495	2.3218	2.4125
	c 3	-3.1536	-3.0230	-3.1552	-3.0246	-3.1596	-3.0300	-3.0771	-2.9522
	c 4	0.4854	0.4947	0.4860	0.4955	0.4853	0.4945	0.4787	0.4891
	c 5	0.0555	0.0382	0.0554	0.0379	0.0570	0.0401	0.0591	0.0416
	c 6	0.4437	0.4957	0.4532	0.5005	0.4275	0.4909	0.4737	0.5321
	c 7	0.7130	0.8212	0.7246	0.8278	0.6803	0.8038	0.7363	0.8549
	c 8	0.3695	0.2143	0.3629	0.2066	0.3833	0.2315	0.3662	0.2006
	τ	0.2670	0.3191	0.2704	0.3228	0.2493	0.2992	0.2399	0.2982
	σ soil	0.9360	0.8948	0.9360	0.8950	0.9455	0.9055	0.9478	0.9128
	σ rock	1.1213	1.0705	1.1122	1.0629	1.1477	1.0961	1.1422	1.0853
	σ T,soil	0.9733	0.9500	0.9742	0.9514	0.9778	0.9536	0.9777	0.9603
	σ T,rock	1.1526	1.1170	1.1446	1.1108	1.1744	1.1362	1.1671	1.1255

Table 5

 5 Regression coefficients for the four considered component definitions and two distance metrics for Model 2, Eq. 12

	Component	Arithmetic mean	Geometric mean	Maximum comp.	Random comp.
		(AM)		(GM)		(MX)		(RN)
	Distance metric	r rup	r jb	r rup	r jb	r rup	r jb	r rup	r jb
	c 1	-5.6006	-6.7243	-5.6618	-6.7726	-5.3564	-6.5032	-5.7091	-6.8138
	c 2	2.5653	2.6639	2.5729	2.6709	2.5502	2.6495	2.5136	2.6120
	c 3	-3.4648	-3.3059	-3.4648	-3.3067	-3.4738	-3.3137	-3.3966	-3.2393
	c 4	0.4939	0.5051	0.4946	0.5058	0.4934	0.5045	0.4872	0.4991
	c 5	0.0603	0.0416	0.0602	0.0413	0.0624	0.0439	0.0655	0.0462
	c 6	0.5014	0.5495	0.5101	0.5538	0.4878	0.5461	0.5362	0.5900
	c 7	0.2258	0.4061	0.2334	0.4087	0.2146	0.4104	0.2166	0.4061
	c 8	-0.1680	-0.1473	-0.1675	-0.1472	-0.1705	-0.1480	-0.1788	-0.1578
	c 9	0.3586	0.2104	0.3521	0.2029	0.3683	0.2241	0.3493	0.1915
	τ	0.2945	0.3337	0.2975	0.3369	0.2800	0.3164	0.2709	0.3153
	σ soil	0.9055	0.8711	0.9057	0.8714	0.9135	0.8808	0.9127	0.8852
	σ rock	1.0990	1.0591	1.0901	1.0516	1.1249	1.0848	1.1188	1.0731
	σ T,soil	0.9522	0.9328	0.9533	0.9342	0.9555	0.9359	0.9520	0.9397
	σ T,rock	1.1378	1.1104	1.1300	1.1042	1.1592	1.1300	1.1512	1.1185

Table 6

 6 Regression coefficients for the four considered component definitions and two distance metrics for Model 3, Eq. 13

	Component	Arithmetic mean	Geometric mean	Maximum comp.	Random comp.
		(AM)		(GM)		(MX)		(RN)
	Distance metric	r rup	r jb	r rup	r jb	r rup	r jb	r rup	r jb
	c 1	0.2934	0.2486	0.3034	0.2667	0.4233	0.3714	0.1166	0.0491
	c 2	1.5569	1.5835	1.5585	1.5836	1.5534	1.5832	1.5505	1.5731
	c 3	-0.4445	-0.4037	-0.4583	-0.4170	-0.4051	-0.3626	-0.4581	-0.4183
	c 4	-3.2566	-3.2395	-3.2625	-3.2443	-3.2571	-3.2451	-3.2176	-3.1938
	c 5	24.9412	28.1067	25.1254	28.2828	24.8510	28.1051	25.0436	27.9669
	c 6	0.0568	0.0437	0.0568	0.0436	0.0581	0.0454	0.0609	0.0479
	c 7	0.4348	0.4771	0.4425	0.4805	0.4222	0.4753	0.4602	0.5100
	c 8	0.7210	0.8165	0.7304	0.8214	0.6941	0.8036	0.7431	0.8502
	c 9	0.2790	0.1360	0.2719	0.1280	0.2898	0.1513	0.2680	0.1228
	τ	0.4346	0.4590	0.4378	0.4621	0.4218	0.4454	0.4071	0.4328
	σ soil	0.9217	0.8768	0.9217	0.8772	0.9310	0.8864	0.9294	0.8911
	σ rock	1.0864	1.0250	1.0766	1.0164	1.1137	1.0522	1.1080	1.0409
	σ T,soil	1.0190	0.9897	1.0204	0.9914	1.0221	0.9920	1.0147	0.9906
	σ T,rock	1.1702	1.1231	1.1622	1.1165	1.1909	1.1426	1.1804	1.1273

Table 7

 7 Regression coefficients for the four considered component definitions and two distance metrics for Model 4, Eq. 14

	Component	Arithmetic mean	Geometric mean	Maximum comp.	Random comp.
		(AM)		(GM)		(MX)		(RN)
	Distance metric	r rup	r jb	r rup	r jb	r rup	r jb	r rup	r jb
	c 1	4.4537	4.0566	4.4444	4.0481	4.5775	4.1823	4.2224	3.8539
	c 2	-4.4454	-3.9741	-4.6038	-4.1342	-3.9811	-3.4629	-4.6126	-4.1664
	c 3	38.4435	35.5813	39.4813	36.6183	35.4177	32.2812	39.4721	36.7613
	c 4	-2.3683	-2.2610	-2.3701	-2.2621	-2.3688	-2.2651	-2.3248	-2.2243
	c 5	11.0675	10.9351	11.1298	10.9791	11.0504	10.9796	10.8620	10.8030
	c 6	0.0598	0.0469	0.0597	0.0468	0.0610	0.0485	0.0637	0.0509
	c 7	0.4537	0.4951	0.4612	0.4983	0.4416	0.4936	0.4799	0.5290
	c 8	0.7442	0.8409	0.7538	0.8461	0.7174	0.8278	0.7670	0.8753
	c 9	0.3095	0.1829	0.3023	0.1752	0.3199	0.1979	0.2987	0.1677
	τ	0.4419	0.4687	0.4455	0.4729	0.4259	0.4500	0.4159	0.4445
	σ soil	0.9330	0.8977	0.9327	0.8974	0.9442	0.9097	0.9385	0.9092
	σ rock	1.0964	1.0474	1.0869	1.0388	1.1235	1.0754	1.1194	1.0636
	σ T,soil	1.0324	1.0127	1.0337	1.0144	1.0359	1.0149	1.0266	1.0120
	σ T,rock	1.1821	1.1474	1.1747	1.1414	1.2015	1.1657	1.1942	1.1527

Table 8

 8 Summary statistics for the final four models presented inEqs. 11, 12, 13 and 14

	Model	df	Joyner-Boore distance, r jb		Rupture distance, r rup
			AIC	BIC	Log-likelihood	AIC	BIC	Log-likelihood
	1	11	1,111.425	1,154.911	-544.7126	1,018.636	1,060.883	-498.3180
	2	12	1,100.067	1,147.506	-538.0334	1,004.779	1,050.867	-490.3897
	3	12	1,107.305	1,154.744	-541.6527	1,023.910	1,069.998	-499.9550
	4	12	1,124.624	1,172.063	-550.3121	1,031.664	1,077.751	-503.8318

Acknowledgements

The authors would like to acknowledge the generous financial support provided by Solid Energy New Zealand Limited, without whom, this research would not have been possible. We would also like to acknowledge the New Zealand GeoNet project and its sponsors EQC, GNS Science and FRST, for providing strong-motion data used in this study. Graeme McVerry is also thanked for providing some metadata for New Zealand records -in particular, some of the site classifications used in this study were checked against those provided by him.