# A continuous map of near-surface S-wave attenuation in New Zealand

Chris van Houtte, Olga-Joan Ktenidou, Tam Larkin, Caroline Holden

<figure><img src="image_1.png" /></figure>

## To cite this version:

Chris van Houtte, Olga-Joan Ktenidou, Tam Larkin, Caroline Holden. A continuous map of near-surface S-wave attenuation in New Zealand. Geophysical Journal International, 2018, 213, pp.408-425. 10.1093/gji/ggx559. insu-03595913

HAL Id: insu-03595913

https://insu.hal.science/insu-03595913

Submitted on 3 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# A continuous map of near-surface S-wave attenuation in New Zealand

Chris Van Houtte,<fnref n="1" />,<fnref n="2" />,* Olga-Joan Ktenidou,<fnref n="3" />,<fnref n="4" /> Tam Larkin<fnref n="2" /> and Caroline Holden<fnref n="1" />

Accepted 2017 December 22. Received 2017 December 18; in original form 2017 March 2

## SUMMARY

Quantifying the near-surface attenuation of seismic waves at a given location can be important for seismic hazard analysis of high-frequency ground motion. This study calculates the site attenuation parameter, $κ_0$, at 41 seismograph locations in New Zealand. Combined with results of a previous study, a total of 46 $κ_0$ values are available across New Zealand. The results compare well with previous t* studies, revealing high attenuation in the volcanic arc and forearc ranges, and low attenuation in the South Island. However, for site-specific seismic hazard analyses, there is a need to calculate $κ_0$ at locations away from a seismograph location. For these situations, it is common to infer $κ_0$ from weak correlations with the shear wave velocity in the top 30 m, $V_{S30}$, or to adopt an indicative regional value. This study attempts to improve on this practice. Geostatistical models of the station-specific $κ_0$ data are developed, and continuous maps are derived using ordinary kriging. The obtained $κ_0$ maps can provide a median $κ_0$ and its uncertainty for any location in New Zealand, which may be useful for future site-specific seismic hazard analyses.

**Key words:** Body waves; Earthquake ground motions; Earthquake hazards; Seismic attenuation; Site effects.

## 1 INTRODUCTION

The expected level of rock-site ground motion due to a given earthquake rupture is often needed for use in site-specific hazard studies and to develop building codes. In many cases, an empirical ground motion model is used to predict the behaviour of an intensity measure from some source, path, and site response parameters. For the site response, the effect of the soil column has been extensively studied, however the site response of a rock site can also greatly affect surface ground motions (Steidl et al. 1996; Boore & Joyner 1997; Ktenidou & Abrahamson 2016). In the case of soil sites, modern empirical ground motion models such as the Next Generation Attenuation West-1 and West-2 models characterize the site response using the time-averaged shear wave velocity in the first 30 m below ground surface ($V_{S30}$) and the depth to a shear wave velocity horizon of 1 km s⁻<fnref n="1" /> ($Z_{1.0}$) as the predictive parameters. For rock-site response, empirical models typically use a single site factor for all rock sites (e.g. Abrahamson & Silva 1997; McVerry et al. 2006), or more recently, consider the rock and soil site response as a continuous function of $V_{S30}$ (Abrahamson et al. 2014; Boore et al. 2014; Campbell & Bozorgnia 2014; Chiou & Youngs 2014).

However, rock-site response is not currently well-constrained in empirical ground motion models. First, rock site $V_{S30}$ measurements are often unavailable or unreliable, both for model developers and for users. Second, research by Anderson et al. (1996), Douglas et al. (2009) and Laurendeau et al. (2013) shows that rock-site response cannot be well-modelled using $V_{S30}$ alone, as it also depends on the site attenuation parameter, $κ_0$, which models the attenuation of approximately vertically propagating seismic waves beneath the ground surface. $κ_0$ is often derived from the spectral decay parameter, $κ$ (Anderson & Hough 1984), which controls the steepness of the decay of earthquake Fourier amplitudes in the commonly used high-frequency filter model,

$$A(f) = A_0 \exp(-\pi \kappa f), f > f_e, \quad (1)$$

where $A_0$ is a source- and path-dependent Fourier amplitude of acceleration, f is the frequency, $f_e$ is the frequency above which the decay is approximately linear on a plot of log(A) against f and $κ$ controls the steepness of the spectral decay. In many studies, $f_e$ is chosen at around 10 Hz. General consensus in the seismological community is that the predominant physical mechanism behind $κ$ is attenuation along the wave propagation path, and thus $κ$ is typically modelled as a function of a distance variable, r, and site variable, s (Anderson 1991), that is,

$$\kappa(r, s) = \kappa_0(s) + \bar{\kappa}(r), \quad (2)$$

<footnote n="1">¹GNS Science, 1 Fairway Drive, Lower Hutt 5010, New Zealand. E-mail: c.vanhoutte@gns.cri.nz</footnote>

<footnote n="2">²Department of Civil Engineering, University of Auckland, 20 Symonds Street, Auckland 1010, New Zealand</footnote>

<footnote n="3">³Department of Engineering Science, University of Greenwich, Medway Campus, Central Avenue, Chatham Maritime, Kent ME4 4TB, UK</footnote>

<footnote n="4">⁴Institut des Sciences de la Terre, ISTerre, Université de Grenoble 1, CNRS, F-38041 Grenoble, France</footnote>

<footnote>*Pacific Region Office, GJI.</footnote>
where $\tilde{\kappa}(r)$ is the distance-dependence of $\kappa$ and represents attenuation along the ray path before seismic waves enter the near-surface material, and $\kappa_0(s)$ is the site attenuation parameter. Papageorgiou & Aki (1983) initially suggested that $\kappa_0$ may be a source effect, rather than site effect, however later Abercrombie (1997) used borehole data to show that near-surface attenuation is likely to be the dominant contributor to $\kappa_0(s)$.

If it is assumed that the spectrum of Fourier amplitudes released by the source, the seismic quality factor $Q$, and the site amplification are frequency-independent for frequencies greater than $f_e$, that the observed high-frequency decay is due to the direct wave, and that $Q$ is laterally homogenous, then eq. (2) can be rewritten as,

$$ \kappa(r, s) = \int_{\text{path}} \frac{dr}{Q(z)V(z)}, \quad (3) $$

(Hough & Anderson 1988), where $z$ is the depth and $V$ is wave velocity. Conditioned on these assumptions, $\kappa(r, s)$, hereafter denoted as $\kappa$, can be interpreted in terms of site and path attenuation alone.

The validity of these assumptions is rather contentious. First, the assumption of frequency-independent source spectra above $f_e$ may not be a good representation of reality. An $\omega$-square source model (e.g. Housner 1947; Haskell 1964; Aki 1967; Brune 1970), for which the displacement amplitudes decay according to $\omega^{-2}$ for large $\omega$, will result in flat amplitudes of the source spectrum of acceleration at frequencies greater than the source corner frequency. While analyses of the Canterbury earthquakes, for example, have found the $\omega$-square model to be, on average, an excellent representation of sources in that region (Oth & Kaiser 2014), this may not hold for all events and stations. It has been shown that the falloff rate of displacement spectra can vary strongly across the focal sphere (Madariaga 1976; Boatwright 1980; Kaneko & Shearer 2014, 2015), which means the assumptions of the $\kappa$ model may not hold for all azimuths and take-off angles. Despite the uncertain nature of the high-frequency shape of source spectra, this study assumes that an $\omega$-square source model applies to all events and stations used to calculate $\kappa$ values.

For the $\kappa$ model to represent attenuation, it is also necessary to assume that site effects above $f_e$ are negligible (Silva & Darragh 1995; Anderson et al. 1996; Ktenidou et al. 2017). Site effects can be both narrow band, for example, from a topographic effect or an impedance contrast in the subsurface geology, or broad-band amplification from an impedance gradient. Narrow band amplification can to some extent be identified by spectral ratios, hence this study calculates the horizontal-to-vertical spectral ratio from earthquake recordings (HVSR; Lermo & Chávez-García 1993) at each station before computing $\kappa$. If high-frequency peaks in the HVSR are observed, it was assumed that the model in eqs (2) and (3) was not appropriate for the particular site, and the site was subsequently discarded from the analysis. A broad-band amplification is much more difficult to detect, because it can typically only be inferred from detailed, site-specific, subsurface data. To apply the $\kappa$ model in New Zealand, where this information is unavailable, it is therefore necessary to equivocally assume that broad-band amplification at New Zealand rock sites is insignificant at high frequencies.

Another assumption of the $\kappa$ model is the frequency-independence of $Q$ above $f_e$. The belief that $Q$ does depend on frequency has a long history in seismology (Aki 1980; Lay & Wallace 1995). It is typically assumed that $Q^{-1}$ (defined as the fractional loss of energy per cycle) can be divided into two parts, an intrinsic component that is related to complex microscopic processes known collectively as ‘internal friction’, and a component arising from scattering of energy due to heterogeneities distributed within a medium (Dainty

1981). It is rare for direct-wave studies to separate the mechanisms, and $Q$ is usually modelled as an effective attenuation operator $Q_e^{-1}$. It is commonly observed that $Q_e$ depends strongly on frequency, at least between around 1 and 10 Hz. For this reason, $Q_e$ is often modelled using the parametric function

$$ Q_e(f) = Q_0 f^{\eta}, \quad (4) $$

where $Q_0$ is the value of $Q_e$ at a reference frequency, such as 1 Hz, and $\eta$ is an exponent controlling the degree of frequency-dependence. Atkinson (2012) and Mereu et al. (2013) compile values of $\eta$ for numerous attenuation studies in eastern North America, and find that most studies have $\eta$ values that range from 0.2 to 0.7. This indicates a strong dependence on frequency. However, several studies have also found that the frequency-dependence of $Q_e$ is more complicated than the simple model in eq. (4). In particular, coda $Q$ studies using deep borehole data, which are subject to fewer trade-offs than surface data, have shown that $Q_e$ has a weaker dependence on frequency above 10 Hz (Leary & Abercrombie 1994; Adams & Abercrombie 1998; Abercrombie 1998; Yoshimoto & Okada 2009), although $\eta$ values are still around 0.3–0.5. While the relationship between coda $Q$ and direct-wave $Q$ is not straight-forward (Yogomida & Benites 1995), it is nevertheless possible that direct-wave $Q$ in the lithosphere has a more mild frequency-dependence above 10 Hz, than for mid-range frequencies. Whether the ‘10 Hz transition’ is related to $Q$ attenuation or geometric attenuation is a matter for debate (Morozov 2008), however regardless of the origin, the consequences of the frequency-dependent attenuation on $\kappa$ values will be the same. Given the trade-offs between source, path and site effects, testing the sensitivity of $\kappa$ values (as defined by eq. 1) to the value of $\eta$ is difficult. $\kappa$ values calculated using eq. (1) are sensitive only to spectral shapes and not absolute amplitudes. If a model with $\eta > 0$ is fit to the high-frequency spectrum using eq. (1), the fit to the data will not change significantly, but the obtained $Q_0$ value will differ greatly from an $\eta = 0$ model (see Supporting Information Fig. S1 for a graphical example of this effect). However, even if a preferred value of $\eta > 0$ is propagated into the $\kappa_0$ calculations, extrapolating the $\kappa_0$ results to lower frequencies depends strongly on how one models the 10 Hz transition. Given these complications, this study follows the precedence of Anderson & Hough (1984) and many subsequent researchers, and assumes that $\eta = 0$ for the purposes of deriving $\kappa_0$. The recent results of Perron et al. (2017) suggest the practical implications of this assumption may be minimal. These authors show that calculating $\kappa$ in high- and low-frequency bands (on acceleration and displacement spectra respectively) yields different values of $\tilde{\kappa}(r)$ but similar values of hard-rock $\kappa_0$, for two sites in France.

