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

A continuous map of near-surface S-wave attenuation in New Zealand 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, VS30, 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. 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 (VS30) and the depth to a shear wave velocity horizon of 1 km s−1 (Z1.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 VS30 (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 VS30 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 VS30 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,$$ (1)where A0 is a source- and path-dependent Fourier amplitude of acceleration, f  is the frequency, fe 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, fe 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) + \tilde{\kappa }(r),$$ (2)where $$\tilde{\kappa }(r)$$ is the distance-dependence of κ and represents attenuation along the ray path before seismic waves enter the near-surface material, and κ0(s) is the site attenuation parameter. Papageorgiou & Aki (1983) initially suggested that κ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 κ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 fe, 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) = \displaystyle \int _{\text{path}}{\frac{\text{d}r}{Q(z)V(z)}}\ ,$$ (3)(Hough & Anderson 1988), where z is the depth and V is wave velocity. Conditioned on these assumptions, κ(r, s), hereafter denoted as κ, 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 fe may not be a good representation of reality. An ω-square source model (e.g. Housner 1947; Haskell 1964; Aki 1967; Brune 1970), for which the displacement amplitudes decay according to ω−2 for large ω, 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 ω-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 κ 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 ω-square source model applies to all events and stations used to calculate κ values. For the κ model to represent attenuation, it is also necessary to assume that site effects above fe 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 κ. 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 κ 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 κ model is the frequency-indepen-dence of Q above fe. 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 Qe depends strongly on frequency, at least between around 1 and 10 Hz. For this reason, Qe is often modelled using the parametric function   $$Q_e(f)=Q_0f^\eta ,$$ (4)where Q0 is the value of Qe at a reference frequency, such as 1 Hz, and η is an exponent controlling the degree of frequency-dependence. Atkinson (2012) and Mereu et al. (2013) compile values of η for numerous attenuation studies in eastern North America, and find that most studies have η 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 Qe 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 Qe has a weaker dependence on frequency above 10 Hz (Leary & Abercrombie 1994; Adams & Abercrombie 1998; Abercrombie 1998; Yoshimoto & Okada 2009), although η 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 κ values will be the same. Given the trade-offs between source, path and site effects, testing the sensitivity of κ values (as defined by eq. 1) to the value of η is difficult. κ values calculated using eq. (1) are sensitive only to spectral shapes and not absolute amplitudes. If a model with η > 0 is fit to the high-frequency spectrum using eq. (1), the fit to the data will not change significantly, but the obtained Q0 value will differ greatly from an η = 0 model (see Supporting Information Fig. S1 for a graphical example of this effect). However, even if a preferred value of η > 0 is propagated into the κ0 calculations, extrapolating the κ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 η = 0 for the purposes of deriving κ0. The recent results of Perron et al. (2017) suggest the practical implications of this assumption may be minimal. These authors show that calculating κ 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 κ0, for two sites in France. Despite these seemingly limiting assumptions, the κ model has become a popular method for calculating site attenuation for engineering applications. The κ 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),$$ (5)where P(f ) is the attenuation function, A0 is a reference amplitude, f is the frequency and R is the distance along the ray path. However, the site attenuation parameter κ0(s) (from eq. 2), subsequently denoted as κ0, has become a widely used parameter for modelling near-surface attenuation at rock sites. κ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 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 VS30. 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–VS30 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–VS30 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. Figure 1. View largeDownload slide 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. Figure 1. View largeDownload slide 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. 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 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, fc, 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 fc 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 fc on the calculated κ slope. The available data span a magnitude range from Mw2.5 to Mw6.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. Figure 2. View largeDownload slide Magnitude and distance distribution of earthquakes that form the data set for this study. Figure 2. View largeDownload slide Magnitude and distance distribution of earthquakes that form the data set for this study. 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 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, fe and fx. fe 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, fx is selected to always be lower than 40 Hz, although in some cases fx is limited by the SNR criterion. A minimum value of fx − fe = 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. Figure 3. View largeDownload slide 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) fc,  fe and fx picks for calculating κ. Figure 3. View largeDownload slide 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) fc,  fe and fx picks for calculating κ. 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{\kappa }(r)$$ term from eq. 2). The distance-dependence is approximately related to the intrinsic Q for S-waves (QS) and VS structure along the horizontal wave-propagation path from source to site, as shown in eq. (3). For a medium with constant QS and VS, κ 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{\kappa }(r)$$ have also been proposed to better represent the geological structure in the region of interest, this study uses a constant QS and VS model,   $$\tilde{\kappa }(r) \approx \frac{R}{Q_SV_S}\ ,$$ (6)where R is the epicentral distance. A more complex form for $$\tilde{\kappa }(r)$$ was considered, but was limited by depth uncertainties for many recorded events, making it too difficult to reliably trace the QS 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 (labelled as ‘free Q’ in Fig. 4), the QS values implied by the $$\tilde{\kappa }(r)$$ slope range from 900-1600, assuming constant VS of 3.5 km s−1. These QS values can be compared to other Q values obtained in an independent study. The 3-D South Island attenuation model of Eberhart-Phillips et al. (2008) calculates the frequency-independent Q structure for P-waves, QP, in Otago. In the typical frequency bandwidth 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 QP values ranging from 900–1100. Given that Eberhart-Phillips et al. (2008) calculate QP and this study calculates QS, it is not possible to directly compare the results. Lay & Wallace (1995) and Castro et al. (1999) suggest QP/QS is around 2 and 1 respectively, although it is not clear how these values should depend on frequency. Yoshimoto et al. (1993) show that QP/QS 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 QS values from this study compare well with the general regional trends of the QP structure derived by Eberhart-Phillips et al. (2008). Figure 4. View largeDownload slide (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. Figure 4. View largeDownload slide (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. 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{\kappa }(r)$$. In this case, the combined regression yields an ‘average regional QS’ 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: DCZ and PYZ, in Fiordland; WEL and BHW, two stations at the southern tip of the North Island; and, 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. 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 to distance and other factors (the R-squared for a linear model of κobserved − κpredicted versus ML 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. Figure 5. View largeDownload slide κ 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 κ data from multiple neighbouring stations. Calculated κ0 values are also indicated in the plots. Figure 5. View largeDownload slide κ 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 κ data from multiple neighbouring stations. Calculated κ0 values are also indicated in the plots. Figure 6. View largeDownload slide κ data for North Island stations in this study. Refer to Fig. 1 for station locations. As with Figure Dark lines indicate an unconstrained regression (free Q), while light lines indicate fixed regional Q slopes calculated using κ data from adjacent stations. Dashed lines represent 5–95 per cent confidence intervals from the unconstrained regression. Calculated κ0 values are indicated in the plots. Figure 6. View largeDownload slide κ data for North Island stations in this study. Refer to Fig. 1 for station locations. As with Figure Dark lines indicate an unconstrained regression (free Q), while light lines indicate fixed regional Q slopes calculated using κ data from adjacent stations. Dashed lines represent 5–95 per cent confidence intervals from the unconstrained regression. Calculated κ0 values are indicated in the plots. Figure 7. View largeDownload slide (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. Figure 7. View largeDownload slide (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. Table 1. Summary of κ data from New Zealand. Refer to Fig. 1 for station locations. Station  κ0 (s)  Slope dκ/dR  Number of recordings  Fiordland      3.62 × 10−4    DCZ  0.013    41      3.34 × 10−4a        3.11 × 10−4    PYZ  0.014    22      3.34 × 10−4a    Otago  WHZ  0.009  1.79 × 10−4  28      1.72 × 10−4    SYZ  0.012    18      2.51 × 10−4a        2.54 × 10−4    TUZ  0.006    50      2.51 × 10−4a        3.15 × 10−4    OPZ  0.010    24      2.51 × 10−4a    Canterbury  EAZ  0.016  3.33 × 10−4  50  JCZ  0.027  2.94 × 10−4  45  LBZ  0.022  2.14 × 10−5  37  RPZ  0.021  2.93 × 10−4  26  LTZ  0.016  3.85 × 10−4  22  KHZ  0.031  5.16 × 10−4  44  Banks Peninsulab  CRLZ  0.032  –  –  MQZ  0.030  –  –  AKSS  0.033  –  –  D14C  0.025  –  –  MTPS  0.039  –  –  West Coast  FOZ  0.019  1.63 × 10−4  30  WVZ  0.009  4.08 × 10−4  19  INZ  0.022  3.91 × 10−4  16  DSZ  0.023  1.91 × 10−4  28  THZ  0.012  4.18 × 10−4  50  NNZ  0.021  5.78 × 10−4  21  Forearc / East Coast Ranges  PAWZ  0.052  5.93 × 10−4  11  BFZ  0.050  5.27 × 10−4  11  PXZ  0.044  3.87 × 10−4  9  KNZ  0.042  9.67 × 10−4  9  CKHZ  0.053  7.35 × 10−4  15  PUZ  0.048  1.61 × 10−3  39  MWZ  0.050  2.82 × 10−4  18  MXZ  0.030  1.49 × 10−3  15  Arc Ranges      4.10 × 10−4    WEL  0.032    48      3.73 × 10−4a        4.77 × 10−4    BHW  0.034    13      3.73 × 10−4a    BKZ  0.029  3.38 × 10−4  29  URZ  0.030  3.64 × 10−4  25  Taupo Volcanic Zone  WHTZ  0.055  4.34 × 10−4  12  UTU  0.053  4.46 × 10−4  10  MARZ  0.052  6.59 × 10−4  17  Northern Districts      1.88 × 10−4    WAZ  0.019    76      2.10 × 10−4a        3.13 × 10−4    VRZ  0.017  2.10 × 10−4a  14  HIZ  0.015  2.38 × 10−4  19  TLZ  0.016  3.24 × 10−4  14  TGRZ  0.024  2.80 × 10−4  29  TOZ  0.027  1.63 × 10−4  26  WIAZ  0.030  2.20 × 10−4  13  WCZ  0.029  5.39 × 10−4  7  Station  κ0 (s)  Slope dκ/dR  Number of recordings  Fiordland      3.62 × 10−4    DCZ  0.013    41      3.34 × 10−4a        3.11 × 10−4    PYZ  0.014    22      3.34 × 10−4a    Otago  WHZ  0.009  1.79 × 10−4  28      1.72 × 10−4    SYZ  0.012    18      2.51 × 10−4a        2.54 × 10−4    TUZ  0.006    50      2.51 × 10−4a        3.15 × 10−4    OPZ  0.010    24      2.51 × 10−4a    Canterbury  EAZ  0.016  3.33 × 10−4  50  JCZ  0.027  2.94 × 10−4  45  LBZ  0.022  2.14 × 10−5  37  RPZ  0.021  2.93 × 10−4  26  LTZ  0.016  3.85 × 10−4  22  KHZ  0.031  5.16 × 10−4  44  Banks Peninsulab  CRLZ  0.032  –  –  MQZ  0.030  –  –  AKSS  0.033  –  –  D14C  0.025  –  –  MTPS  0.039  –  –  West Coast  FOZ  0.019  1.63 × 10−4  30  WVZ  0.009  4.08 × 10−4  19  INZ  0.022  3.91 × 10−4  16  DSZ  0.023  1.91 × 10−4  28  THZ  0.012  4.18 × 10−4  50  NNZ  0.021  5.78 × 10−4  21  Forearc / East Coast Ranges  PAWZ  0.052  5.93 × 10−4  11  BFZ  0.050  5.27 × 10−4  11  PXZ  0.044  3.87 × 10−4  9  KNZ  0.042  9.67 × 10−4  9  CKHZ  0.053  7.35 × 10−4  15  PUZ  0.048  1.61 × 10−3  39  MWZ  0.050  2.82 × 10−4  18  MXZ  0.030  1.49 × 10−3  15  Arc Ranges      4.10 × 10−4    WEL  0.032    48      3.73 × 10−4a        4.77 × 10−4    BHW  0.034    13      3.73 × 10−4a    BKZ  0.029  3.38 × 10−4  29  URZ  0.030  3.64 × 10−4  25  Taupo Volcanic Zone  WHTZ  0.055  4.34 × 10−4  12  UTU  0.053  4.46 × 10−4  10  MARZ  0.052  6.59 × 10−4  17  Northern Districts      1.88 × 10−4    WAZ  0.019    76      2.10 × 10−4a        3.13 × 10−4    VRZ  0.017  2.10 × 10−4a  14  HIZ  0.015  2.38 × 10−4  19  TLZ  0.016  3.24 × 10−4  14  TGRZ  0.024  2.80 × 10−4  29  TOZ  0.027  1.63 × 10−4  26  WIAZ  0.030  2.20 × 10−4  13  WCZ  0.029  5.39 × 10−4  7  aRegionally fixed path dependence. bValues calculated in Van Houtte et al. (2014). View Large 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. Figure 8. View largeDownload slide Calculated κ0 values for the stations analysed in this study, superimposed on a geological map. The length of the bars is proportional to the κ0 values, with a scale included for reference. Dashed lines indicate boundaries of regions with similar κ0 values (discussed in detail within the text), while the solid black line is the boundary between the Australian and Pacific tectonic plates. Figure 8. View largeDownload slide Calculated κ0 values for the stations analysed in this study, superimposed on a geological map. The length of the bars is proportional to the κ0 values, with a scale included for reference. Dashed lines indicate boundaries of regions with similar κ0 values (discussed in detail within the text), while the solid black line is the boundary between the Australian and Pacific tectonic plates. Figure 9. View largeDownload slide Regional differences in κ0. Each point corresponds to an individual station κ0 value, grouped into the regions indicated by Fig. 8. Figure 9. View largeDownload slide Regional differences in κ0. Each point corresponds to an individual station κ0 value, grouped into the regions indicated by Fig. 8. 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 QP. κ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 QS 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 QP 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 QS structure near the Alpine Fault, and on the West Coast. 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 κ0 values in Van Houtte et al. (2014). Five of these stations are deemed to be good representations of rock sites, and κ0 values from these stations are indicated in Figs 8 and 9 and Table 1. Observed κ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 κ0 values in the country. Sedimentary rock sites located on the coastal ranges and in the forearc basin have κ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 QP in the 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 QP 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 QS 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 QP model of Eberhart-Phillips et al. (2005) suggests that either side of the high-QP arc ranges, the forearc and volcanic fronts have very low QP 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 New Zealand κ0 data in ground motion modelling, and subsequently seismic hazard, a continuous map is developed. In this section, the discrete station-specific κ0 values are spatially smoothed and interpolated, to develop a continuous rock-site κ0 map for New Zealand. A number of spatial interpolation schemes were considered as candidates for developing the κ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 κ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 xi, that is,   $$K_{0,i} = S(x_i)+\epsilon _i\ :\ i = 1,\ldots ,n,$$ (7)where K0, i is a vector of n observed κ0 values, and εi are normally distributed errors with mean of 0 and variance of τ2. S(x), known as the ‘signal’, has mean μ, variance σ2 and a correlation function ρ(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 κ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). 6.1 Empirical spatial behaviour A useful tool for investigating the behaviour of geospatial data is the semi-variogram, γ, which describes the variance of the difference between two realizations of the signal S(x) at locations x and x΄,   $$\gamma (x,x^{\prime }) = \frac{1}{2}\text{Var}[S(x)-S(x^{\prime })]\ .$$ (8)In the case of a stationary Gaussian process, the semi-variogram can also be expressed as   $$\gamma (x,x^\prime ) = \gamma (u) = \tau ^2+\sigma ^2[(1-\rho (u)]\ .$$ (9)A typical function for ρ(u) monotonically decreases as the distance u increases, hence γ(u) usually monotonically increases with u. Semi-variograms tend to have the following features. When u = 0, the intercept of the semi-variogram is τ2, which is known as the ‘nugget variance’, or simply ‘nugget’. Also, when ρ(u) = 0 (no correlation between observations beyond a given distance), the semi-variance becomes τ2 + σ2, called the ‘sill’. The ‘range’ of the semi-variogram is the distance at which γ(u) equals the sill. Fig. 10 shows a schematic example of a semi-variogram with ρ(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. Figure 10. View largeDownload slide An example spatial correlation function (a), and its subsequent semi-variogram (b). Explanation of functions and terms is included within the text. Figure 10. View largeDownload slide An example spatial correlation function (a), and its subsequent semi-variogram (b). Explanation of functions and terms is included within the text. The solid line in Fig. 11 is the empirical semi-variogram for the New Zealand κ0 data, determined using eq. (8), substituting the observations κ0 and $$\kappa _0^\prime$$ 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,\ldots ,n.$$ (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 κ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 κ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 11. View largeDownload slide Empirical semi-variograms of the logarithm of the New Zealand κ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. Figure 11. View largeDownload slide Empirical semi-variograms of the logarithm of the New Zealand κ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. 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 κ0 data. The trend surface consists of a single indicator variable for the Taupo Volcanic Zone, FTVZ, such that the mean function is   $$\mu (x) = \beta _0 + \beta _1\cdot F_{\text{TVZ}}$$ (11)where FTVZ 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 κ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 ρ(u) that best models the empirical κ0 semi-variances. 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)=\big [2^{\theta -1}\Gamma (\theta )\big ]^{-1}\bigg (\frac{u}{\phi }\bigg )^\theta K_\theta \bigg (\frac{u}{\phi }\bigg )\ .$$ (12)Kθ denotes a modified Bessel function with non-negative order θ, ϕ is a non-negative scale parameter with units of distance and Γ 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 θ dictates the behaviour of the correlation structure where the station separation distance is small, while ϕ controls the degree of correlation between stations with large separation distance. Fig. 12 gives an example of the effect of different values of θ and ϕ on ρ(u). Note that for θ = 0.5, ρ(u) reduces to the exponential function exp (u/ϕ). Figure 12. View largeDownload slide (a) The effect of θ on the Matérn correlation function, with a fixed value of ϕ = 0.2. (b) The effect of ϕ, for a θ value of 1.5. Figure 12. View largeDownload slide (a) The effect of θ on the Matérn correlation function, with a fixed value of ϕ = 0.2. (b) The effect of ϕ, for a θ value of 1.5. The Matérn parametric form is used to model the spatial correlation of the κ0 data. Treating the transformed spatial κ0 data as Gaussian, the model can be written as   $$\log K_0 \sim N(Y\beta ,\sigma ^2B(\theta ,\phi )+\tau ^2I)\ ,$$ (13)where K0 = (κ0, 1, …, κ0, n), Y is a matrix of covariates and β the corresponding matrix of coefficients, B is an n × n correlation matrix, which in this case depends on the Matérn parameters θ and ϕ, and I is the identity matrix. The corresponding log-likelihood function is   \begin{eqnarray} \mathcal {L}(\beta ,\tau ,\sigma ,\theta ,\phi ) &=& -0.5\big \lbrace (n\log (2\pi )+\log \big |\sigma ^2B(\theta ,\phi )+\tau ^2I\big |\nonumber\\ &&+\,(\log \kappa _0-Y\beta )^T\big [\sigma ^2B(\theta ,\phi )+\tau ^2I\big ]^{-1}\nonumber\\ &&\times \,(\log \kappa _0-Y\beta )\big \rbrace . \end{eqnarray} (14)Maximizing eq. (14) gives the model parameters. However, the Matérn parameters θ and ϕ tend to be strongly correlated, so a common, alternative method for parameter estimation is to fix a discrete set of θ, then maximize eq. (14) to determine β, τ2, σ2 and ϕ (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 κ0 uncertainty. Guided by the κ(r) regression results in Figs 5 and 6, an approximate standard deviation value of 50 per cent (0.18 log10 units) was assumed as the average uncertainty of κ0, and hence τ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 θ = 0.5 is considered, as models with higher values of θ are overparametrized and give unstable solutions. 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. Model 1 - trend surface for TVZ and free τ2  Order  β0  β1  σ2  τ2  ϕ    LL  AIC  θ = 0.5  −1.630  0.355  0.055  0.003  274.2    21.81  −33.62  θ = 1.5  −1.624  0.350  0.045  0.008  105.5    21.80  −33.60  θ = 2.5  −1.622  0.344  0.044  0.009  76.5    21.82  −33.64  Model 1 - trend surface for TVZ and free τ2  Order  β0  β1  σ2  τ2  ϕ    LL  AIC  θ = 0.5  −1.630  0.355  0.055  0.003  274.2    21.81  −33.62  θ = 1.5  −1.624  0.350  0.045  0.008  105.5    21.80  −33.60  θ = 2.5  −1.622  0.344  0.044  0.009  76.5    21.82  −33.64  Model 2 - trend surface for TVZ and fixed τ2 = 0.03  Order  β0  β1  σ2    ϕ    LL  AIC  θ = 0.5  −1.654  0.294  0.045    646.0    14.31  −20.62  Model 2 - trend surface for TVZ and fixed τ2 = 0.03  Order  β0  β1  σ2    ϕ    LL  AIC  θ = 0.5  −1.654  0.294  0.045    646.0    14.31  −20.62  View Large The first point from Table 2 is that ϕ usually decreases as θ increases, illustrating the correlation between the parameters. However, the most important effect of the θ parameter is to increase the nugget variance τ2. The models with different θ values have similar log-likelihoods, which indicates that τ2 is not well-constrained by the data. The models with free τ2 have higher log-likelihoods than the model with fixed τ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 θ values at short station-separation distances. The more complicated form of the θ = 1.5 and θ = 2.5 models is not well justified by the data, hence the model with θ = 0.5 is preferred in this study. Figure 13. View largeDownload slide The binned empirical variogram for the logarithm of the New Zealand κ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. Figure 13. View largeDownload slide The binned empirical variogram for the logarithm of the New Zealand κ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. 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 κ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}\displaystyle \sum _{i=1}^n[(\hat{S}(x)-S(x_i))]^2 \ ,$$ (15)where $$\hat{S}(x)$$ is the prediction of S(x) at an unobserved location. The form of $$\hat{S}(x)$$ that minimizes $$\text{MSPE}[\hat{S}(x)]$$ is,   $$\hat{S}(x)=\mu +b^\prime V^{-1}(K_0-\mu ),$$ (16)where μ 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 + (τ2/σ2)I. In ordinary kriging, μ is assumed to be unknown and is estimated by generalized least squares, given by   $$\hat{\mu }=(\mathbf {1^\prime }V^{-1}\mathbf {1})^{-1}\mathbf {1^\prime }V^{-1}K_0.$$ (17)The prediction variance is given by   $$\text{Var}(S(x)|K_0)=\sigma (1-b^\prime V^{-1}b),$$ (18)Full derivations of these expressions are readily available in the literature. As the models contain the FTVZ 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 κ0 for models 1 and 2, both with θ equal to 0.5. Due to the larger τ2, model 2 predicts a smoother median surface. The effect of increasing τ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 τ2. This calculation means that the maps in Figs 14(b) and (d) provide the prediction standard deviation for an unobserved κ0 value. Figure 14. View largeDownload slide (a) Median κ0 map for Model 1 with θ = 0.5 and (b) its standard deviation in log10(s). (c) The median model for Model 2 with θ = 0.5 and (d) its standard deviation. Figure 14. View largeDownload slide (a) Median κ0 map for Model 1 with θ = 0.5 and (b) its standard deviation in log10(s). (c) The median model for Model 2 with θ = 0.5 and (d) its standard deviation. 6.4 Consideration of κ0 ‘measurement uncertainty’ Fig. 14 clearly demonstrates that the nugget is an important parameter when considering forward application of the κ0 map. While the seismograph network in New Zealand is insufficiently dense to derive the nugget parameter from the kriging analysis alone, estimates of κ0 measurement uncertainty are available from the linear regression of the κ(r) data with distance (Figs 5 and 6). This is the purpose of model 2 in Table 2. The fixed τ2 value of 0.03 for model 5 is only an indicative value of the κ0 uncertainty across all stations. A limitation of this assumption is that some of the κ0 data are better constrained than others, as evidenced by the confidence intervals in Figs 5 and 6. Ideally, a fixed vector of τ2 with each station’s κ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 κ0 map. The station-specific κ0 are derived by regression of κ data with distance on a linear scale, while a lognormal distribution was assumed to derive the κ map. The combination of distributions prevents an elegant determination of the spatially varying nugget, hence the fixed nugget of 0.03 (± 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 κ0 to constrain the nugget, which are mostly unavailable in New Zealand. For model 1, the predicted standard deviation on κ0 is around 0.1 log10 units, while for model 2 it is around 0.2 log10 units. Note that this 0.2 is only slightly larger than the value that was initially assumed for the uncertainty in the station-specific κ0 values, which suggests that the uncertainty in calculating κ0 at a given station location is much larger than the uncertainty introduced by kriging. 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 QP 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 physical 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. https://doi.org/10.1007/s000240050204 Google Scholar CrossRef Search ADS   Abrahamson N., Silva W., 1997. Empirical response spectral attenuation relations for shallow crustal earthquakes, Seismol. Res. Lett. , 68( 1), 94– 127. https://doi.org/10.1785/gssrl.68.1.94 Google Scholar CrossRef Search ADS   Abrahamson N., Silva W., Kamai R., 2014. Summary of the ASK14 ground motion relation for active crustal regions, Earthq. Spectra , 30( 3), 1025– 1055. https://doi.org/10.1193/070913EQS198M Google Scholar CrossRef Search ADS   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. https://doi.org/10.1029/98JB01757 Aki K., 1967. Scaling law of seismic spectrum, J. geophys. Res. , 72( 4), 1217– 1231. https://doi.org/10.1029/JZ072i004p01217 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1016/0031-9201(80)90019-9 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1785/gssrl.83.1.166 Google Scholar CrossRef Search ADS   Atkinson G., Boore D., 2006. Earthquake ground-motion prediction equations for eastern North America, Bull. seism. Soc. Am. , 96( 6), 2181– 2205. https://doi.org/10.1785/0120050245 Google Scholar CrossRef Search ADS   Bannister S., Gledhill K., 2012. Evolution of the 2010–2012 Canterbury earthquake sequence, N.Z. J. Geol. Geophys. , 55( 3), 295– 304. https://doi.org/10.1080/00288306.2012.680475 Google Scholar CrossRef Search ADS   Benjamin D. et al.  , 2018. Redefine statistical significance, Nat. Hum. Behav. , 2( 1), 6– 10. Google Scholar CrossRef Search ADS   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. https://doi.org/10.1007/PL00012553 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1193/070113EQS184M Google Scholar CrossRef Search ADS   Brune J., 1970. Tectonic stress and the spectra of seismic shear waves from earthquakes, J. geophys. Res. , 75( 26), 4997– 5009. https://doi.org/10.1029/JB075i026p04997 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1785/0120020002 Google Scholar CrossRef Search ADS   Campbell K., 2009. Estimates of shear-wave Q and κ0 for unconsolidated and semiconsolidated sediments in eastern North America, Bull. seism. Soc. Am. , 99( 4), 2365– 2392. https://doi.org/10.1785/0120080116 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1193/062913EQS175M Google Scholar CrossRef Search ADS   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. https://doi.org/10.1016/S0040-1951(98)00277-7 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1016/j.soildyn.2006.02.010 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1193/072813EQS219M Google Scholar CrossRef Search ADS   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. https://doi.org/10.1007/s10950-005-9006-7 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1029/GL008i011p01126 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1142/S1363246906002566 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1785/0120080075 Google Scholar CrossRef Search ADS   Douglas J., Gehl P., Bonilla L. F., Gélis C., 2010. A κ model for mainland France, Pure appl. Geophys. , 167( 11), 1303– 1315. https://doi.org/10.1007/s00024-010-0146-5 Google Scholar CrossRef Search ADS   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. https://doi.org/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. https://doi.org/10.1785/0120030036 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1111/j.1365-246X.2005.02530.x Google Scholar CrossRef Search ADS   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. https://doi.org/10.1029/2007JB005359 Eberhart-Phillips D., Reyners M., Bannister S., 2015. A 3D QP attenuation model for all of New Zealand, Seismol. Res. Lett. , 86( 6), 1655– 1663. https://doi.org/10.1785/0220150124 Google Scholar CrossRef Search ADS   Edwards B., Ktenidou O.-J., Cotton F., Abrahamson N., Van Houtte C., Fäh D., 2015. Epistemic uncertainty and limitations of the κ0 model for near-surface attenuation at hard rock sites, Geophys. J. Int. , 202( 3), 1627– 1645. https://doi.org/10.1093/gji/ggv222 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1785/0120100057 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1093/gji/ggu030 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1002/2014JB011642 Google Scholar CrossRef Search ADS   Ktenidou O.-J., Abrahamson N., 2016. Empirical estimation of high-frequency ground motion on hard rock, Seismol. Res. Lett. , 87( 6), 1465– 1478. https://doi.org/10.1785/0220160075 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1785/0120120093 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1785/0220130027 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1093/gji/ggv315 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1785/0120150301 Google Scholar CrossRef Search ADS   Laurendeau A., Cotton F., Ktenidou O.-J., Bonilla L.-F., Hollender F., 2013. Rock and stiff soil site amplification: dependency on vS30 and kappa (κ0), Bull. seism. Soc. Am. , 103( 6), 3131– 3148. https://doi.org/10.1785/0120130020 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1029/94GL00977 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1785/0220130009 Google Scholar CrossRef Search ADS   Mooney H., 1970. Upper mantle inhomogeneity beneath New Zealand: seismic evidence, J. geophys. Res. , 75( 2), 285– 309. https://doi.org/10.1029/JB075i002p00285 Google Scholar CrossRef Search ADS   Morozov I., 2008. Geometrical attenuation, frequency dependence of Q, and the absorption band problem, Geophys. J. Int. , 175( 1), 239– 252. https://doi.org/10.1111/j.1365-246X.2008.03888.x Google Scholar CrossRef Search ADS   Motazedian D., Atkinson G., 2005. Stochastic finite-fault modeling based on a dynamic corner frequency, Bull. seism. Soc. Am. , 95( 3), 995– 1010. https://doi.org/10.1785/0120030207 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1007/s00024-013-0751-1 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1785/0120100135 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1785/0120080029 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1785/0120160374 Google Scholar CrossRef Search ADS   Petersen T., Gledhill K., Chadwick M., Gale N., Ristau J., 2011. The New Zealand national seismograph network, Seismol. Res. Lett. , 82( 1), 9– 20. https://doi.org/10.1785/gssrl.82.1.9 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1785/0120100345 Google Scholar CrossRef Search ADS   Van Houtte C., Ktenidou O., Larkin T., Holden C., 2014. Hard-site κ0 (kappa) calculations for Christchurch, and comparison with local ground-motion prediction models, Bull. seism. Soc. Am. , 104( 4), 1899– 1913. https://doi.org/10.1785/0120130271 Google Scholar CrossRef Search ADS   Wessel P., Smith W., Scharroo R., Luis J., Wobbe F., 2013. Generic mapping tools: improved version released, Eos Trans. AGU , 94( 45), 409– 410. https://doi.org/10.1785/0120130271 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1016/0377-0273(95)00006-G Google Scholar CrossRef Search ADS   Yogomida K., Benites R., 1995. Relation between direct wave Q and coda Q: a numerical approach, Geophys. J. Int. , 123( 2), 471– 483. https://doi.org/10.1111/j.1365-246X.1995.tb06866.x Google Scholar CrossRef Search ADS   Yoshimoto K., Okada M., 2009. Frequency-dependent attenuation of S-waves in the Kanto region, Japan, Earth, Planets Space , 61( 9), 1067– 1075. https://doi.org/10.1186/BF03352957 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1111/j.1365-246X.1993.tb01476.x Google Scholar CrossRef Search ADS   SUPPORTING INFORMATION Supplementary data are available at GJIRAS online. Figure S1. Effect of changing η values when fitting the high-frequency spectral shape. Note that large variations in Q0 and η are possible with little effect on the quality of the fit. Figure S2. Illustration of the degree of variability in the mean κ values (red points), as shown in the ±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 κ0 values, one could either (i) separate data into two regions, calculate two slopes and κ0 values, then take the average κ0, or (ii) combine the data together with a single slope and κ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’ κ0 map essentially smooths over the station-specific κ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. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

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

, Volume 213 (1) – Apr 1, 2018
18 pages

Loading next page...

/lp/ou_press/a-continuous-map-of-near-surface-s-wave-attenuation-in-new-zealand-JWqovDb0jG
Publisher
Oxford University Press
Copyright
© The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx559
Publisher site
See Article on Publisher Site

### Abstract

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, VS30, 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. 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 (VS30) and the depth to a shear wave velocity horizon of 1 km s−1 (Z1.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 VS30 (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 VS30 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 VS30 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,$$ (1)where A0 is a source- and path-dependent Fourier amplitude of acceleration, f  is the frequency, fe 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, fe 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) + \tilde{\kappa }(r),$$ (2)where $$\tilde{\kappa }(r)$$ is the distance-dependence of κ and represents attenuation along the ray path before seismic waves enter the near-surface material, and κ0(s) is the site attenuation parameter. Papageorgiou & Aki (1983) initially suggested that κ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 κ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 fe, 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) = \displaystyle \int _{\text{path}}{\frac{\text{d}r}{Q(z)V(z)}}\ ,$$ (3)(Hough & Anderson 1988), where z is the depth and V is wave velocity. Conditioned on these assumptions, κ(r, s), hereafter denoted as κ, 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 fe may not be a good representation of reality. An ω-square source model (e.g. Housner 1947; Haskell 1964; Aki 1967; Brune 1970), for which the displacement amplitudes decay according to ω−2 for large ω, 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 ω-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 κ 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 ω-square source model applies to all events and stations used to calculate κ values. For the κ model to represent attenuation, it is also necessary to assume that site effects above fe 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 κ. 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 κ 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 κ model is the frequency-indepen-dence of Q above fe. 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 Qe depends strongly on frequency, at least between around 1 and 10 Hz. For this reason, Qe is often modelled using the parametric function   $$Q_e(f)=Q_0f^\eta ,$$ (4)where Q0 is the value of Qe at a reference frequency, such as 1 Hz, and η is an exponent controlling the degree of frequency-dependence. Atkinson (2012) and Mereu et al. (2013) compile values of η for numerous attenuation studies in eastern North America, and find that most studies have η 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 Qe 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 Qe has a weaker dependence on frequency above 10 Hz (Leary & Abercrombie 1994; Adams & Abercrombie 1998; Abercrombie 1998; Yoshimoto & Okada 2009), although η 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 κ values will be the same. Given the trade-offs between source, path and site effects, testing the sensitivity of κ values (as defined by eq. 1) to the value of η is difficult. κ values calculated using eq. (1) are sensitive only to spectral shapes and not absolute amplitudes. If a model with η > 0 is fit to the high-frequency spectrum using eq. (1), the fit to the data will not change significantly, but the obtained Q0 value will differ greatly from an η = 0 model (see Supporting Information Fig. S1 for a graphical example of this effect). However, even if a preferred value of η > 0 is propagated into the κ0 calculations, extrapolating the κ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 η = 0 for the purposes of deriving κ0. The recent results of Perron et al. (2017) suggest the practical implications of this assumption may be minimal. These authors show that calculating κ 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 κ0, for two sites in France. Despite these seemingly limiting assumptions, the κ model has become a popular method for calculating site attenuation for engineering applications. The κ 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),$$ (5)where P(f ) is the attenuation function, A0 is a reference amplitude, f is the frequency and R is the distance along the ray path. However, the site attenuation parameter κ0(s) (from eq. 2), subsequently denoted as κ0, has become a widely used parameter for modelling near-surface attenuation at rock sites. κ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 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 VS30. 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–VS30 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–VS30 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. Figure 1. View largeDownload slide 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. Figure 1. View largeDownload slide 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. 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 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, fc, 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 fc 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 fc on the calculated κ slope. The available data span a magnitude range from Mw2.5 to Mw6.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. Figure 2. View largeDownload slide Magnitude and distance distribution of earthquakes that form the data set for this study. Figure 2. View largeDownload slide Magnitude and distance distribution of earthquakes that form the data set for this study. 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 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, fe and fx. fe 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, fx is selected to always be lower than 40 Hz, although in some cases fx is limited by the SNR criterion. A minimum value of fx − fe = 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. Figure 3. View largeDownload slide 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) fc,  fe and fx picks for calculating κ. Figure 3. View largeDownload slide 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) fc,  fe and fx picks for calculating κ. 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{\kappa }(r)$$ term from eq. 2). The distance-dependence is approximately related to the intrinsic Q for S-waves (QS) and VS structure along the horizontal wave-propagation path from source to site, as shown in eq. (3). For a medium with constant QS and VS, κ 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{\kappa }(r)$$ have also been proposed to better represent the geological structure in the region of interest, this study uses a constant QS and VS model,   $$\tilde{\kappa }(r) \approx \frac{R}{Q_SV_S}\ ,$$ (6)where R is the epicentral distance. A more complex form for $$\tilde{\kappa }(r)$$ was considered, but was limited by depth uncertainties for many recorded events, making it too difficult to reliably trace the QS 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 (labelled as ‘free Q’ in Fig. 4), the QS values implied by the $$\tilde{\kappa }(r)$$ slope range from 900-1600, assuming constant VS of 3.5 km s−1. These QS values can be compared to other Q values obtained in an independent study. The 3-D South Island attenuation model of Eberhart-Phillips et al. (2008) calculates the frequency-independent Q structure for P-waves, QP, in Otago. In the typical frequency bandwidth 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 QP values ranging from 900–1100. Given that Eberhart-Phillips et al. (2008) calculate QP and this study calculates QS, it is not possible to directly compare the results. Lay & Wallace (1995) and Castro et al. (1999) suggest QP/QS is around 2 and 1 respectively, although it is not clear how these values should depend on frequency. Yoshimoto et al. (1993) show that QP/QS 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 QS values from this study compare well with the general regional trends of the QP structure derived by Eberhart-Phillips et al. (2008). Figure 4. View largeDownload slide (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. Figure 4. View largeDownload slide (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. 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{\kappa }(r)$$. In this case, the combined regression yields an ‘average regional QS’ 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: DCZ and PYZ, in Fiordland; WEL and BHW, two stations at the southern tip of the North Island; and, 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. 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 to distance and other factors (the R-squared for a linear model of κobserved − κpredicted versus ML 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. Figure 5. View largeDownload slide κ 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 κ data from multiple neighbouring stations. Calculated κ0 values are also indicated in the plots. Figure 5. View largeDownload slide κ 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 κ data from multiple neighbouring stations. Calculated κ0 values are also indicated in the plots. Figure 6. View largeDownload slide κ data for North Island stations in this study. Refer to Fig. 1 for station locations. As with Figure Dark lines indicate an unconstrained regression (free Q), while light lines indicate fixed regional Q slopes calculated using κ data from adjacent stations. Dashed lines represent 5–95 per cent confidence intervals from the unconstrained regression. Calculated κ0 values are indicated in the plots. Figure 6. View largeDownload slide κ data for North Island stations in this study. Refer to Fig. 1 for station locations. As with Figure Dark lines indicate an unconstrained regression (free Q), while light lines indicate fixed regional Q slopes calculated using κ data from adjacent stations. Dashed lines represent 5–95 per cent confidence intervals from the unconstrained regression. Calculated κ0 values are indicated in the plots. Figure 7. View largeDownload slide (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. Figure 7. View largeDownload slide (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. Table 1. Summary of κ data from New Zealand. Refer to Fig. 1 for station locations. Station  κ0 (s)  Slope dκ/dR  Number of recordings  Fiordland      3.62 × 10−4    DCZ  0.013    41      3.34 × 10−4a        3.11 × 10−4    PYZ  0.014    22      3.34 × 10−4a    Otago  WHZ  0.009  1.79 × 10−4  28      1.72 × 10−4    SYZ  0.012    18      2.51 × 10−4a        2.54 × 10−4    TUZ  0.006    50      2.51 × 10−4a        3.15 × 10−4    OPZ  0.010    24      2.51 × 10−4a    Canterbury  EAZ  0.016  3.33 × 10−4  50  JCZ  0.027  2.94 × 10−4  45  LBZ  0.022  2.14 × 10−5  37  RPZ  0.021  2.93 × 10−4  26  LTZ  0.016  3.85 × 10−4  22  KHZ  0.031  5.16 × 10−4  44  Banks Peninsulab  CRLZ  0.032  –  –  MQZ  0.030  –  –  AKSS  0.033  –  –  D14C  0.025  –  –  MTPS  0.039  –  –  West Coast  FOZ  0.019  1.63 × 10−4  30  WVZ  0.009  4.08 × 10−4  19  INZ  0.022  3.91 × 10−4  16  DSZ  0.023  1.91 × 10−4  28  THZ  0.012  4.18 × 10−4  50  NNZ  0.021  5.78 × 10−4  21  Forearc / East Coast Ranges  PAWZ  0.052  5.93 × 10−4  11  BFZ  0.050  5.27 × 10−4  11  PXZ  0.044  3.87 × 10−4  9  KNZ  0.042  9.67 × 10−4  9  CKHZ  0.053  7.35 × 10−4  15  PUZ  0.048  1.61 × 10−3  39  MWZ  0.050  2.82 × 10−4  18  MXZ  0.030  1.49 × 10−3  15  Arc Ranges      4.10 × 10−4    WEL  0.032    48      3.73 × 10−4a        4.77 × 10−4    BHW  0.034    13      3.73 × 10−4a    BKZ  0.029  3.38 × 10−4  29  URZ  0.030  3.64 × 10−4  25  Taupo Volcanic Zone  WHTZ  0.055  4.34 × 10−4  12  UTU  0.053  4.46 × 10−4  10  MARZ  0.052  6.59 × 10−4  17  Northern Districts      1.88 × 10−4    WAZ  0.019    76      2.10 × 10−4a        3.13 × 10−4    VRZ  0.017  2.10 × 10−4a  14  HIZ  0.015  2.38 × 10−4  19  TLZ  0.016  3.24 × 10−4  14  TGRZ  0.024  2.80 × 10−4  29  TOZ  0.027  1.63 × 10−4  26  WIAZ  0.030  2.20 × 10−4  13  WCZ  0.029  5.39 × 10−4  7  Station  κ0 (s)  Slope dκ/dR  Number of recordings  Fiordland      3.62 × 10−4    DCZ  0.013    41      3.34 × 10−4a        3.11 × 10−4    PYZ  0.014    22      3.34 × 10−4a    Otago  WHZ  0.009  1.79 × 10−4  28      1.72 × 10−4    SYZ  0.012    18      2.51 × 10−4a        2.54 × 10−4    TUZ  0.006    50      2.51 × 10−4a        3.15 × 10−4    OPZ  0.010    24      2.51 × 10−4a    Canterbury  EAZ  0.016  3.33 × 10−4  50  JCZ  0.027  2.94 × 10−4  45  LBZ  0.022  2.14 × 10−5  37  RPZ  0.021  2.93 × 10−4  26  LTZ  0.016  3.85 × 10−4  22  KHZ  0.031  5.16 × 10−4  44  Banks Peninsulab  CRLZ  0.032  –  –  MQZ  0.030  –  –  AKSS  0.033  –  –  D14C  0.025  –  –  MTPS  0.039  –  –  West Coast  FOZ  0.019  1.63 × 10−4  30  WVZ  0.009  4.08 × 10−4  19  INZ  0.022  3.91 × 10−4  16  DSZ  0.023  1.91 × 10−4  28  THZ  0.012  4.18 × 10−4  50  NNZ  0.021  5.78 × 10−4  21  Forearc / East Coast Ranges  PAWZ  0.052  5.93 × 10−4  11  BFZ  0.050  5.27 × 10−4  11  PXZ  0.044  3.87 × 10−4  9  KNZ  0.042  9.67 × 10−4  9  CKHZ  0.053  7.35 × 10−4  15  PUZ  0.048  1.61 × 10−3  39  MWZ  0.050  2.82 × 10−4  18  MXZ  0.030  1.49 × 10−3  15  Arc Ranges      4.10 × 10−4    WEL  0.032    48      3.73 × 10−4a        4.77 × 10−4    BHW  0.034    13      3.73 × 10−4a    BKZ  0.029  3.38 × 10−4  29  URZ  0.030  3.64 × 10−4  25  Taupo Volcanic Zone  WHTZ  0.055  4.34 × 10−4  12  UTU  0.053  4.46 × 10−4  10  MARZ  0.052  6.59 × 10−4  17  Northern Districts      1.88 × 10−4    WAZ  0.019    76      2.10 × 10−4a        3.13 × 10−4    VRZ  0.017  2.10 × 10−4a  14  HIZ  0.015  2.38 × 10−4  19  TLZ  0.016  3.24 × 10−4  14  TGRZ  0.024  2.80 × 10−4  29  TOZ  0.027  1.63 × 10−4  26  WIAZ  0.030  2.20 × 10−4  13  WCZ  0.029  5.39 × 10−4  7  aRegionally fixed path dependence. bValues calculated in Van Houtte et al. (2014). View Large 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. Figure 8. View largeDownload slide Calculated κ0 values for the stations analysed in this study, superimposed on a geological map. The length of the bars is proportional to the κ0 values, with a scale included for reference. Dashed lines indicate boundaries of regions with similar κ0 values (discussed in detail within the text), while the solid black line is the boundary between the Australian and Pacific tectonic plates. Figure 8. View largeDownload slide Calculated κ0 values for the stations analysed in this study, superimposed on a geological map. The length of the bars is proportional to the κ0 values, with a scale included for reference. Dashed lines indicate boundaries of regions with similar κ0 values (discussed in detail within the text), while the solid black line is the boundary between the Australian and Pacific tectonic plates. Figure 9. View largeDownload slide Regional differences in κ0. Each point corresponds to an individual station κ0 value, grouped into the regions indicated by Fig. 8. Figure 9. View largeDownload slide Regional differences in κ0. Each point corresponds to an individual station κ0 value, grouped into the regions indicated by Fig. 8. 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 QP. κ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 QS 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 QP 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 QS structure near the Alpine Fault, and on the West Coast. 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 κ0 values in Van Houtte et al. (2014). Five of these stations are deemed to be good representations of rock sites, and κ0 values from these stations are indicated in Figs 8 and 9 and Table 1. Observed κ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 κ0 values in the country. Sedimentary rock sites located on the coastal ranges and in the forearc basin have κ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 QP in the 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 QP 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 QS 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 QP model of Eberhart-Phillips et al. (2005) suggests that either side of the high-QP arc ranges, the forearc and volcanic fronts have very low QP 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 New Zealand κ0 data in ground motion modelling, and subsequently seismic hazard, a continuous map is developed. In this section, the discrete station-specific κ0 values are spatially smoothed and interpolated, to develop a continuous rock-site κ0 map for New Zealand. A number of spatial interpolation schemes were considered as candidates for developing the κ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 κ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 xi, that is,   $$K_{0,i} = S(x_i)+\epsilon _i\ :\ i = 1,\ldots ,n,$$ (7)where K0, i is a vector of n observed κ0 values, and εi are normally distributed errors with mean of 0 and variance of τ2. S(x), known as the ‘signal’, has mean μ, variance σ2 and a correlation function ρ(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 κ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). 6.1 Empirical spatial behaviour A useful tool for investigating the behaviour of geospatial data is the semi-variogram, γ, which describes the variance of the difference between two realizations of the signal S(x) at locations x and x΄,   $$\gamma (x,x^{\prime }) = \frac{1}{2}\text{Var}[S(x)-S(x^{\prime })]\ .$$ (8)In the case of a stationary Gaussian process, the semi-variogram can also be expressed as   $$\gamma (x,x^\prime ) = \gamma (u) = \tau ^2+\sigma ^2[(1-\rho (u)]\ .$$ (9)A typical function for ρ(u) monotonically decreases as the distance u increases, hence γ(u) usually monotonically increases with u. Semi-variograms tend to have the following features. When u = 0, the intercept of the semi-variogram is τ2, which is known as the ‘nugget variance’, or simply ‘nugget’. Also, when ρ(u) = 0 (no correlation between observations beyond a given distance), the semi-variance becomes τ2 + σ2, called the ‘sill’. The ‘range’ of the semi-variogram is the distance at which γ(u) equals the sill. Fig. 10 shows a schematic example of a semi-variogram with ρ(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. Figure 10. View largeDownload slide An example spatial correlation function (a), and its subsequent semi-variogram (b). Explanation of functions and terms is included within the text. Figure 10. View largeDownload slide An example spatial correlation function (a), and its subsequent semi-variogram (b). Explanation of functions and terms is included within the text. The solid line in Fig. 11 is the empirical semi-variogram for the New Zealand κ0 data, determined using eq. (8), substituting the observations κ0 and $$\kappa _0^\prime$$ 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,\ldots ,n.$$ (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 κ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 κ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 11. View largeDownload slide Empirical semi-variograms of the logarithm of the New Zealand κ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. Figure 11. View largeDownload slide Empirical semi-variograms of the logarithm of the New Zealand κ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. 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 κ0 data. The trend surface consists of a single indicator variable for the Taupo Volcanic Zone, FTVZ, such that the mean function is   $$\mu (x) = \beta _0 + \beta _1\cdot F_{\text{TVZ}}$$ (11)where FTVZ 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 κ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 ρ(u) that best models the empirical κ0 semi-variances. 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)=\big [2^{\theta -1}\Gamma (\theta )\big ]^{-1}\bigg (\frac{u}{\phi }\bigg )^\theta K_\theta \bigg (\frac{u}{\phi }\bigg )\ .$$ (12)Kθ denotes a modified Bessel function with non-negative order θ, ϕ is a non-negative scale parameter with units of distance and Γ 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 θ dictates the behaviour of the correlation structure where the station separation distance is small, while ϕ controls the degree of correlation between stations with large separation distance. Fig. 12 gives an example of the effect of different values of θ and ϕ on ρ(u). Note that for θ = 0.5, ρ(u) reduces to the exponential function exp (u/ϕ). Figure 12. View largeDownload slide (a) The effect of θ on the Matérn correlation function, with a fixed value of ϕ = 0.2. (b) The effect of ϕ, for a θ value of 1.5. Figure 12. View largeDownload slide (a) The effect of θ on the Matérn correlation function, with a fixed value of ϕ = 0.2. (b) The effect of ϕ, for a θ value of 1.5. The Matérn parametric form is used to model the spatial correlation of the κ0 data. Treating the transformed spatial κ0 data as Gaussian, the model can be written as   $$\log K_0 \sim N(Y\beta ,\sigma ^2B(\theta ,\phi )+\tau ^2I)\ ,$$ (13)where K0 = (κ0, 1, …, κ0, n), Y is a matrix of covariates and β the corresponding matrix of coefficients, B is an n × n correlation matrix, which in this case depends on the Matérn parameters θ and ϕ, and I is the identity matrix. The corresponding log-likelihood function is   \begin{eqnarray} \mathcal {L}(\beta ,\tau ,\sigma ,\theta ,\phi ) &=& -0.5\big \lbrace (n\log (2\pi )+\log \big |\sigma ^2B(\theta ,\phi )+\tau ^2I\big |\nonumber\\ &&+\,(\log \kappa _0-Y\beta )^T\big [\sigma ^2B(\theta ,\phi )+\tau ^2I\big ]^{-1}\nonumber\\ &&\times \,(\log \kappa _0-Y\beta )\big \rbrace . \end{eqnarray} (14)Maximizing eq. (14) gives the model parameters. However, the Matérn parameters θ and ϕ tend to be strongly correlated, so a common, alternative method for parameter estimation is to fix a discrete set of θ, then maximize eq. (14) to determine β, τ2, σ2 and ϕ (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 κ0 uncertainty. Guided by the κ(r) regression results in Figs 5 and 6, an approximate standard deviation value of 50 per cent (0.18 log10 units) was assumed as the average uncertainty of κ0, and hence τ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 θ = 0.5 is considered, as models with higher values of θ are overparametrized and give unstable solutions. 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. Model 1 - trend surface for TVZ and free τ2  Order  β0  β1  σ2  τ2  ϕ    LL  AIC  θ = 0.5  −1.630  0.355  0.055  0.003  274.2    21.81  −33.62  θ = 1.5  −1.624  0.350  0.045  0.008  105.5    21.80  −33.60  θ = 2.5  −1.622  0.344  0.044  0.009  76.5    21.82  −33.64  Model 1 - trend surface for TVZ and free τ2  Order  β0  β1  σ2  τ2  ϕ    LL  AIC  θ = 0.5  −1.630  0.355  0.055  0.003  274.2    21.81  −33.62  θ = 1.5  −1.624  0.350  0.045  0.008  105.5    21.80  −33.60  θ = 2.5  −1.622  0.344  0.044  0.009  76.5    21.82  −33.64  Model 2 - trend surface for TVZ and fixed τ2 = 0.03  Order  β0  β1  σ2    ϕ    LL  AIC  θ = 0.5  −1.654  0.294  0.045    646.0    14.31  −20.62  Model 2 - trend surface for TVZ and fixed τ2 = 0.03  Order  β0  β1  σ2    ϕ    LL  AIC  θ = 0.5  −1.654  0.294  0.045    646.0    14.31  −20.62  View Large The first point from Table 2 is that ϕ usually decreases as θ increases, illustrating the correlation between the parameters. However, the most important effect of the θ parameter is to increase the nugget variance τ2. The models with different θ values have similar log-likelihoods, which indicates that τ2 is not well-constrained by the data. The models with free τ2 have higher log-likelihoods than the model with fixed τ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 θ values at short station-separation distances. The more complicated form of the θ = 1.5 and θ = 2.5 models is not well justified by the data, hence the model with θ = 0.5 is preferred in this study. Figure 13. View largeDownload slide The binned empirical variogram for the logarithm of the New Zealand κ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. Figure 13. View largeDownload slide The binned empirical variogram for the logarithm of the New Zealand κ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. 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 κ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}\displaystyle \sum _{i=1}^n[(\hat{S}(x)-S(x_i))]^2 \ ,$$ (15)where $$\hat{S}(x)$$ is the prediction of S(x) at an unobserved location. The form of $$\hat{S}(x)$$ that minimizes $$\text{MSPE}[\hat{S}(x)]$$ is,   $$\hat{S}(x)=\mu +b^\prime V^{-1}(K_0-\mu ),$$ (16)where μ 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 + (τ2/σ2)I. In ordinary kriging, μ is assumed to be unknown and is estimated by generalized least squares, given by   $$\hat{\mu }=(\mathbf {1^\prime }V^{-1}\mathbf {1})^{-1}\mathbf {1^\prime }V^{-1}K_0.$$ (17)The prediction variance is given by   $$\text{Var}(S(x)|K_0)=\sigma (1-b^\prime V^{-1}b),$$ (18)Full derivations of these expressions are readily available in the literature. As the models contain the FTVZ 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 κ0 for models 1 and 2, both with θ equal to 0.5. Due to the larger τ2, model 2 predicts a smoother median surface. The effect of increasing τ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 τ2. This calculation means that the maps in Figs 14(b) and (d) provide the prediction standard deviation for an unobserved κ0 value. Figure 14. View largeDownload slide (a) Median κ0 map for Model 1 with θ = 0.5 and (b) its standard deviation in log10(s). (c) The median model for Model 2 with θ = 0.5 and (d) its standard deviation. Figure 14. View largeDownload slide (a) Median κ0 map for Model 1 with θ = 0.5 and (b) its standard deviation in log10(s). (c) The median model for Model 2 with θ = 0.5 and (d) its standard deviation. 6.4 Consideration of κ0 ‘measurement uncertainty’ Fig. 14 clearly demonstrates that the nugget is an important parameter when considering forward application of the κ0 map. While the seismograph network in New Zealand is insufficiently dense to derive the nugget parameter from the kriging analysis alone, estimates of κ0 measurement uncertainty are available from the linear regression of the κ(r) data with distance (Figs 5 and 6). This is the purpose of model 2 in Table 2. The fixed τ2 value of 0.03 for model 5 is only an indicative value of the κ0 uncertainty across all stations. A limitation of this assumption is that some of the κ0 data are better constrained than others, as evidenced by the confidence intervals in Figs 5 and 6. Ideally, a fixed vector of τ2 with each station’s κ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 κ0 map. The station-specific κ0 are derived by regression of κ data with distance on a linear scale, while a lognormal distribution was assumed to derive the κ map. The combination of distributions prevents an elegant determination of the spatially varying nugget, hence the fixed nugget of 0.03 (± 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 κ0 to constrain the nugget, which are mostly unavailable in New Zealand. For model 1, the predicted standard deviation on κ0 is around 0.1 log10 units, while for model 2 it is around 0.2 log10 units. Note that this 0.2 is only slightly larger than the value that was initially assumed for the uncertainty in the station-specific κ0 values, which suggests that the uncertainty in calculating κ0 at a given station location is much larger than the uncertainty introduced by kriging. 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 QP 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 physical 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. https://doi.org/10.1007/s000240050204 Google Scholar CrossRef Search ADS   Abrahamson N., Silva W., 1997. Empirical response spectral attenuation relations for shallow crustal earthquakes, Seismol. Res. Lett. , 68( 1), 94– 127. https://doi.org/10.1785/gssrl.68.1.94 Google Scholar CrossRef Search ADS   Abrahamson N., Silva W., Kamai R., 2014. Summary of the ASK14 ground motion relation for active crustal regions, Earthq. Spectra , 30( 3), 1025– 1055. https://doi.org/10.1193/070913EQS198M Google Scholar CrossRef Search ADS   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. https://doi.org/10.1029/98JB01757 Aki K., 1967. Scaling law of seismic spectrum, J. geophys. Res. , 72( 4), 1217– 1231. https://doi.org/10.1029/JZ072i004p01217 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1016/0031-9201(80)90019-9 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1785/gssrl.83.1.166 Google Scholar CrossRef Search ADS   Atkinson G., Boore D., 2006. Earthquake ground-motion prediction equations for eastern North America, Bull. seism. Soc. Am. , 96( 6), 2181– 2205. https://doi.org/10.1785/0120050245 Google Scholar CrossRef Search ADS   Bannister S., Gledhill K., 2012. Evolution of the 2010–2012 Canterbury earthquake sequence, N.Z. J. Geol. Geophys. , 55( 3), 295– 304. https://doi.org/10.1080/00288306.2012.680475 Google Scholar CrossRef Search ADS   Benjamin D. et al.  , 2018. Redefine statistical significance, Nat. Hum. Behav. , 2( 1), 6– 10. Google Scholar CrossRef Search ADS   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. https://doi.org/10.1007/PL00012553 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1193/070113EQS184M Google Scholar CrossRef Search ADS   Brune J., 1970. Tectonic stress and the spectra of seismic shear waves from earthquakes, J. geophys. Res. , 75( 26), 4997– 5009. https://doi.org/10.1029/JB075i026p04997 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1785/0120020002 Google Scholar CrossRef Search ADS   Campbell K., 2009. Estimates of shear-wave Q and κ0 for unconsolidated and semiconsolidated sediments in eastern North America, Bull. seism. Soc. Am. , 99( 4), 2365– 2392. https://doi.org/10.1785/0120080116 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1193/062913EQS175M Google Scholar CrossRef Search ADS   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. https://doi.org/10.1016/S0040-1951(98)00277-7 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1016/j.soildyn.2006.02.010 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1193/072813EQS219M Google Scholar CrossRef Search ADS   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. https://doi.org/10.1007/s10950-005-9006-7 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1029/GL008i011p01126 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1142/S1363246906002566 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1785/0120080075 Google Scholar CrossRef Search ADS   Douglas J., Gehl P., Bonilla L. F., Gélis C., 2010. A κ model for mainland France, Pure appl. Geophys. , 167( 11), 1303– 1315. https://doi.org/10.1007/s00024-010-0146-5 Google Scholar CrossRef Search ADS   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. https://doi.org/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. https://doi.org/10.1785/0120030036 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1111/j.1365-246X.2005.02530.x Google Scholar CrossRef Search ADS   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. https://doi.org/10.1029/2007JB005359 Eberhart-Phillips D., Reyners M., Bannister S., 2015. A 3D QP attenuation model for all of New Zealand, Seismol. Res. Lett. , 86( 6), 1655– 1663. https://doi.org/10.1785/0220150124 Google Scholar CrossRef Search ADS   Edwards B., Ktenidou O.-J., Cotton F., Abrahamson N., Van Houtte C., Fäh D., 2015. Epistemic uncertainty and limitations of the κ0 model for near-surface attenuation at hard rock sites, Geophys. J. Int. , 202( 3), 1627– 1645. https://doi.org/10.1093/gji/ggv222 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1785/0120100057 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1093/gji/ggu030 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1002/2014JB011642 Google Scholar CrossRef Search ADS   Ktenidou O.-J., Abrahamson N., 2016. Empirical estimation of high-frequency ground motion on hard rock, Seismol. Res. Lett. , 87( 6), 1465– 1478. https://doi.org/10.1785/0220160075 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1785/0120120093 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1785/0220130027 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1093/gji/ggv315 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1785/0120150301 Google Scholar CrossRef Search ADS   Laurendeau A., Cotton F., Ktenidou O.-J., Bonilla L.-F., Hollender F., 2013. Rock and stiff soil site amplification: dependency on vS30 and kappa (κ0), Bull. seism. Soc. Am. , 103( 6), 3131– 3148. https://doi.org/10.1785/0120130020 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1029/94GL00977 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1785/0220130009 Google Scholar CrossRef Search ADS   Mooney H., 1970. Upper mantle inhomogeneity beneath New Zealand: seismic evidence, J. geophys. Res. , 75( 2), 285– 309. https://doi.org/10.1029/JB075i002p00285 Google Scholar CrossRef Search ADS   Morozov I., 2008. Geometrical attenuation, frequency dependence of Q, and the absorption band problem, Geophys. J. Int. , 175( 1), 239– 252. https://doi.org/10.1111/j.1365-246X.2008.03888.x Google Scholar CrossRef Search ADS   Motazedian D., Atkinson G., 2005. Stochastic finite-fault modeling based on a dynamic corner frequency, Bull. seism. Soc. Am. , 95( 3), 995– 1010. https://doi.org/10.1785/0120030207 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1007/s00024-013-0751-1 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1785/0120100135 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1785/0120080029 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1785/0120160374 Google Scholar CrossRef Search ADS   Petersen T., Gledhill K., Chadwick M., Gale N., Ristau J., 2011. The New Zealand national seismograph network, Seismol. Res. Lett. , 82( 1), 9– 20. https://doi.org/10.1785/gssrl.82.1.9 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1785/0120100345 Google Scholar CrossRef Search ADS   Van Houtte C., Ktenidou O., Larkin T., Holden C., 2014. Hard-site κ0 (kappa) calculations for Christchurch, and comparison with local ground-motion prediction models, Bull. seism. Soc. Am. , 104( 4), 1899– 1913. https://doi.org/10.1785/0120130271 Google Scholar CrossRef Search ADS   Wessel P., Smith W., Scharroo R., Luis J., Wobbe F., 2013. Generic mapping tools: improved version released, Eos Trans. AGU , 94( 45), 409– 410. https://doi.org/10.1785/0120130271 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1016/0377-0273(95)00006-G Google Scholar CrossRef Search ADS   Yogomida K., Benites R., 1995. Relation between direct wave Q and coda Q: a numerical approach, Geophys. J. Int. , 123( 2), 471– 483. https://doi.org/10.1111/j.1365-246X.1995.tb06866.x Google Scholar CrossRef Search ADS   Yoshimoto K., Okada M., 2009. Frequency-dependent attenuation of S-waves in the Kanto region, Japan, Earth, Planets Space , 61( 9), 1067– 1075. https://doi.org/10.1186/BF03352957 Google Scholar CrossRef Search ADS   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. https://doi.org/10.1111/j.1365-246X.1993.tb01476.x Google Scholar CrossRef Search ADS   SUPPORTING INFORMATION Supplementary data are available at GJIRAS online. Figure S1. Effect of changing η values when fitting the high-frequency spectral shape. Note that large variations in Q0 and η are possible with little effect on the quality of the fit. Figure S2. Illustration of the degree of variability in the mean κ values (red points), as shown in the ±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 κ0 values, one could either (i) separate data into two regions, calculate two slopes and κ0 values, then take the average κ0, or (ii) combine the data together with a single slope and κ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’ κ0 map essentially smooths over the station-specific κ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. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Apr 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off