PROBABILISTIC SEISMIC HAZARD ANALYSIS OF PEAK GROUND ACCELERATION FOR MAJOR REGIONAL NEW ZEALAND LOCATIONS

This paper presents site-specific probabilistic seismic hazard analysis (PSHA) results at 24 locations throughout New Zealand (NZ). Specifically, peak ground acceleration (PGA) hazard curves for two generic soft soil conditions are considered. For specific return periods of interest, seismic hazard disaggregation is used to obtain the percentage contributions of various seismic sources to the hazard, including metrics such as mean earthquake magnitude used for simplified geotechnical calculations. The seismic hazard analyses utilise concensus models for seismic source and ground-motion characterisation, including consideration of alternative ground-motion models. The analyses therefore represent an appreciable improvement relative to the science that underpin current loading standards [e.g., 1,2]. Consequently, we advocate the use of these results as a scientific basis for potential revisions to standards and guidance documents that characterise seismic hazard via PGA.


INTRODUCTION
The PSHA method has formed the internationally-accepted scientific approach for the consideration of seismic hazards for the past 60 years [3]. Seismic design standards in New Zealand (NZ) [e.g., 1,2] are based PSHA outputs [4], with adjustments in numerical values based on adopted approaches to facilitate practical implementation.
Of the myriad of different ground-motion intensity measures to quantify seismic hazard in a simple manner, PGA remains a commonly utilised parameter, particularly in geotechnical earthquake engineering problems, for example, simplified liquefaction triggering assessments [5]. For the majority of structures and infrastructure in NZ, PGA values for seismic design (and assessment) are specified via standards and guidance documents, with a small proportion of cases making direct use of a sitespecific PSHA [6]. The two most common approaches adopted in NZ practice over the past decade are to obtain PGA values from: (i) the NZS1170.5 [1] design response spectra (SA) for a vibration period of T = 0 (i.e., SA(T = 0) = PGA), or (ii) directly from the NZ Transport Agency Bridge Manual (NZTA-BM) [2], as described in the NZ Geotechnical Society (NZGS) guidelines [7].
Comparisons of the numerical values obtained from these two guidance 1 documents for magnitude-weighted PGA illustrates significant variations across NZ [8], despite the fact that they should be based on the same underlying PSHA approach. Additionally, the scientific models for earthquake sources and ground-motion shaking in these documents are also nearly 20 years old [4], and therefore do not reflect current understanding in the scientific field of seismic hazard analysis.
Given the aforementioned issues, this paper documents research that was undertaken to obtain the PGA seismic hazard using site-specific PSHA for multiple locations throughout NZ. The analyses were undertaken in the context of obtaining results for direct use in geotechnical earthquake engineering calculations; specifically simplified liquefaction triggering assessments. That said, these results are also applicable (keeping in mind the subsequently stated assumptions and limitations) in a general sense to understand the spatial variation in seismic hazard across NZ as well as the deviation between present-day understanding and that of current guidance documents, whose scientific basis is nearly two decades old.
The subsequent section discusses the seismic hazard analysis details, specifically the PSHA and disaggregation calculations and underpinning seismic source and ground-motion models adopted. The 24 selected locations are then discussed, and hazard analysis results presented at several example locations. This is proceeded by a high-level discussion of the aggregated results for all 24 locations. Finally, the implications of this study are summarised.
A companion paper [8] provides further interpretation of the results from this hazard study for use in geotechnical assessment and design.