Despite these seemingly limiting assumptions, the $\kappa$ model has become a popular method for calculating site attenuation for engineering applications. The $\kappa$ value itself is not of primary interest, because it is usually preferred to model path attenuation using $Q$ directly with the seismic wave velocity $V$, according to

$$ P(f) = A_0 \exp(-\pi f R/QV), \quad (5) $$

where $P(f)$ is the attenuation function, $A_0$ is a reference amplitude, $f$ is the frequency and $R$ is the distance along the ray path. However, the site attenuation parameter $\kappa_0(s)$ (from eq. 2), subsequently denoted as $\kappa_0$, has become a widely used parameter for modelling near-surface attenuation at rock sites. $\kappa_0$ is used in ground motion simulations using the stochastic method (e.g. Boore 2003; Motazedian & Atkinson 2005; Graves & Pitarka 2010), and as a ground motion prediction equation (GMPE) adjustment parameter in the host-to-target method, to account for regional differences in rock
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Map of continuous recording stations utilized from the New Zealand National Seismic Network (triangles), and supplementary stations from volcanic or regional networks (squares) that were also used in this study. Circles and lines indicate events and event-station paths, respectively.</figcaption></figure>

site attenuation between the host and target regions (Campbell 2003; Cotton *et al.* 2006; Douglas *et al.* 2006; Van Houtte *et al.* 2011). For simplicity, these applications almost always assume that S-waves are the dominant wave phase for peak motions, hence $κ_0$, as well as Q and V, correspond to attenuation, quality factor and velocity for S-waves. Recently, Douglas *et al.* (2009) and Laurendeau *et al.* (2013) corroborate the earlier work of Anderson *et al.* (1996) and show that the rock site response is potentially modelled better using $κ_0$ in conjunction with $V_{S30}$. Therefore, including $κ_0$ in empirical ground motion models may improve rock site predictions.

Although $κ_0$ is a widely used parameter in a multitude of ground-motion prediction applications, obtaining reliable data for forward prediction is often problematic. For site-specific hazard studies, $κ_0$ cannot typically be calculated from earthquake data due to the lack of instrumentation, and hence is often inferred from weak global $κ_0-V_{S30}$ correlations (Silva & Darragh 1995; Chandler *et al.* 2006; Van Houtte *et al.* 2011). While this method may have some benefit for ground motion prediction compared to excluding $κ_0$ altogether, the benefits are likely to be minimal. Before $κ_0$ can be widely implemented in empirical ground motion prediction for a given region, it is essential that it is well-defined and understood in the intended area of application. It is therefore the objective of this study to provide nationwide $κ_0$ values for New Zealand, to help facilitate the uptake of $κ_0$ as a parameter in empirical ground motion models. $κ_0$ values are calculated at individual stations across New Zealand,

before kriging the data with a spatially varying mean, to derive a continuous $κ_0$ map. The map can be used to infer $κ_0$ values for any rock site in New Zealand, which is likely to be a more reliable method of inference than the use of a $κ_0-V_{S30}$ correlation.

## 2 DATA

The majority of the data for this study come from GeoNet's New Zealand National Seismic Network (NZNSN; Petersen *et al.* 2011). At the time of writing, the NZNSN consists of 52 sites with paired strong-motion and broad-band instruments, located throughout New Zealand at approximately 100 km spacing. From this network, there are many continuously recording stations on rock sites that have been operational for up to ten years, providing a wealth of data from which $κ_0$ can be calculated. The data set for this study is obtained from 33 rock sites in the NZNSN, the locations of which are shown as triangles in Fig. 1.

To supplement the network in the North Island, where the station spacing of the NZNSN was deemed insufficient, eight stations from regional networks are also included in this study (shown as squares in Fig. 1). The NZNSN sites in this study comprise co-located continuous and strong-motion sensors, however for this study only data from the continuous broad-band seismometers are utilized. The regional networks use continuously recording short-period
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Magnitude and distance distribution of earthquakes that form the data set for this study.</figcaption></figure>

instruments. All stations in the data set sample at 100 Hz. The only site information available is NZS1170.5:2004 site classifications (Standards New Zealand 2004), according to which some stations are classified either as hard rock (class A) or soft rock (class B). However, these classifications are inferred, and no direct measurements of site properties exist for any sites in this study. In addition to the 41 stations analysed here, the Christchurch $κ_0$ results from Van Houtte et al. (2014) are used in the development of the continuous $κ_0$ map.

In accordance with the minimum magnitude recommendations of Van Houtte et al. (2014), only events with magnitudes greater than 2.5 are considered in this study. The purpose of this limit is to provide sufficient bandwidth above the source corner frequency, $f_c$, and below the maximum usable frequency, to allow the calculated $κ$ slope to be attributed solely to the site and regional attenuation without contamination by source effects. Although $f_c$ also depends on the stress parameter $Δσ$, which may exhibit systematic regional differences as well as event-to-event variability, this minimum magnitude requirement is considered sufficient for avoiding strong biases of $f_c$ on the calculated $κ$ slope. The available data span a magnitude range from $M_u$ 2.5 to $M_u$ 6.6. Only crustal events with depths less than 15 km are analysed in this study to ensure that the wave propagation paths from all events are sampling similar portions of the seismogenic crust. Additionally, source-to-site distances are limited to less than 150 km, so that path attenuation does not overshadow the site attenuation effect. With these criteria, a total of 1050 recordings form the data set for this study. The magnitudes and epicentral distances for these recordings are shown in Fig. 2.

## 3 METHOD

The Anderson & Hough (1984) method is adopted to calculate $κ_0$ (rather than other calculation methods), where the high-frequency slope of the acceleration spectra for recorded earthquake data is empirically fit with the low-pass filter model in eq. (1). The slope usually steepens as distance from the source increases, and removal of the trend with distance yields a site-specific $κ_0$. This method was denoted by Ktenidou et al. (2014), as $κ_{0,AS}$. This choice of method is to retain consistency between how $κ_0$ is calculated and how it is intended to be applied, in empirical ground motion modelling.

The details of the data processing are very similar to the method detailed in Van Houtte et al. (2014), and are only briefly repeated

here. For each recording, S-wave windows of five second duration are used and are used to calculate the Fourier amplitude spectrum (FAS). The instrument response is removed using the transfer function in the Seismic Analysis Code software (SAC; Goldstein & Snoke 2005). Five second pre-event noise windows are also calculated, and the recordings are only analysed for frequencies where the signal-to-noise ratio (SNR) is greater than three. $κ$ values are calculated by fitting the logarithm of the high-frequency acceleration spectrum, using linear regression between two frequencies, $f_c$ and $f_s$. $f_c$ is selected to be well beyond a theoretical source corner frequency (from the Brune (1970) model with $Δσ = 2$ MPa), and in general is greater than 10 Hz. Given the 100 Hz sampling of the data, $f_s$ is selected to always be lower than 40 Hz, although in some cases $f_s$ is limited by the SNR criterion. A minimum value of $f_s - f_c = 10$ Hz is applied to obtain reasonably stable regression.

To identify narrow band site effects that may bias $κ$ measurements, average HVSRs are calculated at each station using multiple earthquakes. If high-frequency amplification was detected at a given station, it was considered to violate the assumptions for calculating $κ$, and removed from the data set. An example $κ$ calculation, for a recording with a relatively low $κ$ value, is shown in Fig. 3. To obtain a single horizontal $κ$ value per recording, the horizontal component orientation of the recording is rotated at 5° increments through 180°, and $κ$ is calculated as detailed above for each rotation increment. The mean $κ$ value across all orientations is adopted as the $κ$ value for that particular recording. More details on this procedure can be found in Van Houtte et al. (2014). An example showing the degree of variability of $κ$ values due to component orientation is shown in Supporting Information Fig. S2.

## 4 CONSIDERATION OF DISTANCE AND DIRECTION DEPENDENCE

Once $κ$ values are calculated for all events recorded at a station, the site attenuation parameter, $κ_0$, is then calculated by extrapolating the trend of $κ$ with epicentral distance to zero distance (i.e. removing the $\tilde{κ}(r)$ term from eq. 2). The distance-dependence is approximately related to the intrinsic $Q$ for S-waves ($Q_S$) and $V_S$ structure along the horizontal wave-propagation path from source to site, as shown in eq. (3). For a medium with constant $Q_S$ and $V_S$, $κ$ will increase linearly with distance, and several previous studies adopt this simple form for the distance-dependence (e.g. Hough et al. 1988; Douglas et al. 2010; Van Houtte et al. 2011; Ktenidou et al. 2013). While many alternative forms for $\tilde{κ}(r)$ have also been proposed to better represent the geological structure in the region of interest, this study uses a constant $Q_S$ and $V_S$ model,

$$ \tilde{κ}(r) \approx \frac{R}{Q_S V_S}, \qquad (6) $$

where R is the epicentral distance. A more complex form for $\tilde{κ}(r)$ was considered, but was limited by depth uncertainties for many recorded events, making it too difficult to reliably trace the $Q_S$ structure being sampled by the ray paths. Therefore, a linear distance-dependence of $κ$ is adopted for simplicity and consistency.

For three stations in the Otago region of New Zealand, Fig. 4 shows the locations of events used to calculate $κ$, and plots the $κ$ results against epicentral distance. If the $κ$ data at each station are fitted using unconstrained linear least-squares regression (labeled as ‘free $Q$’ in Fig. 4), the $Q_S$ values implied by the $\tilde{κ}(r)$ slope range from 900-1600, assuming constant $V_S$ of 3.5 km s$^{-1}$. These $Q_S$ values can be compared to other $Q$ values obtained in an independent study. The 3-D South Island attenuation model of
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Example calculation of a low $κ$ value. (a) S-wave and noise windows from an earthquake recorded by the OPZ station, (b) the HVSR for OPZ averaged over 23 events, where dashed lines represent ±1 standard deviation from the median (solid line) and shaded area indicates the region where the amplification can be considered insignificant, (c) the FAS for the noise and S-wave windows in (a), and (d) $f_c$, $f_e$ and $f_x$ picks for calculating $κ$.</figcaption></figure>

Eberhart-Phillips et al. (2008) calculates the frequency-independent $Q$ structure for P-waves, $Q_P$, in Otago. In the typical frequency band-width for $κ$ calculations, 10–40 Hz, and a similar depth range to the data in this study, less than 10 km, Eberhart-Phillips et al. (2008) calculate $Q_P$ values ranging from 900–1100. Given that Eberhart-Phillips et al. (2008) calculate $Q_P$ and this study calculates $Q_S$, it is not possible to directly compare the results. Lay & Wallace (1995) and Castro et al. (1999) suggest $Q_P/Q_S$ is around 2 and 1 respectively, although it is not clear how these values should depend on frequency. Yoshimoto et al. (1993) show that $Q_P/Q_S$ does depend on frequency, and can range between 0.5–2.5 for frequencies greater than 1 Hz. The variability in the ratio means that the absolute $Q$ values are not comparable, but qualitatively the relative $Q_S$ values from this study compare well with the general regional trends of the $Q_P$ structure derived by Eberhart-Phillips et al. (2008).

The $κ_0$ values in Otago inferred from the ‘free $Q$’ regression are also indicated in Fig. 4. However, the station-specific $κ_0$ values may be sensitive to constraints associated with the linear distance-dependence model (Edwards et al. 2015; Ktenidou et al. 2017). With the assumption that the ray paths sample similar underlying $Q$ structure for all three Otago stations, their combined $κ$ data may be regressed together to calculate a better-constrained regional $\tilde{κ}(r)$. In this case, the combined regression yields an ‘average regional $Q_S$’ of 1100. If the slope of the linear trend is constrained at each station to fit this $Q$ model (indicated in Fig. 4 as ‘fixed regional $Q$’), the obtained $κ_0$ values may be different from those of the unconstrained, ‘free $Q$’ case. For the example of the SYZ station, adopting either ‘free $Q$’ or ‘fixed regional $Q$’ distance-dependence models causes changes in $κ_0$ that are nearly a factor of two, but for most stations the difference is within a factor of two.

For groups of stations where sufficient data have propagation paths sampling similar areas of the crust, such as the three Otago

stations in Fig. 4, the overall $κ_0$ is taken as the average between a ‘free $Q$’ and ‘fixed $Q$’ distance model. This accounts for some of the epistemic uncertainty in $κ_0$ due to the adopted distance-dependence model. In addition to the three Otago stations, the following stations were also combined to calculate a regionally fixed distance dependence:

(i) DCZ and PYZ, in Fiordland;

(ii) WEL and BHW, two stations at the southern tip of the North Island; and,

(iii) WAZ and VRZ, near the west coast of the North Island.

The other regions in New Zealand have insufficient spatial coverage of data and stations to allow the computation of robust regionally fixed $Q$ models, and for stations in these regions, $κ_0$ is calculated using only a ‘free $Q$’ linear-regression model.

To assist with identifying and interpreting strong $Q$ contrasts that are affecting $κ$, the results of this study are compared with available crustal $Q$ models for New Zealand (Eberhart-Phillips & Chadwick 2002; Eberhart-Phillips et al. 2005, 2008; Eberhart-Phillips & Bannister 2010; Eberhart-Phillips et al. 2015). Linear regression of the $κ$ data is only applied to events where the wave propagation paths sample relatively constant $Q$, and events that are likely to be affected by highly heterogeneous $Q$ structures are not considered when calculating $κ_0$. For example, Eberhart-Phillips et al. (2008) observe a localized, shallow, low $Q$ region in the central South Island, which they believe is likely to be due to the high seismic activity rate in this part of the brittle crust. Given that this $Q$ anomaly greatly affects the path-dependence of $κ$, events with ray paths passing through this region are omitted from the $κ_0$ calculations. The justification for this decision is shown in Supporting Information Fig. S3.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** (a) Locations of three stations in Otago, with events and event-station paths indicated by circles and lines, respectively, and (b)–(d) their $κ$ data. Dark lines in (b)–(d) indicate a ‘free *Q*’ distance-dependence model derived separately for each individual station, while lighter lines indicate a regionally fixed *Q* = 1100 model, common for all stations. Calculated $κ_0$ values based on free and fixed *Q* assumptions are indicated in the plots.</figcaption></figure>

## 5 CALCULATION OF $κ_0$ AT STATION LOCATIONS

For 33 stations in the NZNSN, there were sufficient reliable data to furnish $κ_0$ values. Additionally, $κ_0$ is calculated at eight stations from regional networks. Combined with five rock-site $κ_0$ values near the city of Christchurch, calculated in Van Houtte *et al.* (2014), this gives a total of 46 rock-site $κ_0$ values throughout the country. Table 1 lists the $κ_0$ values corresponding to each station, the gradient $dκ/dR$ and the number of events used to calculate $κ_0$. It also indicates stations for which $κ_0$ was calculated as the average of ‘free *Q*’ and ‘fixed regional *Q*’ regressions. Figs 5 and 6 show the individual station $κ$ data for the South and North Islands respectively, along with the derived $κ_0$ values and their 5–95 per cent

confidence intervals. Fig. 7 shows the residuals between the New Zealand $κ$ data and the linear distance-dependence model in eq. (2). The residuals do not have a significant trend with distance, which suggests the adopted constant *Q* assumption is a reasonable model for this data set. The histogram and normal *Q*–*Q* plots suggest that the residuals do not deviate significantly from a normal distribution. There is a slight increasing trend between the $κ$ residuals and magnitude, as shown in Fig. 7(d). If a linear trend is fit to the residual data, the *p*-value of its slope, 0.00575, would be considered statistically significant in traditional statistics, although perhaps not according to modern thinking (Benjamin *et al.* 2018; McShane *et al.* 2017). While there may be some effect of magnitude on the obtained $κ$ estimates, it is clear that the effect is very small compared
<table><caption>**Table 1.** Summary of $κ$ data from New Zealand. Refer to Fig. 1 for station locations.</caption><thead><tr><th>Station</th><th>κ<sub>0</sub> (s)</th><th>Slope dκ/dR</th><th>Number of recordings</th></tr></thead><tbody><tr><td colspan="4"><em>Fiordland</em></td></tr><tr><td rowspan="2">DCZ</td><td rowspan="2">0.013</td><td>3.62 × 10<sup>-4</sup></td><td rowspan="2">41</td></tr><tr><td>3.34 × 10<sup>-4</sup> a</td></tr><tr><td rowspan="2">PYZ</td><td rowspan="2">0.014</td><td>3.11 × 10<sup>-4</sup></td><td rowspan="2">22</td></tr><tr><td>3.34 × 10<sup>-4</sup> a</td></tr><tr><td colspan="4"><em>Otago</em></td></tr><tr><td rowspan="2">WHZ</td><td rowspan="2">0.009</td><td>1.79 × 10<sup>-4</sup></td><td rowspan="2">28</td></tr><tr><td>1.72 × 10<sup>-4</sup></td></tr><tr><td rowspan="2">SYZ</td><td rowspan="2">0.012</td><td>2.51 × 10<sup>-4</sup> a</td><td rowspan="2">18</td></tr><tr><td>2.54 × 10<sup>-4</sup></td></tr><tr><td rowspan="2">TUZ</td><td rowspan="2">0.006</td><td>2.51 × 10<sup>-4</sup> a</td><td rowspan="2">50</td></tr><tr><td>3.15 × 10<sup>-4</sup></td></tr><tr><td rowspan="2">OPZ</td><td rowspan="2">0.010</td><td>2.51 × 10<sup>-4</sup> a</td><td rowspan="2">24</td></tr><tr><td>3.15 × 10<sup>-4</sup></td></tr><tr><td colspan="4"><em>Canterbury</em></td></tr><tr><td>EAZ</td><td>0.016</td><td>3.33 × 10<sup>-4</sup></td><td>50</td></tr><tr><td>JCZ</td><td>0.027</td><td>2.94 × 10<sup>-4</sup></td><td>45</td></tr><tr><td>LBZ</td><td>0.022</td><td>2.14 × 10<sup>-5</sup></td><td>37</td></tr><tr><td>RPZ</td><td>0.021</td><td>2.93 × 10<sup>-4</sup></td><td>26</td></tr><tr><td>LTZ</td><td>0.016</td><td>3.85 × 10<sup>-4</sup></td><td>22</td></tr><tr><td>KHZ</td><td>0.031</td><td>5.16 × 10<sup>-4</sup></td><td>44</td></tr><tr><td colspan="4"><em>Banks Peninsula<sup>b</sup></em></td></tr><tr><td>CRLZ</td><td>0.032</td><td>—</td><td>—</td></tr><tr><td>MQZ</td><td>0.030</td><td>—</td><td>—</td></tr><tr><td>AKSS</td><td>0.033</td><td>—</td><td>—</td></tr><tr><td>D14C</td><td>0.025</td><td>—</td><td>—</td></tr><tr><td>MTPS</td><td>0.039</td><td>—</td><td>—</td></tr><tr><td colspan="4"><em>West Coast</em></td></tr><tr><td>FOZ</td><td>0.019</td><td>1.63 × 10<sup>-4</sup></td><td>30</td></tr><tr><td>WVZ</td><td>0.009</td><td>4.08 × 10<sup>-4</sup></td><td>19</td></tr><tr><td>INZ</td><td>0.022</td><td>3.91 × 10<sup>-4</sup></td><td>16</td></tr><tr><td>DSZ</td><td>0.023</td><td>1.91 × 10<sup>-4</sup></td><td>28</td></tr><tr><td>THZ</td><td>0.012</td><td>4.18 × 10<sup>-4</sup></td><td>50</td></tr><tr><td>NNZ</td><td>0.021</td><td>5.78 × 10<sup>-4</sup></td><td>21</td></tr><tr><td colspan="4"><em>Forearc / East Coast Ranges</em></td></tr><tr><td>PAWZ</td><td>0.052</td><td>5.93 × 10<sup>-4</sup></td><td>11</td></tr><tr><td>BFZ</td><td>0.050</td><td>5.27 × 10<sup>-4</sup></td><td>11</td></tr><tr><td>PXZ</td><td>0.044</td><td>3.87 × 10<sup>-4</sup></td><td>9</td></tr><tr><td>KNZ</td><td>0.042</td><td>9.67 × 10<sup>-4</sup></td><td>9</td></tr><tr><td>CKHZ</td><td>0.053</td><td>7.35 × 10<sup>-4</sup></td><td>15</td></tr><tr><td>PUZ</td><td>0.048</td><td>1.61 × 10<sup>-3</sup></td><td>39</td></tr><tr><td>MWZ</td><td>0.050</td><td>2.82 × 10<sup>-4</sup></td><td>18</td></tr><tr><td>MXZ</td><td>0.030</td><td>1.49 × 10<sup>-3</sup></td><td>15</td></tr><tr><td colspan="4"><em>Arc Ranges</em></td></tr><tr><td rowspan="2">WEL</td><td rowspan="2">0.032</td><td>4.10 × 10<sup>-4</sup></td><td rowspan="2">48</td></tr><tr><td>3.73 × 10<sup>-4</sup> a</td></tr><tr><td rowspan="2">BHW</td><td rowspan="2">0.034</td><td>4.77 × 10<sup>-4</sup></td><td rowspan="2">13</td></tr><tr><td>3.73 × 10<sup>-4</sup> a</td></tr><tr><td>BKZ</td><td>0.029</td><td>3.38 × 10<sup>-4</sup></td><td>29</td></tr><tr><td>URZ</td><td>0.030</td><td>3.64 × 10<sup>-4</sup></td><td>25</td></tr><tr><td colspan="4"><em>Taupo Volcanic Zone</em></td></tr><tr><td>WHTZ</td><td>0.055</td><td>4.34 × 10<sup>-4</sup></td><td>12</td></tr><tr><td>UTU</td><td>0.053</td><td>4.46 × 10<sup>-4</sup></td><td>10</td></tr><tr><td>MARZ</td><td>0.052</td><td>6.59 × 10<sup>-4</sup></td><td>17</td></tr><tr><td colspan="4"><em>Northern Districts</em></td></tr><tr><td rowspan="3">WAZ</td><td rowspan="3">0.019</td><td>1.88 × 10<sup>-4</sup></td><td rowspan="3">76</td></tr><tr><td>2.10 × 10<sup>-4</sup> a</td></tr><tr><td>3.13 × 10<sup>-4</sup></td></tr><tr><td>VRZ</td><td>0.017</td><td>2.10 × 10<sup>-4</sup> a</td><td>14</td></tr><tr><td>HIZ</td><td>0.015</td><td>2.38 × 10<sup>-4</sup></td><td>19</td></tr><tr><td>TLZ</td><td>0.016</td><td>3.24 × 10<sup>-4</sup></td><td>14</td></tr></tbody></table>