PSHA Equations
The seismic hazard curves for PGA were obtained from conventional PSHA, which can be summarised via the equation: where λ (PGA > x) is the mean annual rate of a ground motion with PGA exceeding the value x; λ (rup i ) is the rate of occurrence of earthquake rupture rup i , provided by the Seismic Source Model (SSM); P(PGA > x|rup i ) is the probability that PGA > x given rup i , provided by a Ground-Motion Model (GMM); and the summation represents the aggregation over all N Rup ruptures in the SSM that pose a hazard to the site of interest.
Both the SSM and GMM contain uncertainty, and it is common to consider multiple models to address this modelling (epistemic) uncertainty explicitly. It is also common to consider the hazard curve outputs in the form of Return Period (RP), which is related to annual rate of exceedance via RP(PGA Because seismic hazard curves represent the aggregate contribution of all potential earthquake ruptures at a given location, then it is also of interest to disaggregate the total seismic hazard to determine the relative contribution of individual ruptures. The probability of each rupture can be determined from: Further details on seismic hazard analysis principles can be found in Kramer [9] and Baker et al. [10], while Bradley [6] provides discussion on the benefits of its use in a NZ context.

Seismic Source Model (SSM)
The SSM describes the location, characteristics, and rate of occurrence of all earthquake ruptures in the region of interest.
Characteristics include details such as, e.g., tectonic type, magnitude, and focal mechanism, among others. For application in a regional context of NZ, the model of Stirling et al. [11] provides the most recent nationwide consensus model, and was adopted in this study 2 . It is acknowledged that there continues to be advancements in SSM components (e.g., active faulting [13], time-dependent treatment of distributed seismicity [14], and insights from retrospective and prospective testing [15]), and a revision to develop a supersceeding model has recently commenced. Rather than attempt to make incremental changes to this 2010 concensus model given the science advancements in the intervening time period, which would naturally be subject to proponents and opponents, it was preferred to retain the use of a consensus model that already represents an appreciable science advancement over those models used in the development of Standards New Zealand [1], and consequently NZ Transport Agency [2].

Ground-Motion Models (GMMs)
GMMs provide the probability distribution of ground-motion intensity (PGA in the context of this study) at a location of interest resulting from a specific rupture. The variation in available GMMs is known to have an appreciable effect on PSHA results, and therefore it is conventional to consider multiple such models in order to (partially) reflect epistemic (modelling) uncertainty [10]. In this specific study, six different GMMs were considered across the three tectonic environments of relevance. For Active Shallow Crustal tectonic sources, the adopted models were Bradley [16], Abrahamson et al. [17], Campbell and Bozorgnia [18], and Boore et al. [19].  [20] and Abrahamson et al. [21] were adopted, with respective weights of 25% and 75%. The basis for these model weights draws upon the performance of these models against NZ validation data, as well as the recognised extrapolation required beyond the validation data to the seismic scenarios that dominate hazard, as discussed in Bradley [22],Van Houtte [23], among others.

Guidelines and Standards for Comparison
Given the scope of this study, comparison is primarily made with two principal documents -NZTA-BM [2] and NZS1170.5 [1]. The NZTA-BM provides guidance for both structural design of bridges, and also geotechnical-specific guidance for liquefaction and the stability and displacement of soil structures. NZS1170.5 is a structural loading standard, but has also historically been used for determining PGA values for geotechnical applications (see Cubrinovski et al. [8] for further discussion).
While these two guidance documents are routinely used in practice, it is important to note that, in the authors' view, they both contain scientific aspects that make them sub-optimal for the determination of PGA values (and associated earthquake moment magnitudes) for geotechnical applications. Firstly, NZS1170.5 explicitly notes its focus on structural applications, implicitly alluding to the fact that the parametric functions presented for determining response spectral design loadings are not intended to be used for computing PGA values. NZS1170.5 is also based on the 'larger component' definition of horizontal groundmotion intensity [e.g., 24], as opposed to the 'geometric mean' (or rotation-independent 50 th percentile) definition that is the basis for liquefaction triggering models [e.g., 5,25]. Finally, NZS1170.5 adopts a 'magnitude weighting' definition in order to reduce the 'peaked' spectra at short vibration periods that resulted from the adopted GMMs in the hazard calculation at the time [22,26] (see discussion in Cubrinovski et al. [8]). McVerry [4] provides further details on the manner in which NZS1170.5 is obtained from underpinning PSHA calculations.
In a similar vein, NZTA-BM contains 'updated' PSHA results from those in NZS1170.5 that are intended for use in geotechnical applications. Un-weighted PGA values are provided, in order to allow geotechnical analysts to apply magnitude scaling directly. Magnitude values are also provided as a so-called effective magnitude. Unfortunately, the methodology behind the calculation of the seismic hazard PGA values and effective magnitudes are not elaborated upon in NZTA-BM, and as we will show through comparison with site-specific hazard analysis values in this paper, these effective magnitude values deviate significantly from the disaggregation mean magnitude values that should ultimately be used for geotechnical calculations (See Appendix).
Despite the limitations outlined in the above two paragraphs, subsequent figures illustrate comparisons between the obtained site-specific PSHA results computed for this study and those of NZTA-BM and NZS1170. 5. In comparing such figures, it is emphasised that the site-specific PSHA and NZTA-BM re-  Figure 1 illustrates the 24 locations at which PGA hazard curves were computed using site-specific PSHA, and Table 1 provides the specific parameter values used in NZS1170.5 and NZTA-BM for computing the equivalent PGA hazard values based on these documents.

Site Conditions
In order to conduct the PSHA at each of the locations in Table 1, definitions of the surficial soil conditions are required. The majority of the GMMs utilise the 30 m average shear-wave velocity, V S, 30 , to describe the shallow site conditions. The Zhao et al. model uses a discrete site classification scheme, which was mapped to corresponding V S,30 values as described by Bradley [16].
The spatial variation of V S, 30 values in each region based on available models [27] largely resides between V S,30 = 200−300 m/s. In addition, considering the primary application (geotechnical assessment) of this study, two different reference site conditions of V S, 30 = 200 and 300 m/s were selected for the PSHA calculations. The subsequently presented figures illustrate that there are relatively minor differences between the PSHA results for these two different V S, 30 values, relative to the differences of the PHSA results to those of the existing guidance documents. Therefore, primary focus is given to the V S, 30 = 300 m/s results, with further justification provided in Cubrinovski et al. [8].

Spatial Discretisation of Locations
Another aspect considered was the degree to which the seismic hazard varies over spatial distances on the order of a few km -of relevance for the footprint of NZ's largest cities. Spatial variations in site-specific seismic hazard are affected by: (i) the variation in soil conditions that occur; and (ii) the variation in source-to-site distances to earthquake ruptures (i.e., faults). The prior subsection already addressed the consideration of two different generic descriptors of soil conditions, and therefore the remaining effect is that associated with varying source-to-site distances. Sensitivity analyses illustrate that over spatial distances of 1-2 km, without changing soil conditions, the seismic hazard varies by a maximum of a few percent for the highseismicity locations considered, and significantly less for lowseismicity locations that are not proximal to modelled active faults. Therefore, in the context of the focus of this study, a single geographic coordinate was adopted for each of the 24 locations of interest. This approach is also consistent with the definition of PGA values in the guidance documents used here for comparison, but this illustrates one of several differences between the generic PSHA results presented here and those of truly site-specific PSHA studies, as discussed subsequently.

HAZARD RESULTS
This section firstly presents example PSHA results at several locations to provide insight and introduce several terms to quantify the results. It is followed by an aggregate analysis of the results for all 24 sites.
In discussing results in this section, use of the phrases underestimation and over-estimation are made for the guidance document-based PGA hazard values with respect to the sitespecific PSHA results. This terminology is adopted on the basis that the site-specific PSHA results are considered the best-estimate of the true seismic hazard at the present time.       NZS1170.5 (as previously discussed).
Comparison across the three sites illustrates that, in both Napier and Wellington, the mean hazard from the site-specific PSHA results are higher than both guidance document-based hazard values, and that the NZTA-BM results are lower than those of NZS1170.5. The under-estimation of the site-specific haz-ard from the guidance documents is more pronounced at the 500 year RP, particularly in Wellington, where numerous largemagnitude sources with frequent rates of occurrence exist. In contrast, for Nelson, the site-specific and guidance documents are very similar. Other locations, that are subsequently discussed in the aggregated results, do exhibit site-specific PSHA results that are lower than those prescribed in these guidance documents. There is a relatively strong correlation between over-/under-prediction and the level of seismicity, with high-seismic hazard locations (e.g., Wellington, Napier) being under-estimated, and low-seismic hazard locations being overestimated in current guidance documents 4 .
The uncertainty due to alternative GMMs, reflected in the individual hazard branches, also has an appreciable effect on the seismic hazard. For example, the 500 year RP hazard in Wellington has hazard branch values that range from approximately 0.55 g to 0.75 g, with a mean of 0.68 g.
The seismic hazard disaggregation results 5 in panels (b) and (c) of Figures 2-4 provide the percentage contribution of different seismic sources binned into discrete magnitude and distance ranges, with the mean magnitude, M, noted in the figure caption. There is a general trend that the disaggregation for the 25 year RP is comprised of a greater number of contributing sources than the 500 year RP, particularly distributed seismicity sources 6 Figure 5 illustrates the mean site-specific hazard curves for the three previously discussed locations and two different V S,30 conditions of 200 and 300 m/s. In all three instances, the difference between the two hazard curves are relatively small, particularly for RP < 500 yrs. From a physics perspective, at the two higher seismic hazard sites of Napier and Wellington, the curves for 200 and 300 m/s can be seen to 'cross over', reflecting the fact that nonlinear soil response at large ground-motion intensities leads to near-surface deamplification of ground motion, and a subsequently smaller PGA hazard for V S,30 = 200 m/s than for V S,30 = 300 m/s. The effect of this nonlinear soil behaviour is less notable in Nelson, because of the smaller PGA values for a given exceedance rate of interest. As a result of these small differences, and attempting to balance practical simplicity, attention is focused on the results for V S,30 = 300 m/s in examining the results for all 24 sites in the next section (with further discussion of the relative effect of V S,30 in Cubrinovski et al. [8]). 4 In many such locations the design levels are purposefully greater than those from site-specific PSHA, and equal to some minimum threshold, in order to achieve a minimum level of resilience. 5 The coloring reflects the 'epsilon' value of the ground-motion intensity from each seismic source, indicating the exceedance probability from the ground-motion distribution (see Baker et al. [10,Chapter 7]). 6 Which reflect potential earthquakes on faults that are not explicitly modelled as individual sources.

Aggregate Results
The format of results discussed in the previous subsection were obtained for each of the 24 locations given in Figure 1. The paragraphs below discuss several summative metrics that describe the salient features of the seismic hazard curves and disaggregation in order to concisely discuss all sites. Figure 6 presents the 25 and 500 year return period PGA values from the site-specific PSHA compared to those from NZS1170.5 and NZTA-BM. Additionally, Figure 7 provides a direct one-to-one comparison of the 25 and 500 year return period PGA values from site-specific and NZTA-BM; NZS1170.5 PGA values are not explicitly shown because they implicitly already have magnitude weighting (see Cubrinovski et al. [8]). In Figure 7, points below the 1:1 line indicate site-specific values that are greater than the NZTA-BM values, and vice versa for those above the line. With Arthur's Pass as the one exception, Figure 7 illustrates the trend of NZTA-BM-based values underestimating the site-specific seismic hazard in high seismicity regions, and over-estimating it in low seismicity regions. This under-estimation is particularly significant, and further exacerbated by, the effective magnitudes prescribed being appreciably lower than the mean magnitudes obtained from site-specific PSHA, discussed briefly below and in Cubrinovski et al. [8].

PGA Hazard at RPs of 25 and 500 Years
Because of the multitude of factors that influence seismic hazard, combined with differences that result from the codification of direct PSHA results, it is difficult to isolate all the reasons for the lower values from NZS1170.5 and NZTA-BM compared to the PSHA results from this study in high seismicity regions. One notable factor is that many of these instances occur along the east coast of the North Island, proximal to the Hikurangi subduction interface. A comparison of the 2012-era source model used in this study [11], with the 2002-era model on which NZS1170.5 and NZTA-BM are based [28], illustrates that the total seismic moment rate (Ṁ O ) for all modelled events on the Hikurangi subduction interface areṀ 2012 19 Nm, respectively. This 50% larger seismic moment rate in the 2012-era model has a first-order effect on the higher seismic hazard predicted for sites where this subduction interface makes an appreciable contribution to the seismic hazard -which are east coast North Island, and upper South Island locations.   As alluded to in the context of the three sites examined in detail in the prior sub-section, there is a general trend of increasing mean magnitude for higher seismic hazard locations, and also with increasing RP. The principal observation in Figure 8 is the effective magnitudes for high-seismicity sites are substantially lower than the mean magnitudes from the site-specific hazard analysis. For example, Wellington has effective magnitudes of [6.2,7.1] from NZTA-BM, but mean magnitudes of [6.47,7.71] from the sitespecific hazard analysis (for the 25 and 500 year RPs). The individual contributions of earthquake ruptures to the seismic hazard for Wellington for the 500 year return period is illustrated in Figure 4c, with the Wellington fault (Hutt valley segment) (M7.53), Wairarapa fault (M8.17), and Hikurangi interface sources (M8.2+) contributing 72% of the total hazard. The same trend of NZTA-BM effective magnitudes being well below the site-specific mean magnitudes is true for Whanganui, Palmerston North, Nelson, Napier, Gisborne, and Blenheim.

Mean Magnitude (M)
As previously stated, the absence of specific documentation for the methodology by which the NZTA-BM effective magnitudes have been calculated precludes further understanding of the reasons why these effective magnitudes are so low; however, the implications for magnitude-weighted PGA values for assessment and design are significant, as discussed in Cubrinovski et al. [8].

Shape of the Seismic Hazard Curve
Documents such as NZS1170.5 use a Return Period Factor, R, to describe the variation in seismic hazard with RP (or, equivalently, its recriprocal -exceedance rate). R is a constant value for all geographic locations, for a given RP. The site-specific PSHA results in this study provide a means to examine the appropriateness of this approximation [4]. Figure 9 illustrates the PGA hazard for the 24 sites normalised by the 500-year RP value, PGA 500 , as compared with the NZS1170.5 R Factor. In essence, the figure provides a comparison of the shape of the normalized hazard curve. For all exceedance rates, the NZS1170.5 R factor is similar to the median (50 th percentile) of the 24 sites, indicating its accuracy in an average sense. However, there is a significant variation in the normalized PGA values across the 24 sites for any given exceedance probability. For λ = 0.04 (25 year RP), the normalized PGA values range from less than 0.1 to approximately 0.3; compared with the NZS1170.5 R factor of 0.25. Similarly, for λ = 4 × 10 −4 (2500 year RP), the normalised PGA values range from approximately 1.5 to 2.9, as compared to R = 1.8 in NZS1170.5.
This site-specific variation in the 'shape' of the seismic hazard curve is well-recognised [29], however, given the factors of 2-3 difference in the normalised PGA values across the 24 sites for the 25 and 2500 year RPs, the potential errors in using the simple R-factor on a site-by-site basis are substantial.

Sensitivity of PGA Hazard to Simplifying Assumptions
The previous sections have alluded to some important limitations of simplifying assumptions that are often made in prescriptive guidance documents stipulating seismic hazard values. Figure 10 provides an illustration of the potential for error in the specified PGA hazard (expressed as a multiplicative factor) as a result of three different attributes as described below.
The first factor is associated with the difference in the seismic hazard for the two V S30 site conditions of V S30 = 200 and 300 m/s. Figure 5 illustrated the effect of V S30 , and Figure 10 illustrates that an approximate 20% error (based on the ratio of the hazard for these two site conditions across the 24 sites for the 500 year return period) can be introduced by not discriminating between these two conditions. As noted in Cubrinovski et al. [8], this error can be significantly reduced through the use of a site amplification that is a function of the ground-motion intensity.
The second factor is associated with uncertainty due to groundmotion modelling. Figure 2- Figure 4 illustrated the uncertainty in the PGA hazard (for a given exceedance rate or return period) can vary by approximately 30% as a result of the different GMMs that comprise the consideration of GMM epistemic uncertainty. This is reflected in Figure 10 for the specific case of the Wellington PGA hazard for the 500 year return period.
The third factor is associated with the assumed 'shape' of the hazard curve. Figure 10 illustrates the range of error in the PGA hazard, using a single R-factor for the shape of the hazard curve, based on the values in Figure 9 for the 2500 year return period.
It is acknowledged that the results in Figure 10 do not represent a rigorous sensitivity analysis, and the relative sizes of the three variables are not directly compatible. That is, the site amplification range is based on 24 sites for RP = 500 yr, the GMM uncertainty range is for one site (Wellington) for RP = 500 yr, and the R factor range is for 24 sites for RP = 2500 yr (since the R factor is normalized about RP = 500 yr). Nonetheless, Figure 10 provides important context for the consideration of simplifying assumptions when PSHA results are used to subsequently derive parametric expressions for using seismic hazard values in seismic assessment and design prescriptions. considering GMM epistemic uncertainty, and using a constant R-factor instead of the site-specific shape of the hazard curve.

DISCUSSION AND CONCLUSIONS
This study has undertaken site-specific probabilistic seismic hazard analysis (PSHA) for two generic soft soil conditions and 24 locations across NZ to determine PGAs for use in geotechnical earthquake engineering applications, such as simplified liquefaction triggering assessments, slope stability analyses and retaining wall design. The adopted seismic source and groundmotion models reflect a notable advance in scientific understanding compared with those that underpin the existing standards and guidelines that are used in present practice in NZ (i.e., NZS1170.5 and NZTA-BM). The site-specific PGA values were compared against these guideline-based values for different locations, and trends identified with return period and level of seismicity.
In general, it can be concluded that the existing guidelines under-represent the true seismic hazard at locations in NZ with high seismicity, and also under-estimate the mean magnitude of the causative earthquake ruptures for the 500 year return period. Consequently, the design ground-motion loading for geotech-nical applications is appreciably under-estimated at these locations. The results presented in this study provide a scientific basis that could guide revisions to these guidance documents to more accurately reflect the accepted scientific understanding in NZ's highest seismicity regions.
The results in this study for 24 locations (sites) considered generic soil conditions, and ground-motion and seismic source models that were not tailored for the calculation of the seismic hazard at any one specific site. As such, while the results are site-specific in the sense that they are literally at a specific (set of) site(s), the rigour in geotechnical site characterisation, ground-motion and seismic source modelling are at one end of a continuum, and generally not at the same level as would be anticipated for PSHA results that are used to provide alternative solutions to the use of prescriptive guidance documents.