<table><caption>**Table 1 – Continued.**</caption><thead><tr><th>Station</th><th>κ<sub>0</sub> (s)</th><th>Slope dκ/dR</th><th>Number of recordings</th></tr></thead><tbody><tr><td>TGRZ</td><td>0.024</td><td>2.80 × 10<sup>-4</sup></td><td>29</td></tr><tr><td>TOZ</td><td>0.027</td><td>1.63 × 10<sup>-4</sup></td><td>26</td></tr><tr><td>WIAZ</td><td>0.030</td><td>2.20 × 10<sup>-4</sup></td><td>13</td></tr><tr><td>WCZ</td><td>0.029</td><td>5.39 × 10<sup>-4</sup></td><td>7</td></tr></tbody></table>

<sup>a</sup>Regionally fixed path dependence.

<sup>b</sup>Values calculated in Van Houtte et al. (2014).

to distance and other factors (the R-squared for a linear model of $κ_{observed} - κ_{predicted}$ versus $M_L$ is only 0.007). The trend may also be a result of underlying correlations in the data set, for example large earthquakes mostly occurring in regions of high fracture density and hence low $Q$.

The calculated $κ_0$ values are illustrated in Fig. 8 over a map of New Zealand’s surface geology, with the tectonic plate boundary between the Pacific and Australian plates also included for reference. Dashed lines are the boundaries between regions with distinctly different site attenuation characteristics, which are to be discussed in detail. There appears to be a strong regional dependence of $κ_0$, the variation of which is likely due to near-surface geology (i.e. depth less than 5 km) and the tectonic setting. Fig. 9 shows the range of $κ_0$ values observed for each region that is delineated in Fig. 8. These regions are addressed individually hereafter.

## 5.1 Otago

The $κ_0$ values in the southernmost regions of New Zealand are consistently lower than other regions of the country. The Otago region is geologically characterized by hard rock (schist) sites, and is identified by Eberhart-Phillips et al. (2008) as a region of high $Q_p$. $κ_0$ values from this region range from 0.006 to 0.012 s, similar to those observed in Eastern North America (Silva & Darragh 1995; Ktenidou & Abrahamson 2016; Ktenidou et al. 2016). This result might be significant for some engineering purposes, as the low $κ_0$ values are likely to result in stronger high-frequency ground motion in Otago compared to other regions in New Zealand.

## 5.2 Fiordland

$κ_0$ values are obtained at two stations in Fiordland, PYZ and DCZ, both of which are hard rock sites. $κ_0$ for PYZ and DCZ are 0.014 and 0.013 s, respectively, which are slightly higher than the values in Otago, but still lower than other regions around New Zealand. The implied $Q_s$ values of the $dκ/dR$ slope are also slightly lower than in Otago, a trend that is similarly observed by Eberhart-Phillips et al. (2008) in their 3-D South Island $Q_p$ model.

## 5.3 Canterbury and West Coast

For Canterbury and West Coast regions of the South Island, $κ_0$ values range from 0.009 to 0.027 s. Three stations in Canterbury, to the east of the Alpine Fault, are located on greywacke rock, and despite station separations of up to 250 km, have similar $κ_0$ values ranging from 0.016 to 0.022 s. The surface geology is more variable to the west of the Alpine Fault, and this is reflected in the larger range of $κ$ values. The $dκ/dR$ slopes are reasonably consistent to the east of the Alpine Fault, but there appears to be more variable $Q_s$ structure near the Alpine Fault, and on the West Coast.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** $\kappa$ data for South Island stations in this study, ordered roughly south-to-north, and by region (refer to Fig. 1). Dark lines indicate an unconstrained regression per station (free $Q$), while lighter lines indicate fixed regional $Q$ slopes calculated using $\kappa$ data from multiple neighbouring stations. Calculated $\kappa_0$ values are also indicated in the plots.</figcaption></figure>

## 5.4 Banks Peninsula

There are several strong motion stations on the Banks Peninsula, which recorded many events in the Canterbury earthquake sequence (Bannister & Gledhill 2012) and have been assigned $\kappa_0$ values in Van Houtte *et al.* (2014). Five of these stations are deemed to be good representations of rock sites, and $\kappa_0$ values from these stations are indicated in Figs 8 and 9 and Table 1. Observed $\kappa_0$ values are larger than other locations in Canterbury, likely due to stations being located on a highly weathered volcanic outcrop.

## 5.5 East Coast of the North Island

Stations located on the North Island's East Coast have some of the highest $\kappa_0$ values in the country. Sedimentary rock sites located on the coastal ranges and in the forearc basin have $\kappa_0$ values ranging from 0.030 to 0.053, similar to observations in Western North America (Anderson & Hough 1984; Silva & Darragh 1995), although the confidence intervals in Fig. 6 indicate that these values are poorly constrained in some places. This region was identified in Eberhart-Phillips *et al.* (2005) as having very low $Q_p$ in the
<figure><img src="image_9.png" /><figcaption>**Figure 6.** $\kappa$ data for North Island stations in this study. Refer to Fig. 1 for station locations. As with Figure 4, dark lines indicate an unconstrained regression (free $Q$), while light lines indicate fixed regional $Q$ slopes calculated using $\kappa$ data from adjacent stations. Dashed lines represent 5–95 per cent confidence intervals from the unconstrained regression. Calculated $\kappa_0$ values are indicated in the plots.</figcaption></figure>
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 7.** (a) The residuals between the $κ$ data and the model in eq. (2), including mean and standard deviation of 20 km distance bins, indicated as squares. (b) A histogram of the residuals, with the solid line indicating the corresponding Gaussian function. (c) A comparison of the residuals with theoretical normal quantiles. (d) Residuals plotted against magnitude.</figcaption></figure>

upper crust due to the high fracture density of material adjacent to the nearby Hikurangi subduction interface. Further west, stations located on the arc ranges have lower $κ_0$ values than in the forearc, ranging from 0.029 to 0.036 s. A similar increase in shallow $Q_p$ structure in the arc ranges is observed by Eberhart-Phillips *et al.* (2005).

In general, the $dκ/dR$ values are much higher on the North Island's east coast compared to those in the South Island. However, these values are highly variable. This observed variation in $dκ/dR$ may be due to strong $Q_s$ variation, or it could be due to uncertainties in the depth of seismicity, as events occurring within the subducting Pacific plate are likely to have much higher $κ$ values compared to events in the overlying Australian plate.

## 5.6 Taupo Volcanic Zone

The Taupo Volcanic Zone (TVZ), which forms the northern section of the backarc region, has long been recognized as a region of very high attenuation (e.g. Mooney 1970; Cousins *et al.* 1999; Dowrick & Rhoades 1999; Eberhart-Phillips & McVerry 2003; McVerry *et al.* 2006). While there are no NZNSN stations located in the TVZ, there are many short-period instruments from volcanic networks. $κ_0$ values could be estimated at three of these stations. While the results are poorly constrained, particularly at the UTU station, the $κ_0$ values are all greater than 0.050 s and are the largest in the country. These observations have an important influence on the derivation of the continuous $κ_0$ map, which will be addressed in the following section.

The TVZ $κ_0$ results, along with those on the North Island's East Coast, compare well with previous Q attenuation studies. The lower North Island $Q_p$ model of Eberhart-Phillips *et al.* (2005) suggests that either side of the high-$Q_p$ arc ranges, the forearc and volcanic fronts have very low $Q_p$ values. This pattern has also been observed in subduction regions of Japan (Pei *et al.* 2009).

## 5.7 Northern districts

Calculating $κ_0$ values for stations in the northern regions of New Zealand was more difficult than in the South Island, as these regions typically have low seismicity and hence have fewer data available. Four stations located on sedimentary rock sites, WAZ, VRZ, HIZ and TLZ, have very consistent attenuation properties, with $κ_0$ values ranging from 0.015 to 0.019 s. The four northernmost stations in this study, TGRZ, TOZ, WIAZ and WCZ, have very similar $κ_0$ values ranging from 0.024 to 0.030, however for WIAZ and WCZ in particular, these are derived from few events and thus are not well constrained. The $dκ/dR$ slopes are in general lower than the rest of the North Island, and are similar to South Island values.

# 6 DEVELOPMENT OF A CONTINUOUS $κ_0$ MAP

It is often the case in seismic hazard assessment that ground motion estimates are required at a site located away from a recording station where $κ_0$ can be directly computed. To facilitate the use of the
<figure><img src="image_14.png" /><figcaption>**Figure 8.** Calculated $\kappa_0$ values for the stations analysed in this study, superimposed on a geological map. The length of the bars is proportional to the $\kappa_0$ values, with a scale included for reference. Dashed lines indicate boundaries of regions with similar $\kappa_0$ values (discussed in detail within the text), while the solid black line is the boundary between the Australian and Pacific tectonic plates.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 9.** Regional differences in $\kappa_0$. Each point corresponds to an individual station $\kappa_0$ value, grouped into the regions indicated by Fig. 8.</figcaption></figure>

New Zealand $\kappa_0$ data in ground motion modelling, and subsequently seismic hazard, a continuous map is developed. In this section, the discrete station-specific $\kappa_0$ values are spatially smoothed and interpolated, to develop a continuous rock-site $\kappa_0$ map for

New Zealand. A number of spatial interpolation schemes were considered as candidates for developing the $\kappa_0$ map and ultimately, Gaussian process regression, or ‘kriging’, was preferred due to its flexibility for modelling the underlying process. Information on kriging is readily available in the literature, and the texts of Cressie (1993) and Diggle & Ribeiro (2007) are recommended.

The problem is initially formulated as follows. Each $\kappa_0$ observation is assumed to be a realization of a random variable, which has a distribution that depends on the value of an underlying spatially continuous Gaussian process $S(x)$ at the location $x_i$, that is,

$$K_{0,i} = S(x_i) + \epsilon_i : i = 1, \dots, n, \quad (7)$$

where $K_{0,i}$ is a vector of $n$ observed $\kappa_0$ values, and $\epsilon_i$ are normally distributed errors with mean of 0 and variance of $\tau^2$. $S(x)$, known as the ‘signal’, has mean $\mu$, variance $\sigma^2$ and a correlation function $\rho(u)$, where $u$ is the distance between two locations $x$ and $x'$ (Diggle & Ribeiro 2007). The objective of this study is to provide a model for $S(x)$ based on the New Zealand $\kappa_0$ data in Table 1. Developing this model is a three step process. First, the empirical spatial relationship of the data is examined, second the parameters of a theoretical spatial correlation model are estimated, and then finally the prediction of $S(x)$.
<figure><img src="image_16.png" /><figcaption>Figure 10. An example spatial correlation function (a), and its subsequent semi-variogram (b). Explanation of functions and terms is included within the text.</figcaption></figure>

## 6.1 Empirical spatial behaviour

A useful tool for investigating the behaviour of geospatial data is the semi-variogram, $\gamma$, which describes the variance of the difference between two realizations of the signal $S(x)$ at locations $x$ and $x'$,

$$ \gamma(x, x') = \frac{1}{2} \mathrm{Var}[S(x) - S(x')]. \qquad (8) $$

In the case of a stationary Gaussian process, the semi-variogram can also be expressed as

$$ \gamma(x, x') = \gamma(u) = \tau^2 + \sigma^2[(1 - \rho(u))]. \qquad (9) $$

A typical function for $\rho(u)$ monotonically decreases as the distance $u$ increases, hence $\gamma(u)$ usually monotonically increases with $u$. Semi-variograms tend to have the following features. When $u = 0$, the intercept of the semi-variogram is $\tau^2$, which is known as the ‘nugget variance’, or simply ‘nugget’. Also, when $\rho(u) = 0$ (no correlation between observations beyond a given distance), the semi-variance becomes $\tau^2 + \sigma^2$, called the ‘sill’. The ‘range’ of the semi-variogram is the distance at which $\gamma(u)$ equals the sill. Fig. 10 shows a schematic example of a semi-variogram with $\rho(u)$ decreasing according to $\exp(-u)$, with a nugget variance of 0.1 and sill of 0.8. Note that, as the correlation function is exponential, it never reaches zero and hence the range is undefined. It is typical to instead define a ‘practical range’ equal to 0.95 times the sill.

The solid line in Fig. 11 is the empirical semi-variogram for the New Zealand $\kappa_0$ data, determined using eq. (8), substituting the observations $\kappa_0$ and $\kappa'_0$ in place of $S(x)$ and $S(x')$. The data are assumed to be log-normally distributed, as it was found that a logarithmic transform of the samples brings the data much closer to Gaussian. eq. (7) is then rewritten as

$$ \log K_{0,i} = S(x_i) + \epsilon_i : i = 1, \dots, n. \qquad (10) $$

Semi-variance samples are combined into 50 km bins. There are several points of interest in this plot. First, there appears to be a small, but non-zero, nugget variance. The nugget variance has a dual interpretation as representing either the measurement error for each observation, or the short-distance spatial variation. If the nugget parameter is calculated from these $\kappa_0$ semi-variance data, then its value will be primarily determined by data from the Banks Peninsula, where inter-station distances are around 5–10 km. The increasing semi-variance plateaus at around 400 km, before significantly increasing beyond 800 km. The second increase reflects the large differences between $\kappa_0$ values at long inter-station distances, for example between Otago and the TVZ. This indicates that $S(x)$ is likely to have a mean that depends on location, known as a trend surface, which is a departure from a stationary Gaussian process.

<figure><img src="image_17.png" /><figcaption>Figure 11. Empirical semi-variograms of the logarithm of the New Zealand $\kappa_0$ data in 50 km distance bins. Solid line represents the semi-variances calculated based on interstation distance alone, while dashed line is calculated using a spatially varying mean, using an indicator variable for stations in the Taupo Volcanic Zone.</figcaption></figure>

Given that the spatially varying mean is likely due to complex variation of crustal properties, this trend surface is difficult to model geostatistically. To minimize the effect of this trend on the spatial prediction, particularly at short inter-station distances, a simple trend surface is applied to the spatial $\kappa_0$ data. The trend surface consists of a single indicator variable for the Taupo Volcanic Zone, $F_{TVZ}$, such that the mean function is

$$ \mu(x) = \beta_0 + \beta_1 \cdot F_{TVZ} \qquad (11) $$

where $F_{TVZ}$ is 1 for sites located within the TVZ, and 0 otherwise. The justification for this trend is that the TVZ is well-known for having very high near-surface attenuation properties, and its geographical extent is well-defined in Wilson et al. (1995). The dashed line in Fig. 11 shows the new empirical semi-variogram with this trend surface applied. Not only does this model lower the semi-variogram at large distances, but the semi-variances are also lower at short distances. This is because the trend surface allows stronger correlations between $\kappa_0$ from the TVZ sites and nearby stations outside the volcanic arc.

## 6.2 Semi-variogram parameter estimation

Before spatial prediction can occur, a model for the empirical data in Fig. 11 must be developed. A first step is to choose a parametric function for $\rho(u)$ that best models the empirical $\kappa_0$ semi-variances.
<figure><img src="image_18.png" /><figcaption>**Figure 12.** (a) The effect of $\theta$ on the Matérn correlation function, with a fixed value of $\phi = 0.2$. (b) The effect of $\phi$, for a $\theta$ value of 1.5.</figcaption></figure>

<table><caption>**Table 2.** Results for the two semi-variogram models described in the text. LL represents the model log-likelihood and AIC represents the Akaike Information Criterion.</caption><thead><tr><th colspan="7">Model 1 - trend surface for TVZ and free τ<sup>2</sup></th></tr><tr><th>Order</th><th>β<sub>0</sub></th><th>β<sub>1</sub></th><th>σ<sup>2</sup></th><th>τ<sup>2</sup></th><th>φ</th><th>LL</th><th>AIC</th></tr></thead><tbody><tr><td>θ = 0.5</td><td>-1.630</td><td>0.355</td><td>0.055</td><td>0.003</td><td>274.2</td><td>21.81</td><td>-33.62</td></tr><tr><td>θ = 1.5</td><td>-1.624</td><td>0.350</td><td>0.045</td><td>0.008</td><td>105.5</td><td>21.80</td><td>-33.60</td></tr><tr><td>θ = 2.5</td><td>-1.622</td><td>0.344</td><td>0.044</td><td>0.009</td><td>76.5</td><td>21.82</td><td>-33.64</td></tr></tbody><tbody><tr><th colspan="7">Model 2 - trend surface for TVZ and fixed τ<sup>2</sup> = 0.03</th></tr><tr><th>Order</th><th>β<sub>0</sub></th><th>β<sub>1</sub></th><th>σ<sup>2</sup></th><th>φ</th><th></th><th>LL</th><th>AIC</th></tr><tr><td>θ = 0.5</td><td>-1.654</td><td>0.294</td><td>0.045</td><td></td><td>646.0</td><td>14.31</td><td>-20.62</td></tr></tbody></table>

There are many different correlation functions that can be adopted to model different types of spatial correlation behaviour. In this study, a widely used family of correlation functions is adopted, known as the Matérn (1960) family, which have the general form

$$ \rho(u) = [2^{\theta-1}\Gamma(\theta)]^{-1} \left(\frac{u}{\phi}\right)^{\theta} K_{\theta}\left(\frac{u}{\phi}\right). \quad (12) $$

$K_{\theta}$ denotes a modified Bessel function with non-negative order $\theta$, $\phi$ is a non-negative scale parameter with units of distance and $\Gamma$ is the gamma function. As with previous notation, $u$ is the distance between two locations $x$ and $x'$. The advantage of this form is that it is very flexible, because $\theta$ dictates the behaviour of the correlation structure where the station separation distance is small, while $\phi$ controls the degree of correlation between stations with large separation distance. Fig. 12 gives an example of the effect of different values of $\theta$ and $\phi$ on $\rho(u)$. Note that for $\theta = 0.5$, $\rho(u)$ reduces to the exponential function $\exp(u/\phi)$.

The Matérn parametric form is used to model the spatial correlation of the $\kappa_0$ data. Treating the transformed spatial $\kappa_0$ data as Gaussian, the model can be written as

$$ \log K_0 \sim N(Y\beta, \sigma^2 B(\theta, \phi) + \tau^2 I), \quad (13) $$

where $K_0 = (\kappa_{0,1}, \dots, \kappa_{0,n})$, $Y$ is a matrix of covariates and $\beta$ the corresponding matrix of coefficients, $B$ is an $n \times n$ correlation matrix, which in this case depends on the Matérn parameters $\theta$ and $\phi$, and $I$ is the identity matrix. The corresponding log-likelihood function is

$$ \mathcal{L}(\beta, \tau, \sigma, \theta, \phi) = -0.5 \left\{ \begin{aligned} & n \log(2\pi) + \log |\sigma^2 B(\theta, \phi) + \tau^2 I| \\ & + (\log \kappa_0 - Y\beta)^T [\sigma^2 B(\theta, \phi) + \tau^2 I]^{-1} \\ & \times (\log \kappa_0 - Y\beta) \end{aligned} \right\}. \quad (14) $$

Maximizing eq. (14) gives the model parameters. However, the Matérn parameters $\theta$ and $\phi$ tend to be strongly correlated, so a

common, alternative method for parameter estimation is to fix a discrete set of $\theta$, then maximize eq. (14) to determine $\beta$, $\tau^2$, $\sigma^2$ and $\phi$ (Diggle & Ribeiro 2007).

Using this method, Table 2 shows the maximum likelihood parameters for two candidate models, and their corresponding log-likelihoods. Both models apply the trend surface from eq. (11), but the first model solves for the nugget variance, while the second model fixes the nugget variance to an indicative value of the station-specific $\kappa_0$ uncertainty. Guided by the $\kappa(r)$ regression results in Figs 5 and 6, an approximate standard deviation value of 50 per cent (0.18 log$_{10}$ units) was assumed as the average uncertainty of $\kappa_0$, and hence $\tau^2$ was fixed to equal 0.03. For the first model, three values of the Matérn order are considered, 0.5, 1.5 and 2.5. These three values are commonly adopted in geostatistical modelling to represent different smoothness of $S(x)$. Specifically, they represent processes that are mean-square continuous, once mean-square differentiable and twice mean-square differentiable respectively. For the second model, only a model with $\theta = 0.5$ is considered, as models with higher values of $\theta$ are overparametrized and give unstable solutions.

The first point from Table 2 is that $\phi$ usually decreases as $\theta$ increases, illustrating the correlation between the parameters. However, the most important effect of the $\theta$ parameter is to increase the nugget variance $\tau^2$. The models with different $\theta$ values have similar log-likelihoods, which indicates that $\tau^2$ is not well-constrained by the data. The models with free $\tau^2$ have higher log-likelihoods than the model with fixed $\tau^2$.

The best-fit theoretical semi-variograms are shown in Fig. 13. The best-fit models fit the empirical semi-variances at small separation distances better than at large distances, which is expected because eq. (14) allows for less precise estimates as distance increases. The inset in Fig. 13(a) shows the effect of the different $\theta$ values at short station-separation distances. The more complicated form of the $\theta = 1.5$ and $\theta = 2.5$ models is not well justified by the data, hence the model with $\theta = 0.5$ is preferred in this study.
<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 13.** The binned empirical variogram for the logarithm of the New Zealand $\kappa_0$ data corresponding to (a) ‘Model 1’ and (b) ‘Model 2’. Solid, dashed and dotted lines correspond to models with Matérn orders of 0.5, 1.5 and 2.5, respectively, and other parameters determined by maximum likelihood estimation. Circles in (a) and (b) correspond to the dashed lines in Fig. 11. The inset in (a) is a more detailed view of the short distance behaviour of Model 1 for different values of the Matérn order.</figcaption></figure>

## 6.3 Spatial prediction

The final step for developing a continuous site attenuation model is to predict the underlying signal $S(x)$ from the observed $\kappa_0$ realizations and the candidate covariance models. This study uses the ordinary kriging algorithm, the theory behind which is well described in the statistical literature and is only briefly repeated here. In short, ordinary kriging, involves minimizing the mean square prediction error

$$ \text{MSPE}[\hat{S}(x)] = \frac{1}{n} \sum_{i=1}^{n} [(\hat{S}(x) - S(x_i))]^2, \qquad (15) $$

where $\hat{S}(x)$ is the prediction of $S(x)$ at an unobserved location. The form of $\hat{S}(x)$ that minimizes MSPE[$\hat{S}(x)$] is,

$$ \hat{S}(x) = \mu + b'V^{-1}(K_0 - \mu), \qquad (16) $$

where $\mu$ is a vector of mean values of $S(x)$, $b$ is a vector of length $n$ with values corresponding to the correlation between the unobserved location and the $n$ observed locations, and $V = B + (\tau^2/\sigma^2)I$.

In ordinary kriging, $\mu$ is assumed to be unknown and is estimated by generalized least squares, given by

$$ \hat{\mu} = (\mathbf{1}'V^{-1}\mathbf{1})^{-1}\mathbf{1}'V^{-1}K_0. \qquad (17) $$

The prediction variance is given by

$$ \text{Var}(S(x)|K_0) = \sigma(1 - b'V^{-1}b), \qquad (18) $$

Full derivations of these expressions are readily available in the literature.

As the models contain the $F_{TVZ}$ indicator variable, it is necessary to define the boundaries of the TVZ in the prediction grid. Wilson et al. (1995) define different boundaries for the TVZ based on its geological history, and their model for the ‘whole TVZ’ was preferred by Cousins et al. (1999) to represent volcanic path attenuation within empirical ground motion models. Like Cousins et al. (1999), the ‘whole TVZ’ model is selected here as the boundary for the high attenuation volcanic region.

Fig. 14 compares the predicted median and standard deviation of $\kappa_0$ for models 1 and 2, both with $\theta$ equal to 0.5. Due to the larger $\tau^2$, model 2 predicts a smoother median surface. The effect of increasing $\tau^2$ is essentially to trade in data reproduction for variance, and hence $\hat{S}(x)$ approaches $\bar{\kappa}_0$. The standard deviations in Fig. 14 are obtained by taking the square root of the sum of the underlying signal variance and $\tau^2$. This calculation means that the maps in Figs 14(b) and (d) provide the prediction standard deviation for an unobserved $\kappa_0$ value.

## 6.4 Consideration of $\kappa_0$ ‘measurement uncertainty’

Fig. 14 clearly demonstrates that the nugget is an important parameter when considering forward application of the $\kappa_0$ map. While the seismograph network in New Zealand is insufficiently dense to derive the nugget parameter from the kriging analysis alone, estimates of $\kappa_0$ measurement uncertainty are available from the linear regression of the $\kappa(r)$ data with distance (Figs 5 and 6). This is the purpose of model 2 in Table 2.

The fixed $\tau^2$ value of 0.03 for model 5 is only an indicative value of the $\kappa_0$ uncertainty across all stations. A limitation of this assumption is that some of the $\kappa_0$ data are better constrained than others, as evidenced by the confidence intervals in Figs 5 and 6. Ideally, a fixed vector of $\tau^2$ with each station’s $\kappa_0$ uncertainty could be propagated into the maximum likelihood step. However, this calculation is complicated by the assumption of two different distributions to derive the $\kappa_0$ map. The station-specific $\kappa_0$ are derived by regression of $\kappa$ data with distance on a linear scale, while a lognormal distribution was assumed to derive the $\kappa$ map. The combination of distributions prevents an elegant determination of the spatially varying nugget, hence the fixed nugget of 0.03 ($\pm$ 50 per cent of the median value) is adopted for simplicity.

The choice of the preferred model is likely a personal preference, however we are inclined to select model 5. While the AIC of model 2 is lower than model 1, the higher value of the nugget is more representative of the data and its underlying uncertainty. Model 1 relies on spatially dense sampling of $\kappa_0$ to constrain the nugget, which are mostly unavailable in New Zealand.

For model 1, the predicted standard deviation on $\kappa_0$ is around 0.1 log$_{10}$ units, while for model 2 it is around 0.2 log$_{10}$ units. Note that this 0.2 is only slightly larger than the value that was initially assumed for the uncertainty in the station-specific $\kappa_0$ values, which suggests that the uncertainty in calculating $\kappa_0$ at a given station location is much larger than the uncertainty introduced by kriging.
<figure><img src="image_21.png" /><figcaption>**Figure 14.** (a) Median $κ_0$ map for Model 1 with $θ = 0.5$ and (b) its standard deviation in $log_{10}(s)$. (c) The median model for Model 2 with $θ = 0.5$ and (d) its standard deviation.</figcaption></figure>
An additional aspect of uncertainty in the $κ_0$ map comes from the small data set. This uncertainty can be demonstrated using jackknife resampling, also known as leave-one-out cross-validation. Maps of the jackknife standard deviation for the $κ_0$ map can be found in Supporting Information Fig. S4, although in general, the jackknife standard deviation is much smaller than the prediction standard deviation shown in Fig. 14.

# 7 SUMMARY AND DISCUSSION

The overarching objective of this study is to provide the best available information on $κ_0$ in New Zealand, to aid efforts to include it in empirical ground-motion prediction, and subsequently in hazard calculations. Regions within New Zealand are shown to have significantly different attenuation properties depending on local geology and tectonic setting, and rock site $κ_0$ can vary from 0.006 to 0.055 s. The general path attenuation trends observed in this study are consistent with 3-D $Q_p$ attenuation models of New Zealand (Eberhart-Phillips et al. 2015). Many stations have $κ_0$ values that are poorly constrained, with 19 values being derived from fewer than 20 recordings each. While it is acknowledged that these poor constraints limit the level of confidence in the interpretation of the $κ_0$ map, the proposed model nevertheless represents the best information currently available on $κ_0$ in New Zealand. The model is a first attempt at a nationwide mapping of the site attenuation parameter, and can be updated and improved as more data become available.

Two geostatistical models of the $κ_0$ data are used to derive continuous maps of $κ_0$. The purpose of presenting two models was to illustrate that very different $κ_0$ maps can be obtained with different kriging assumptions. Our preferred model is model 2. This model accounts for the measurement uncertainty in $κ_0$, albeit in a simplified manner. A spatially varying measurement uncertainty could not be easily applied, because the selection of the best distribution to represent $κ_0$ is still an unresolved issue. Hough et al. (1988), Oth et al. (2011), Edwards et al. (2015) and Van Houtte et al. (2014) suggest that $κ_0$ is normally distributed. However, the assumption of normality generates problems in forward modelling, as it introduces the possibility of unphysical, negative $κ_0$ values. As such, $κ_0$ is typically assumed to be either uniformly or lognormally distributed when it is applied in seismic hazard assessments (e.g. EPRI 1993; Silva & Darragh 1995; Atkinson & Boore 2006; Campbell et al. 2014), to prevent consideration of non-physical model parameters. This study supports these assumptions, because the population of $κ_0$ data in New Zealand are much better represented by a lognormal distribution.

Ktenidou & Abrahamson (2016) suggest that negative $κ_0$ values may be a result of site amplification at high frequencies, although still find some negative $κ_0$ values when correcting for generic ‘crustal amplification’ effects. Edwards et al. (2015) suggest that negative $κ_0$ values are likely to be a result of overly-simplistic assumptions in the derivation of $κ_0$, for example the assumption of $ω^2$ decay of the source displacement spectrum. Should regional earthquakes have displacement spectra that consistently decay with $ω^{n < 2}$, combined with low site attenuation, it is conceivable that a negative $κ_0$ value can be observed, when deriving $κ_0$ empirically from the FAS of acceleration (the method used in this study). Such behaviour of source spectra may bias empirical $κ$ observations, like those derived in this study. However, it is our opinion that if negative $κ_0$ values are being observed in a given region, then the typical assumptions for deriving $κ_0$ values are being violated, and it is necessary to re-examine the spectral behaviour of regional seismic sources. If near-surface attenuation is to be modelled in seismic hazard analysis, it is essential that the applied values have a physi-

cal interpretation. Rather than permitting non-physical attenuation values, and propagating these into seismic hazard calculations, it is preferable to develop a more appropriate physical model for the spectral characteristics.

Given that no negative $κ_0$ values are observed in this study, it is assumed that the $κ_0$ data here can be considered positive-only. It is therefore desirable to adopt a distribution for $κ_0$ that prevents negative values. Unfortunately, this is not a trivial task. A logarithmic transformation of the $κ$ data seems to be an obvious starting point. However, the results in Fig. 7 suggest that a normal distribution is a good representation of $κ$ data. While beyond the scope of this work, resolving this issue would be important for future research. When a satisfactory method for modelling the distribution of $κ_0$ is available, a more robust $κ_0$ map and its uncertainty can be derived. This will make the derived median and standard deviation models more useful in the context of seismic hazard assessment.

Lastly, this model is intended to represent rock site $κ_0$ only. Should a $κ_0$ estimate be required for a soil site, the continuous $κ_0$ model in Fig. 14 can be used in conjunction with a model for the dependence of $κ_0$ with sediment depth. The effect of the sedimentary column on $κ_0$ has not yet been widely investigated, at least in terms of observed spectral amplitudes. To the authors’ knowledge, the only available models are those of Campbell (2009) and Ktenidou et al. (2015). These models are derived using data from Eastern North America and Greece respectively. Given the very different tectonic environment and sediment types in New Zealand, it is unlikely that these models are applicable to New Zealand conditions. Therefore, the $κ_0$ map in this study can only currently be applied to rock sites, and further research is necessary to determine the effect of the soil column on the site attenuation in New Zealand.

# 8 DATA AND SOFTWARE RESOURCES

All data in this study are publicly available from GeoNet (www.geonet.org.nz), which is funded by the Earthquake Commission (EQC), Land Information New Zealand (LINZ) and GNS Science. Codes to fit the geostatistical model and predict the $κ_0$ maps can be found on github (www.github.com/cvanhoutte/kappa). Grd files for all $κ_0$ maps presented in this study can also be found at this web address. Some analysis was undertaken using the geoR package (Ribeiro & Diggle 2001), and some figures were generated using Generic Mapping Tools (GMT; Wessel et al. 2013).

# ACKNOWLEDGEMENTS

This research was partially funded by a University of Auckland Doctoral Scholarship, and by Institut des Sciences de la Terre. Peter Stafford, Gail Atkinson, Zach Eilon and an anonymous reviewer are greatly and genuinely thanked for their encouragement and comments on early versions of this manuscript, which greatly improved the final model.

# REFERENCES

* Abercrombie, R., 1997. Near-surface attenuation and site effects from comparison of surface and deep borehole recordings, *Bull. seism. Soc. Am.*, **87**(3), 731–744.

* Abercrombie, R., 1998. A summary of attenuation measurements from borehole recordings of earthquakes: the 10 Hz transition problem, *Pure appl. Geophys.*, **153**(2–4), 475–487.

* Abrahamson, N. & Silva, W., 1997. Empirical response spectral attenuation relations for shallow crustal earthquakes, *Seismol. Res. Lett.*, **68**(1), 94–127.
Abrahamson, N., Silva, W. & Kamai, R., 2014. Summary of the ASK14 ground motion relation for active crustal regions, *Earthq. Spectra*, **30**(3), 1025–1055.

Adams, D. & Abercrombie, R., 1998. Seismic attenuation above 10 Hz in southern California from coda waves recorded in the Cajon Pass borehole, *J. geophys. Res.*, **103**(B10), 24 257–24 270.

Aki, K., 1967. Scaling law of seismic spectrum, *J. geophys. Res.*, **72**(4), 1217–1231.

Aki, K., 1980. Attenuation of shear-waves in the lithosphere for frequencies from 0.05 to 25 Hz, *Phys. Earth planet. Inter.*, **21**(1), 50–60.

Anderson, J., 1991. A preliminary descriptive model for the distance dependence of the spectral decay parameter in southern California, *Bull. seism. Soc. Am.*, **81**(6), 2186–2193.

Anderson, J. & Hough, S., 1984. A model for the shape of the Fourier amplitude spectrum of acceleration at high frequencies, *Bull. seism. Soc. Am.*, **74**(5), 1969–1993.

Anderson, J., Lee, Y., Zeng, Y. & Day, S., 1996. Control of strong motion by the upper 30 meters, *Bull. seism. Soc. Am.*, **86**(6), 1749–1759.

Atkinson, G., 2012. Evaluation of attenuation models for the northeastern United States/southeastern Canada, *Seismol. Res. Lett.*, **83**(1), 166–178.

Atkinson, G. & Boore, D., 2006. Earthquake ground-motion prediction equations for eastern North America, *Bull. seism. Soc. Am.*, **96**(6), 2181–2205.

Bannister, S. & Gledhill, K., 2012. Evolution of the 2010–2012 Canterbury earthquake sequence, *N.Z. J. Geol. Geophys.*, **55**(3), 295–304.

Benjamin, D. et al., 2018. Redefine statistical significance, *Nat. Hum. Behav.*, **2**(1), 6–10.

Boatwright, J., 1980. A spectral theory for circular seismic sources; simple estimates of source dimension, dynamic stress drop, and radiated seismic energy, *Bull. seism. Soc. Am.*, **70**(1), 1–27.

Boore, D., 2003. Simulation of ground motion using the stochastic method, *Pure appl. Geophys.*, **160**(3–4), 635–676.

Boore, D. & Joyner, W., 1997. Site amplifications for generic rock sites, *Bull. seism. Soc. Am.*, **87**(2), 327–341.

Boore, D., Stewart, J., Seyhan, E. & Atkinson, G., 2014. NGA-West2 equations for predicting PGA, PGV, and 5 per cent damped PSA for shallow crustal earthquakes, *Earthq. Spectra*, **30**(3), 1057–1085.

Brune, J., 1970. Tectonic stress and the spectra of seismic shear waves from earthquakes, *J. geophys. Res.*, **75**(26), 4997–5009.

Campbell, K., 2003. Prediction of strong ground motion using the hybrid empirical method and its use in the development of ground-motion (attenuation) relations in eastern North America, *Bull. seism. Soc. Am.*, **93**(3), 1012–1033.

Campbell, K., 2009. Estimates of shear-wave *Q* and *κ*₀ for unconsolidated and semiconsolidated sediments in eastern North America, *Bull. seism. Soc. Am.*, **99**(4), 2365–2392.

Campbell, K. & Bozorgnia, Y., 2014. NGA-West2 ground motion model for the average horizontal components of PGA, PGV, and 5 per cent damped linear acceleration response spectra, *Earthq. Spectra*, **30**(3), 1087–1115.

Campbell, K., Hashash, Y., Kim, B., Kottke, A., Rathje, E., Silva, W. & Stewart, J., 2014. Reference rock site condition for Central and Eastern North America Part II—attenuation (kappa) definition, PEER Research Report 2014/12, University of California Berkeley, Berkeley, CA.

Castro, R., Monachesi, G., Mucciarelli, M., Trojani, L. & Pacor, F., 1999. *P*- and *S*-wave attenuation in the region of Marche, Italy, *Tectonophysics*, **302**(1), 123–132.

Chandler, A., Lam, N. & Tsang, H., 2006. Near-surface attenuation modelling based on rock shear-wave velocity profile, *Soil Dyn. Earthq. Eng.*, **26**(11), 1004–1014.

Chiou, B. & Youngs, R., 2014. Update of the Chiou and Youngs NGA model for the average horizontal component of peak ground motion and response spectra, *Earthq. Spectra*, **30**(3), 1117–1153.

Cotton, F., Scherbaum, F., Bommer, J. & Bungum, H., 2006. Criteria for selecting and adjusting ground-motion models for specific target regions: application to Central Europe and rock sites, *J. Seismol.*, **10**(2), 137–156.

Cousins, W. J., Zhao, J. & Perrin, N., 1999. A model for the attenuation of peak ground acceleration in New Zealand earthquakes based on seismograph and accelerograph data, *Bull. N.Z. Soc. Earthq. Eng.*, **32**(4), 193–220.

Cressie, N., 1993. *Statistics for Spatial Data*, John Wiley & Sons.

Dainty, A., 1981. A scattering model to explain seismic Q observations in the lithosphere between 1 and 30 Hz, *Geophys. Res. Lett.*, **8**(11), 1126–1128.

Diggle, P. & Ribeiro, P., 2007. *Model-based Geostatistics*, Springer Series in Geostatistics, Springer-Verlag.

Douglas, J., Bungum, H. & Scherbaum, F., 2006. Ground-motion prediction equations for southern Spain and southern Norway obtained using the composite model perspective, *J. Earthq. Eng.*, **10**(1), 33–72.

Douglas, J., Gehl, P., Bonilla, L. F., Scotti, O., Régnier, J., Duval, A.-M. & Bertrand, E., 2009. Making the most of available site information for empirical ground-motion prediction, *Bull. seism. Soc. Am.*, **99**(3), 1502–1520.

Douglas, J., Gehl, P., Bonilla, L. F. & Gélis, C., 2010. A *κ* model for mainland France, *Pure appl. Geophys.*, **167**(11), 1303–1315.

Dowrick, D. & Rhoades, D., 1999. Attenuation of Modified Mercalli Intensity in New Zealand earthquakes, *Bull. N.Z. Soc. Earthq. Eng.*, **32**(2), 55–89.

Eberhart-Phillips, D. & Bannister, S., 2010. 3-D imaging of Marlborough, New Zealand, subducted plate and strike-slip fault systems, *Geophys. J. Int.*, **182**(1), 73–96.

Eberhart-Phillips, D. & Chadwick, M., 2002. Three-dimensional attenuation model of the shallow Hikurangi subduction zone in the Raukumara Peninsula, New Zealand, *J. geophys. Res.*, **107**(B2), doi:10.1029/2000JB000046.

Eberhart-Phillips, D. & McVerry, G., 2003. Estimating slab earthquake response spectra from a 3D Q model, *Bull. seism. Soc. Am.*, **93**(6), 2649–2663.

Eberhart-Phillips, D., Reyners, M., Chadwick, M. & Chiu, J.-M., 2005. Crustal heterogeneity and subduction processes: 3-D Vp, Vp/Vs and Q in the southern North Island, New Zealand, *Geophys. J. Int.*, **162**(1), 270–288.

Eberhart-Phillips, D., Chadwick, M. & Bannister, S., 2008. Three-dimensional attenuation structure of central and southern South Island, New Zealand, from local earthquakes, *J. geophys. Res.*, **113**(B5), doi:10.1029/2007JB005359.

Eberhart-Phillips, D., Reyners, M. & Bannister, S., 2015. A 3D *Q*ₚ attenuation model for all of New Zealand, *Seismol. Res. Lett.*, **86**(6), 1655–1663.

Edwards, B., Ktenidou, O.-J., Cotton, F., Abrahamson, N., Van Houtte, C. & Fäh, D., 2015. Epistemic uncertainty and limitations of the *κ*₀ model for near-surface attenuation at hard rock sites, *Geophys. J. Int.*, **202**(3), 1627–1645.

EPRI, 1993. Method and guidelines for estimating earthquake ground motion in eastern North America, Guidelines for Determining Design Basis Ground Motions, EPRI TR-102293, Vol. 1, Electric Power Research Institute, Palo Alto, CA.

Goldstein, P. & Snoke, A., 2005. SAC availability for the IRIS community, *Incorporated Institutions for Seismology Data Management Center Electronic Newsletter*, 7(1).

Graves, R. & Pitarka, A., 2010. Broadband ground-motion simulation using a hybrid approach, *Bull. seism. Soc. Am.*, **100**(5A), 2095–2123.

Haskell, N., 1964. Total energy and energy spectral density of elastic wave radiation from propagating faults, *Bull. seism. Soc. Am.*, **54**(6A), 1811–1841.

Hough, S. & Anderson, J., 1988. High-frequency spectra observed at Anza, California: implications for *Q* structure, *Bull. seism. Soc. Am.*, **78**(2), 692–707.

Hough, S. et al., 1988. Attenuation near Anza, California, *Bull. seism. Soc. Am.*, **78**(2), 672–691.

Housner, G., 1947. Characteristics of strong-motion earthquakes, *Bull. seism. Soc. Am.*, **37**(1), 19–31.

Kaneko, Y. & Shearer, P., 2014. Seismic source spectra and estimated stress drop derived from cohesive-zone models of circular subshear rupture, *Geophys. J. Int.*, **197**(2), 1002–1015.

Kaneko, Y. & Shearer, P., 2015. Variability of seismic source spectra, estimated stress drop, and radiated energy, derived from cohesive-zone models of symmetrical and asymmetrical circular and elliptical ruptures, *J. geophys. Res.*, **120**(2), 1053–1079.
Ktenidou, O.-J. & Abrahamson, N., 2016. Empirical estimation of high-frequency ground motion on hard rock, *Seismol. Res. Lett.*, **87**(6), 1465–1478.

Ktenidou, O.-J., Gélis, C. & Bonilla, L.-F., 2013. A study on the variability of kappa (κ) in a borehole: implications of the computation process, *Bull. seism. Soc. Am.*, **103**(2A), 1048–1068.

Ktenidou, O.-J., Cotton, F., Abrahamson, N. & Anderson, J., 2014. Taxonomy of κ: A review of definitions and estimation approaches targeted to applications, *Seismol. Res. Lett.*, **85**(1), 135–146.

Ktenidou, O.-J., Abrahamson, N., Drouet, S. & Cotton, F., 2015. Understanding the physics of kappa (κ): insights from a downhole array, *Geophys. J. Int.*, **203**(1), 678–691.

Ktenidou, O.-J., Abrahamson, N., Darragh, R. & Silva, W., 2016. A methodology for the estimation of kappa (κ) for large datasets: example application to rock sites in the NGA-East database, Pacific Earthquake Engineering Research Center Report 2016/01, Berkeley, CA, USA.

Ktenidou, O.-J., Silva, W., Darragh, R., Abrahamson, N. & Kishida, T., 2017. Squeezing kappa (κ) out of the Transportable array: a strategy for using bandlimited data in regions of sparse seismicity, *Bull. seism. Soc. Am.*, **107**(1), 256–275.

Laurendeau, A., Cotton, F., Ktenidou, O.-J., Bonilla, L.-F. & Hollender, F., 2013. Rock and stiff soil site amplification: dependency on $v_{S30}$ and kappa ($\kappa_0$), *Bull. seism. Soc. Am.*, **103**(6), 3131–3148.

Lay, T. & Wallace, T., 1995. *Modern Global Seismology*, Academic Press.

Leary, P. & Abercrombie, R., 1994. Frequency-dependent crustal scattering and absorption at 5–160 Hz from coda decay observed at 2.5 km depth, *Geophys. Res. Lett.*, **21**(11), 971–974.

Lermo, J. & Chávez-García, F., 1993. Site effect evaluation using spectral ratios with only one station, *Bull. seism. Soc. Am.*, **83**(5), 1574–1594.

Madariaga, R., 1976. Dynamics of an expanding circular fault, *Bull. seism. Soc. Am.*, **66**(3), 639–666.

Matérn, B., 1960. Spatial variation, Technical Report, Statens Skogsforsningsinstitut, Stockholm.

McShane, B., Gal, D., Gelman, A., Robert, C. & Tackett, J., 2017. *Abandon statistical significance*, arXiv:1709.07588.

McVerry, G., Zhao, J., Abrahamson, N. & Somerville, P., 2006. New Zealand acceleration response spectrum attenuation relations for crustal and subduction zone earthquakes, *Bull. N.Z. Soc. Earthq. Eng.*, **39**(1), 1–58.

Mereu, R., Dineva, S. & Atkinson, G., 2013. The application of velocity spectral stacking to extract information on source and path effects for small-to-moderate earthquakes in Southern Ontario with evidence for constant-width faulting, *Seismol. Res. Lett.*, **84**(5), 899–916.

Mooney, H., 1970. Upper mantle inhomogeneity beneath New Zealand: seismic evidence, *J. geophys. Res.*, **75**(2), 285–309.

Morozov, I., 2008. Geometrical attenuation, frequency dependence of Q, and the absorption band problem, *Geophys. J. Int.*, **175**(1), 239–252.

Motazedian, D. & Atkinson, G., 2005. Stochastic finite-fault modeling based on a dynamic corner frequency, *Bull. seism. Soc. Am.*, **95**(3), 995–1010.

Oth, A. & Kaiser, A. E., 2014. Stress release and source scaling of the 2010–2011 Canterbury, New Zealand earthquake sequence from spectral inversion of ground motion data, *Pure appl. Geophys.*, **171**(10), 2767–2782.

Oth, A., Bindi, D., Parolai, S. & Di Giacomo, D., 2011. Spectral analysis of K-NET and KiK-net data in Japan, part II: On attenuation characteristics, source spectra, and site response of borehole and surface stations, *Bull. seism. Soc. Am.*, **101**(2), 667–687.

Papageorgiou, A. & Aki, K., 1983. A specific barrier model for the quantitative description of inhomogeneous faulting and the prediction of strong ground motion. I. Description of the model, *Bull. seism. Soc. Am.*, **73**(3), 693–722.

Pei, S. et al., 2009. Structure of the upper crust in Japan from S-wave attenuation tomography, *Bull. seism. Soc. Am.*, **99**(1), 428–434.

Perron, V., Hollender, F., Bard, P.-Y., Gélis, C., Guyonnet-Benaize, C., Hernandez, B. & Ktenidou, O.-J., 2017. Robustness of kappa (κ) measurement in low-to-moderate seismicity areas: insight from a site-specific study in Provence, France, *Bull. seism. Soc. Am.*, **107**(5), 2272–2292.

Petersen, T., Gledhill, K., Chadwick, M., Gale, N. & Ristau, J., 2011. The New Zealand national seismograph network, *Seismol. Res. Lett.*, **82**(1), 9–20.

Ribeiro, P. & Diggle, P., 2001. geoR: a package for geostatistical analysis, *R-NEWS*, **1**(2), 15–18.

Silva, W. & Darragh, R., 1995. Engineering characterization of earthquake strong ground motion recorded at rock sites, Electric Power Research Institute, Report TR-102261, Palo Alto, California.

Standards New, Zealand, 2004. *Structural Design Actions Part 5: Earthquake Actions*, Department of Building and Housing, Wellington, New Zealand.

Steidl, J., Tumarkin, A. & Archuleta, R., 1996. What is a reference site?, *Bull. seism. Soc. Am.*, **86**(6), 1733–1748.

Van Houtte, C., Drouet, S. & Cotton, F., 2011. Analysis of the origins of κ (kappa) to compute hard rock to rock adjustment factors for GMPEs, *Bull. seism. Soc. Am.*, **101**(6), 2926–2941.

Van Houtte, C., Ktenidou, O., Larkin, T. & Holden, C., 2014. Hard-site $\kappa_0$ (kappa) calculations for Christchurch, and comparison with local ground-motion prediction models, *Bull. seism. Soc. Am.*, **104**(4), 1899–1913.

Wessel, P., Smith, W., Scharroo, R., Luis, J. & Wobbe, F., 2013. Generic mapping tools: improved version released, *Eos Trans. AGU*, **94**(45), 409–410.

Wilson, C., Houghton, B., McWilliams, M., Lanphere, M., Weaver, S. & Briggs, R., 1995. Volcanic and structural evolution of Taupo Volcanic Zone, New Zealand: a review, *J. Volcanol. Geotherm. Res.*, **68**(1), 1–28.

Yogomida, K. & Benites, R., 1995. Relation between direct wave Q and coda Q: a numerical approach, *Geophys. J. Int.*, **123**(2), 471–483.

Yoshimoto, K. & Okada, M., 2009. Frequency-dependent attenuation of S-waves in the Kanto region, Japan, *Earth, Planets Space*, **61**(9), 1067–1075.

Yoshimoto, K., Sato, H. & Ohtake, M., 1993. Frequency-dependent attenuation of P and S waves in the Kanto area, Japan, based on the coda-normalization method, *Geophys. J. Int.*, **114**(1), 165–174.

## SUPPORTING INFORMATION

Supplementary data are available at [GJIRAS](https://gjiras.org) online.

**Figure S1.** Effect of changing $\eta$ values when fitting the high-frequency spectral shape. Note that large variations in $Q_0$ and $\eta$ are possible with little effect on the quality of the fit.

**Figure S2.** Illustration of the degree of variability in the mean $\kappa$ values (red points), as shown in the $\pm 1$ standard deviation error bars.

**Figure S3.** Motivation to ignore certain events near a localized region of low $Q$. Orange points are events with ray paths that may pass predominantly through the low $Q$ zone, red points do not. The low $Q$ zone itself is not well-defined, and is a subjective interpretation from Eberhart-Phillips et al. (2008). Its depth dependence is also ignored here. To furnish $\kappa_0$ values, one could either (i) separate data into two regions, calculate two slopes and $\kappa_0$ values, then take the average $\kappa_0$, or (ii) combine the data together with a single slope and $\kappa_0$. Option (i) yields plausible results for RPZ and WVZ, but not INZ, while option (ii) gives plausible results for WVZ, but not RPZ or INZ. Hence our final decision was to ignore the red points from INZ and yellow points from RPZ, as we believe that this is the most pragmatic solution. Note that this decision is unlikely to have an adverse effect on our inference, as our preferred ‘model 2’ $\kappa_0$ map essentially smooths over the station-specific $\kappa_0$ values in the South Island.

**Figure S4.** Results of the leave-one-out cross validation. These maps show the jackknife standard deviation for Model 1 (a) and Model 2 (b).

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.