Hybrid broadband ground-motion simulation validation of small-magnitude subduction earthquakes in New Zealand

Abstract

Physics-based simulation of subduction earthquake ground motions remains less comprehensively validated than for shallow crustal earthquakes, despite subduction events contributing significantly to global seismic hazard. In this study, subduction-specific simulation models were developed and validated for small-magnitude ( $M_w$ 3.5–5) interface and slab earthquakes using hybrid broadband ground-motion simulation. Simulation models were constrained by (1) global empirical ground-motion models, (2) a global database of finite-fault rupture models, and (3) global ground-motion simulation studies. The simulation models include subduction source-specific representations for stress parameter and rupture velocity, with a significant depth dependence of the stress parameter for slab earthquakes. Volcanic backarc effects on anelastic path attenuation are also included through path attenuation-based scaling of the rock quality factors in the simulations. The simulations leverage models for site and basin effects which have been validated using crustal earthquakes, for which there are many recorded ground motions, and the subduction-specific model modification mainly affect the high-frequency component of the ground-motion simulations above 1Hz. Simulation predictions were validated against a compiled dataset of observed subduction earthquake ground-motion records in New Zealand. The subduction-specific modifications significantly improve predictive performance compared to the use of simulation parameter values for active shallow crustal earthquakes and provide comparable accuracy to prior simulations for crustal earthquakes in New Zealand. Despite these promising results, further studies are needed to address prediction residuals at low frequencies (f ≤ 1 Hz), as well as extending validation to larger magnitude earthquakes.

Introduction

Subduction earthquakes contribute significantly to global seismic hazard; however, comprehensive validation for physics-based ground-motion simulation of subduction earthquakes has received less attention compared to active shallow crustal earthquakes (e.g. Goulet et al., 2015; Graves and Pitarka, 2010, 2015, 2016; Olsen and Takedatsu, 2015). A likely contributing factor is that ground-motion prediction for subduction earthquakes, using either empirical ground-motion models (GMMs) or physics-based simulations, is more challenging than for crustal earthquakes. Subduction earthquakes occur at depth and often offshore or along coastlines; therefore, ground-motion records from these events tend to have large rupture distances and poor azimuthal representation. Furthermore, there is greater evidence of region-to-region variability of subduction earthquake-induced ground motions and therefore a simple combination of global data is not appropriate (e.g. Bozorgnia et al., 2022).

It is perhaps due to these aforementioned challenges that subduction earthquake simulation studies have primarily focused on large-magnitude interface earthquakes, such as the 2003 Tokachi-Oki $M_w$ 8.3 (Iwaki et al., 2016a, 2016b; Wirth et al., 2017), 2010 Maule $M_w$ 8.8 (Frankel, 2017), and 2011 Tohoku $M_w$ 9.0 (Ghofrani et al., 2013; Goda et al., 2017) events, and studies of potential rupture scenarios for $M_w$ 9.0 earthquakes on the Cascadia Subduction Zone (Frankel et al., 2018; Wirth et al., 2018) and a $M_w$ 8.6 earthquake on the Hikurangi Subduction Zone (Bayless et al., 2019), among others.

Although forensic examination of significant historical earthquakes serve as useful case studies, and reveal insights for ground-motion simulations of potential future events, the performance of ground-motion simulations must be comprehensively validated using observed ground-motion records prior to the implementation of simulated ground motions in seismic hazard analysis (e.g. Bradley et al., 2017; Graves et al., 2011) or engineering applications (e.g. Bijelić et al., 2018; Galasso et al., 2013). Consideration of small- and moderate-magnitude events offers the possibility to consider many more earthquakes and observed ground motions and thus achieve more predictive confidence. In particular, small-magnitude events ( $M_w\le 5$ ) offer an ideal domain for validation because (1) there are many observed ground-motion records, (2) their ruptures are generally accurately represented with simple point-source models, and (3) they produce relatively low-amplitude ground motions which do not induce a nonlinear site response. The self-similar nature of earthquake ruptures, which is well-established for small-to-moderate magnitudes (e.g. Somerville et al., 1999), allows for models that are validated in the small-magnitude range to be readily extended to moderate magnitudes ( $M_w$ 5–7) where validation rigor is impinged by a lack of data.

New Zealand (NZ) provides a convenient natural laboratory to validate the performance of physics-based ground-motion simulation. It is an active seismic region that experiences shallow crustal earthquakes as well as subduction interface and slab earthquakes along two subduction zones. This seismicity, coupled with a dense network of strong ground-motion instruments, has yielded a wealth of strong ground-motion data (e.g. Hutchinson et al., 2024). As an example, Lee et al. (2020) conducted a simulation validation study for crustal earthquakes using 1896 observed ground-motion records from 148 small-magnitude ( $M_w$ 3.5–5) events in the Canterbury, NZ, region and demonstrated predictive accuracy similar to empirical GMMs. Subsequently, Lee et al. (2022) extended this validation to a comprehensive set of 5218 ground motions recorded at 212 sites from 479 crustal earthquakes across NZ with modified simulation parameter values inferred from previous NZ validation.

Naturally, such a validation and modification process should also address interface and slab earthquakes that warrant different simulation models and parameters than for shallow crustal earthquakes. This paper presents a validation study of hybrid broadband ground-motion simulation for small-magnitude ( $M_w$ 3.5–5) earthquake ground-motion records from NZ with model-specific adjustments for subduction interface and slab earthquakes. Subduction-specific parameter models are developed based on examining the main features of empirical GMMs, analyzing a database of finite-fault source rupture models, and considering studies of subduction earthquakes and subduction earthquake ground-motion simulations by others. Ground-motion simulations are validated using observed ground-motion records from interface and slab earthquakes. Mixed-effects regression and residual partitioning are used to identify systematic rupture characteristics for subduction earthquakes and to validate the performance of the subduction models. Partitioned simulation prediction residuals are examined to identify observed systematic effects of the model-specific adjustments. Residual dependence on causal parameters and geospatial residual trends are then interrogated to support the subduction-specific models. Additional details of the analysis and simulation predictions are provided in the Electronic Supplement along with other supporting material.

Methods

Observed ground-motion data for validation

Observed ground motion data from the 2023 NZ Ground-Motion Database (Hutchinson et al., 2024) was used whereby events were classified using a modified hypocentre-based tectonic classification logic of Bozorgnia et al. (2020) with the Hikurangi (Williams et al., 2013) and Puysegur (Hayes et al., 2018) subduction interface geometries. The classification of each centroid moment tensor (CMT) solution (Ristau, 2008, 2013) was then manually reviewed considering the focal mechanism location and orientation relative to the adjacent subduction interface geometry. $M_w$ 3.5–5 earthquakes were considered in this study; smaller magnitudes were not considered because the associated CMTs are fewer in number, less accurate, and tend to have fewer associated ground-motion records (Ristau, 2018). Simulation validation for $M_w>5$ events was left for a subsequent study. This gave rise to $160$ and $323$ candidate interface and slab events, respectively.

Only observed ground-motion records with two high-quality horizontal components and maximum usable vibration periods of T ≥ 3 s were considered (Dupuis et al., 2023)—hereinafter referred to as a “high-quality” ground-motion record. A minimum of three high-quality ground-motion records were required for each event and site to be included. This is consistent with Lee et al. (2022), and based on balancing the competing objectives of facilitating statistically significant residuals of ground-motion IMs while also accommodating a large geographical area with relatively sparse distribution of available events.

Table 1 shows the number of events, sites, and records for interface and slab earthquakes along the Hikurangi and Puysegur subduction zones, and compares with the crustal ground-motion validation study of Lee et al. (2022). As shown in Figure 1, for interface and slab earthquakes, respectively, there is good geographical coverage of events and sites along the Hikurangi Subduction Zone. The majority of the sites included are located in the forearc region because of increased instrumentation there and because of the higher attenuation for backarc travel paths. There are fewer high-quality ground-motion records associated with the Puysegur Subduction Zone, as shown in Table 1, due its relatively few sites and few CMTs for that region (Ristau, 2008, 2013). Details of the event, site, and record distributions are presented in the Electronic Supplement.

Table 1.

Summary of high-quality observed ground-motion data used for validation

	Interface			Slab
	Hikurangi	Puysegur	Total	Hikurangi	Puysegur	Total	Crustal a
Events	$59$	$3$	$62$	$132$	$15$	$147$	$496$
Sites	$63$	$3$	$66$	$116$	$8$	$124$	$382$
Records	$358$	$9$	$367$	$909$	$66$	$975$	$5218$

a

Reference crustal ground-motion simulations by Lee et al. (2022).

Figure 1.

Small-magnitude subduction earthquakes, ray paths, and sites for high-quality ground-motion records used for validation of ground-motion simulations. Depth contours of the subduction zone interfaces are shown in $10$ km increments with the surface ( $0$ km depth) contour shown in bold.(a) Interface data. (b) Slab data.

Ground-motion simulation approach and inputs

This study uses the hybrid broadband ground-motion simulation approach of Graves and Pitarka (2010, 2015, 2016). Hybrid broadband ground-motion simulation is a physics-based method which combines a comprehensive three-dimensional (3D) wave propagation solution at low frequencies (LF) with a simplified physics solution at high frequencies (HF). The LF and HF components are combined using a set of matched fourth-order Butterworth filters to produce a broadband ground motion. In this study, a finite difference grid spacing of $100$ m is used with a minimum shear wave velocity, VS = 500 m/s, resulting in an LF-HF transition frequency of $1$ Hz (Graves and Pitarka, 2010).

This study considers kinematic point-source rupture models for each CMT, which adequately represent small-magnitude sources, with prescribed slip, rake, and initiation times for both the LF and HF components. For the LF component, the theoretical wave propagation solution from each source rupture is computed using a staggered-grid finite difference method (Graves, 1996) and a 3D velocity model with domain extents optimized for land coverage encompassing sites with high-quality records. Version 2.07 of the NZ Velocity Model (NZVM; Thomson et al., 2020), which builds on the background travel-time tomography-based seismic velocity model of Eberhart-Phillips et al. (2010), was used. For the HF component, S-wave travel times, frequency-dependent attenuation, and quarter-wavelength impedance effects are computed using a one-dimensional (1D) velocity model. A summary of the simulation models pertaining to this specific study is included in the Electronic Supplement and comprehensive details of the hybrid broadband ground-motion simulation approach can be found in the work by Graves and Pitarka (2010, 2015, 2016).

Following Lee et al. (2022), we adopt minor modifications to the hybrid broadband simulation approach to better suit NZ conditions. These changes were motivated by observations from Lee et al. (2020), which found that simulations systematically under-predicted ground motion significant duration and over-predicted short-period response spectra. Modifications include a path duration model of Boore and Thompson (2014), developed for crustal earthquakes, which produces larger predictions of significant duration than the original crustal reference model used in Lee et al. (2020) for all rupture distances. The $V_{S30}$ -based empirical site amplification factor model of Campbell and Bozorgnia (2014) was removed for the LF component because long-period site effects are at least partially included through explicit treatment in the NZVM used for the LF simulations, which includes sedimentary basins and a shallow geotechnical layer to represent weathered rock and sedimentary soils (where sedimentary basins are not explicitly included with embedded models). The appropriate treatment of long-period site effects is an important consideration and is the subject of dedicated studies (e.g. Kuncar et al., 2025) and other forthcoming publications. These modifications are also appropriate for subduction simulations, which use the same NZVM, and for which other measures of duration indicate that although there are differences from crustal events (e.g. Bahrampouri et al., 2021), these are small relative to the aleatory uncertainty (e.g. Baltay et al., 2017).

Partially crossed linear mixed-effects regression

A partially crossed linear mixed-effects regression approach is applied to partition prediction residuals into various components of variability (Bates et al., 2014; Stafford, 2014). The language and notation used is that of Atik et al. (2010) for ground-motion prediction validation which follows the general form of a GMM for an event, $e$ , and site, $s$ , pairing, explained subsequently. The general and expanded forms of the equation are:

$$\mathrm{\Delta }=\ln I{M}_{\mathrm{es}}-f_{\mathrm{es}}=a+\delta B_e+\delta S2S_s+\delta {W^0}_{\mathrm{es}}$$

(1)

where $\mathrm{\Delta }$ is the total prediction residual; $\ln I{M}_{\mathrm{es}}$ is the natural logarithm of the observed intensity measure (IM) for event $e$ and site $s$ ; and $f_{\mathrm{es}}$ is the median (for the case without parameter uncertainty there is a single IM prediction) of the predicted logarithmic IM either from a simulation or empirical GMM. $a$ is the model bias; $\delta B_e$ is the between-event residual for event $e$ , with zero mean and variance ${\tau }^2$ . $\delta S2S_s$ is the systematic site-to-site residual for site $s$ , with zero mean and variance ${{\varphi }_{S2S}}^2$ ; $\delta {W^0}_{\mathrm{es}}$ is the “remaining” within-event residual with zero mean and variance ${{\varphi }_{\mathrm{SS}}}^2$ . The total residual, $\mathrm{\Delta }$ , has a variance:

$${\sigma }^2={\tau }^2+{{\varphi }_{S2S}}^2+{{\varphi }_{\mathrm{SS}}}^2$$

(2)

Although, Equations 1 and 2 are presented with the notation of Atik et al. (2010) for brevity, the mixed-effects method is also applied to analyze a database of finite-fault rupture models in the subsequent Analysis of Global SRCMOD Database of Rupture Models section.

Constraining subduction-specific simulation models and parameters

Modifications to the simulation models for application to subduction earthquakes were constrained by three primary methods: (1) analysis of global empirical GMMs for subduction earthquakes, (2) analysis of a global database of finite-fault rupture models (SRCMOD) created from source inversions of historical earthquakes, and (3) analysis of global ground-motion simulation studies which have considered subduction earthquakes. While each of these methods offers only partial support for individual model refinements, collectively they inform the adopted subduction-specific simulation models.

Key features of empirical GMMs

Global empirical GMMs for crustal and subduction earthquakes were reviewed to elucidate systematic characteristics of subduction earthquake ground motions. The empirical GMMs considered herein are from the NGA-West2 (Bozorgnia et al., 2014) and NGA-Subduction (Bozorgnia et al., 2022) projects, which address crustal and subduction earthquakes, respectively, and include additional GMMs for crustal earthquakes by Abrahamson et al. (2014) and Bradley (2013). Due to the simplified source and path parametrization within empirical GMMs, it is not possible to make a direct link from empirical GMMs to simulation parameters for subduction earthquakes. However, two phenomenological features of subduction earthquakes which are used in ground-motion simulations are closely related to empirical GMM parametrizations: (1) stress parameter is strongly related to source depth and (2) backarc effects on attenuation are typically represented through volcanic arc-delineated regionalization.

Figure 2 compares the median IM predictions from these GMM suites for crustal, interface, and slab earthquakes for pseudo-spectral accelerations (pSAs) at selected periods. Representative reference depths for each earthquake type and a site $V_{S30}$ of $750$ m/s were used; the observed trends were robust to these selections. The comparison indicates that slab ground motions have larger short-period (e.g. T = 0.1 s) pSAs compared to crustal ground motions and smaller long-period (e.g. T = 2.0 s) pSAs; interface ground motions demonstrate comparable pSAs at short periods, except for very large magnitudes, and smaller pSAs at long periods.

Figure 2.

Comparison of (a, b) empirical GMM pSA predictions for crustal, interface, and slab earthquakes and (c, d) the ratio of interface and slab to crustal earthquake predictions by (a, c) $R_{\mathrm{rup}}$ and (b, d) $M_w$ . The median predictions are shown for $R_{\mathrm{rup}}$ between 100– $500$ km and $M_w$ 5–8 for suites of GMMs for crustal (Abrahamson et al. (2014); Bradley (2013); Boore et al. (2014); Campbell and Bozorgnia (2014); and Chiou and Youngs (2014), interface, and slab (Abrahamson et al. (2016); Abrahamson and Gülerce (2020); Kuehn et al. (2020); and Parker et al. (2020)) earthquakes. Reference depths of $15$ km for crustal, $35$ km for interface, and $75$ km for slab are considered with $V_{S30}$ of $750$ m/s. (a) Predicted pSAs by $R_{\mathrm{rup}}$ . (b) Predicted pSAs by $M_w$ . (c) Predictions relative to crustal by $R_{\mathrm{rup}}$ . (d) Predictions relative to crustal by $M_w$ .

Significant differences in the observed interface and slab earthquake ground motions have led to the development of separate source models with unique depth scaling as summarized for selected GMMs in Table 2. For these GMMs, depth scaling was qualitatively characterized through review of their functional forms and examination of prediction trends with source depth. Such GMMs generally exhibit greater depth scaling for slab than for interface earthquakes (e.g. Abrahamson and Gülerce, 2020) and stronger scaling of short-period pSAs than of long-period pSAs (e.g. Kuehn et al., 2020). This depth-scaling behavior is examined in further detail in Figure 3 based on depth to top of rupture, $Z_{\mathrm{TOR}}$ ; depth scaling is most significant for high-frequency IMs and for slab earthquake ground motions, with practically no effect for crustal earthquake ground motions. These differences in depth scaling may result from the influence of disparate geophysical conditions within and adjacent to subducting oceanic slabs compared with continental crust.

Table 2.

Qualitative summary of depth-dependent scaling of selected GMMs for interface and slab earthquakes: AG20 (Abrahamson and Gülerce, 2020), C20 (Chao et al., 2020), HA21 (Hassani and Atkinson, 2021), K20 (Kuehn et al., 2020), P20 (Parker et al., 2020), and S20 (Si et al., 2020)

	Depth-dependent scaling
Ground motion model	Interface	Slab
AG20	None	<50 km: Strong
		50–200 km: Moderate
		>200 km: None
C20	Slight linear	Strong linear
HA21	Slight linear	Strong linear
K20	<30 km: Slight linear>30 km: None	<80 km: Strong linear>80 km: None
P20	Slight linear	<20 km: Strong20–67 km: Moderate>67 km: Slight
S20	None	Strong linear

For all GMMs, depth dependence results in increased IMs for increased depths when all other parameters are held constant.

Figure 3.

Comparison of (a, b) empirical GMM predictions for crustal, interface, and slab earthquakes at various depths to top of rupture, $Z_{\mathrm{TOR}}$ and (c, d) the ratio of these predictions to reference predictions at $Z_{\mathrm{TOR}}$ of $30$ km for (a, c) pSA ( $0.01$ s) and (b, d) pSA ( $0.5$ s). Median predictions are shown for $R_{\mathrm{rup}}$ of $200$ km, $M_w$ 5, and $V_{S30}$ of $750$ m/s for selected GMMs: A14 (Abrahamson et al., 2014), AG20 (Abrahamson and Gülerce, 2020), HA21 (Hassani and Atkinson, 2021), K20 (Kuehn et al., 2020), and P20 (Parker et al., 2020); the mean for each GMM suite is also included. (a) Predicted pSA ( $0.01$ s) by $Z_{\mathrm{TOR}}$ . (b) Predicted pSA ( $0.5$ s) by $Z_{\mathrm{TOR}}$ . (c) Predicted pSA ( $0.01$ s) at various $Z_{\mathrm{TOR}}$ relative to predictions at $Z_{\mathrm{TOR}}$ of $30$ km. (d) Predicted pSA ( $0.5$ s) at various $Z_{\mathrm{TOR}}$ relative to predictions at $Z_{\mathrm{TOR}}$ of $30$ km.

In addition to the NGA-West2 and NGA-Sub models, the GMM of Hassani and Atkinson (2021) includes equivalent point-source models for crustal, interface, and slab earthquakes in Japan and is particularly useful because a source-specific term for stress parameter is explicitly included, although it is not directly equivalent to the stress parameter values used in the simulations performed here (Atkinson and Beresnev, 1997). The models include depth-dependent scaling of stress parameter with greater stress parameter and depth scaling for slab earthquakes than for interface earthquakes, and with a strong effect on the resulting high-frequency pSAs (Figure 3).

For both interface and slab earthquakes, there is greater anelastic attenuation (relative to crustal) for travel paths which pass through the zone of partial melting in the volcanic backarc, especially for deeper events (e.g. Cousins et al., 1999; Hassani and Atkinson, 2021) and significantly greater backarc attenuation (relative to forearc) has been identified in several regions including Japan, Romania, and NZ (e.g. Beauval et al., 2017; Vacareanu et al., 2015), among others. Conversely, the relatively cold and competent subducting slab provides an up-dip and along-strike travel path toward the forearc region with relatively little anelastic attenuation (Skarlatoudis et al., 2013).

Differences in forearc and backarc attenuation are handled differently by different GMMs; some GMMs address the issue by only considering ground motion in the forearc region. For example, a GMM for subduction earthquakes in Japan was developed by Si et al. (2020) with separate treatment of interface and slab earthquakes, which only considers forearc sites. Parker et al. (2020) excluded data recorded at sites in the backarc but hypothesizes that backarc anelastic attenuation was greater than in the forearc and also potentially more variable between regions. Similarly, Abrahamson and Gülerce (2020) developed a set of global and region-specific GMMs using the regionalization of Bozorgnia et al. (2022) which does not consider sites in the backarc region. Some GMMs have applied simplified approaches; Chao et al. (2020) identified very different backarc anelastic attenuation between crustal, interface, and slab earthquakes; however, due to limited data, they derived the same anelastic attenuation coefficients for all subduction earthquakes based on a similar approach by Abrahamson et al. (2016). Phung et al. (2020) refit and validated the GMM of Abrahamson et al. (2016) for Taiwan using data resources from the NGA-Subduction project for Japan and Taiwan. Phung et al. (2020) did not address the difference in attenuation between forearc and backarc regions due to metadata limitations and subsequently identified this as the major limitation of their GMM.

Finally, some GMMs, as shown in Figure 4, have addressed backarc attenuation directly through regionalized or arc crossing-based attenuation terms. Abrahamson et al. (2016) used an arc crossing-based anelastic attenuation term that was applied to lower overall ground-motion amplitudes (simple offset) in the backarc. In the equivalent point-source model of Hassani and Atkinson (2021), period-dependent crustal-, interface-, and slab-specific models for anelastic attenuation are used based on the site location, with greater attenuation for travel paths toward the backarc. Finally, in a GMM for slab earthquakes in Romania (not shown), Sokolov et al. (2008) developed regional attenuation models based on source-to-site azimuth which were found to agree well with observed ground-motion IMs. As shown in Figure 4, both Abrahamson et al. (2016) and Hassani and Atkinson (2021) provide similar treatment for backarc sites, with a greater effect for HF and similar treatment regardless of tectonic classification and source depth.

Figure 4.

Comparison of empirical GMM predictions for interface (a, c) and slab (b, d) earthquakes for sites in the forearc and backarc at various depths to top of rupture, $Z_{\mathrm{TOR}}$ for (a, b) pSA ( $0.02$ s) and (c, d) pSA ( $2.0$ s). Predictions are shown for $R_{\mathrm{rup}}$ of $200$ km, $M_w$ 5, and $V_{S30}$ of $750$ mãs for selected GMMs which have bespoke treatment for backarc and forearc sites: A16 (Abrahamson et al., 2016) and HA21 (Hassani and Atkinson, 2021). The arithmetic mean predictions of the two GMMs (A16 and HA21) are also shown. (a) Predicted pSA ( $0.02$ s) for interface. (b) Predicted pSA ( $0.02$ s) for slab. (c) Predicted interface pSA ( $2.0$ s) for interface. (d) Predicted pSA ( $2.0$ s) for slab.

Analysis of global SRCMOD database of rupture models

To identify systematic features of subduction interface and slab ruptures, the SRCMOD rupture model database (Mai and Thingbaijam, 2014) was analyzed. SRCMOD contained $423$ rupture models which represent the spatio-temporal ruptures of global earthquakes with magnitudes of $M_w$ 4.6–9.2. The rupture models are provided without complete tectonic classifications; therefore, to examine variations between crustal, interface, and slab earthquakes, each model was manually classified based on a review of its centroid location and focal mechanism relative to global subduction interfaces (Hayes et al., 2018). Some earthquakes have multiple associated rupture models from a given author; in these cases, the most recent model from that author was used and the earlier models were discarded—in total, 51 duplicates were removed. The labeled database contained $206$ crustal, $197$ interface, and $20$ slab rupture models.

The models in SRCMOD are based on source inversions compiled from different authors who used various inversion methods. The finite-fault rupture models in the database are not equivalent to the source rupture models used as inputs in physics-based ground-motion simulation which have prescribed kinematic ruptures with spatio-temporal description of slip, rise time, and rake at a significantly greater spatial resolution. Therefore, it is not appropriate to directly extend data synthesis from the database to ground-motion simulation. Similarly, direct comparisons between the average values for each tectonic classification are also imperfect because the average magnitudes and depths differ between the crustal, interface, and slab models. Despite these challenges, it is still possible to compare the aggregated rupture properties of each tectonic classification, as is done in the following sections, using a mixed-effects regression analysis which examines relative differences and trends with magnitude and depth. Although only $11$ of the $423$ models are for small-magnitude ( $M_w\le 5$ ) earthquakes, relative differences between the crustal, interface, and slab earthquakes are anticipated to generalize well to small magnitudes due to the self-similar nature of earthquake ruptures in the small-to-moderate magnitude range (e.g. Somerville et al., 1999).

Various parameters are used to describe the kinematic rupture models in the database: for example, magnitude, area, rise time, rupture velocity, dip, slip, rake, and geospatial location (latitude, longitude, and depth). These are important for finite-fault simulations (e.g. dip, slip, rake, and area); however, this study focuses on small-magnitude earthquakes that are modeled as point sources and so only rise time and rupture velocity were investigated. Rise time is inversely related to stress parameter which controls high-frequency energy release (e.g. Boore, 1983). Rupture velocity is a simulation parameter used for both the high- and low-frequency components in the hybrid broadband ground-motion simulation approach (Graves and Pitarka, 2010).

Crustal-based parameter values for ground-motion simulations (Graves and Pitarka, 2015) of rise time (in seconds), ${\tau }_A$ , and rupture velocity ratio (relative to shear wave velocity, $V_S$ , at the source), $V_R/V_S$ , were used as reference predictions in the mixed-effects regression analysis:

$${\tau }_A={\alpha }_T\times 1.45\times {10}^{-9}\times {M_o}^{1/3}$$

(3)

$$V_R/V_S=0.8$$

(4)

where ${\alpha }_T$ is a mechanism-dependent scaling factor computed based on dip and rake, $M_o$ is the seismic moment (in Nm), and $V_S$ is the shear wave velocity at the source. The regression was performed with and without shallow and deep reduced rupture velocity zones as prescribed in the Graves and Pitarka (2016) model, and the observed trends in the residuals were affected very slightly; thus, only the case without these zones was considered thereafter. A reference rupture velocity prediction of VR = 2.5 km/s was used for the analysis of rupture velocity for all tectonic classifications.

The partially crossed linear mixed-effects regression method (i.e. Equations 1 and 2) was applied to the database with the SRCMOD models treated as observed data for an event, $e$ , and source author, $s$ , pairing. The general and expanded forms of the equation in the context of this analysis of SRCMOD rupture parameters are:

$$\mathrm{\Delta }=\ln P_{\mathrm{es}}-\ln f_{\mathrm{es}}=a+\delta B_e+\delta S2S_s+\delta {W^0}_{\mathrm{es}}$$

(5)

where $\mathrm{\Delta }$ is the total prediction residual; $\ln P_{\mathrm{es}}$ is the natural logarithm of the parameter value from the SRCMOD rupture model; and $\ln f_{\mathrm{es}}$ is the natural logarithm of the reference crustal-based parameter value prediction for event, $e$ , and source author, $s$ . Once decomposed into fixed and random effects, $a$ is the bias relative to the reference crustal-based model; $\delta B_e$ is the between-event residual for event $e$ ; $\delta S2S_s$ is the systematic source author-to-source author residual; $\delta {W^0}_{\mathrm{es}}$ is the “remaining” within-event residual which represents other factors.

Because some inversions may poorly represent reality, models with greater than 2 standard deviations (SDs) from the median value in the database were removed to improve the stability of the inferences. The analysis was done collectively for the entire database with one computed value of bias, $a$ , for all tectonic classifications due to the limited amount of data available. The variations of the between-event residuals, $\delta B_e$ , with magnitude and hypocentre depth for these three tectonic classifications is shown in Figure 5 and 6.

Figure 5.

Between-event residuals, $\delta B_e$ , and moving averages for rise time, ${\tau }_A$ , from the SRCMOD database by (a) hypocentre depth and (b) magnitude. Moving averages from a separate regression with the subduction-specific models (Subduction-Specific Simulation Parameters section) used as predictions for subduction earthquakes instead of the reference crustal models are shown with thin gray lines. (a) ${\tau }_A$ residuals by depth. (b) ${\tau }_A$ residuals by magnitude.

Figure 6.

Between-event residuals, $\delta B_e$ , and moving averages for (a, b) rupture velocity, $V_R$ , from the SRCMOD database (in km/s) and (c, d) rupture velocity as a ratio of shear wave velocity at the source, $V_R/V_S$ by (a, c) hypocentre depth and (b, d) magnitude. Moving averages from a separate regression with the subduction-specific models (Subduction-Specific Simulation Parameters section) used as predictions for subduction earthquakes instead of the reference crustal models are shown with thin gray lines. (a) $V_R$ residuals by depth. (b) $V_R$ residuals by magnitude. (c) $V_R/V_S$ residuals by depth. (d) $V_R/V_S$ residuals by magnitude.

The inversion-based source models are inherently limited in their ability to capture short-wavelength features due to spatial smoothing constraints; thus, the analysis reflects model-based representations rather than direct observations. By comparing between-event residuals of rise time and rupture velocity and their dependence on depth and magnitude, it is possible to infer systematic differences between crustal, interface, and slab ruptures within the SRCMOD database. It is useful to consider magnitude scaling in tandem with depth scaling due to strong multicollinearity of hypocentre depth and magnitude observed for slab models and because of disparities in magnitudes and depths for the crustal, interface, and slab models in the database.

Figure 5 illustrates there is a systematic relationship of the rise time between-event residuals with magnitude for all three tectonic classifications and shows a trend toward over-prediction of the crustal-based simulation rise times relative to the SRCMOD database at larger magnitudes. This trend maps to a similar trend for hypocentre depth due to the multicollinearity of depth and magnitude in the database. Despite the general misfit of the simulation-based predictions for rise times of the finite-fault rupture models, the results indicate a pattern of smaller over-prediction for interface earthquakes and greater over-prediction for slab earthquakes. Although there are relatively few studies of stress parameter for interface ruptures, this result is consistent with studies of global slab ruptures which have identified high stress parameter for slab events (e.g. Chhangte et al., 2021; García et al., 2004; Takeo et al., 1993) and with the equivalent point-source GMM of Hassani and Atkinson (2021) for Japan, for example, which uses smaller values of stress parameter for interface than for slab earthquakes. The results support the implementation of a larger stress parameter for slab and a slightly smaller stress parameter for interface ground-motion simulations relative to crustal-based simulations.

The interpretation of analysis results for rupture velocity is slightly more nuanced since the shear wave velocities for crustal models from the SRCMOD database are systematically low. Therefore, both absolute rupture velocity (km/s) and (relative) rupture velocity ratio were considered, the results for which are shown in Figure 6. For interface earthquakes, the analysis of both absolute and relative rupture velocities indicates over-prediction by the reference crustal-based value compared with the SRCMOD database, that is, supporting slower interface rupture velocities, which is consistent with studies of interface ruptures (e.g. Macias et al., 2008; Mikumo et al., 1998).

For slab earthquakes, the results for relative and absolute rupture velocities show the opposite behavior and indicate that although the rupture velocity (measured in km/s) for slab ruptures may be greater—which are consistent with studies of slab ruptures (e.g. Liu et al., 2020; Takeo et al., 1993)—they may actually have a slower velocity relative to the shear wave velocity of the (deeper) host rock than do crustal earthquakes (i.e. $V_R$ is greater, but $V_R/V_S$ may not be).

Studies of large-magnitude subduction ground-motion simulations

Although studies of subduction earthquakes have focused on a few large-magnitude earthquakes—which may have limited similarity to small magnitudes—they provide insights and form a compliment to the analyses of empirical GMMs and SRCMOD database in the previous sections.

Ground-motion simulation studies of historical interface earthquakes include the 2003 Tokachi-Oki $M_w$ 8.3 (e.g. Wirth et al., 2017), 2010 Maule $M_w$ 8.8 (e.g. Frankel, 2017; Lay et al., 2010), and 2011 Tohoku $M_w$ 9.0 (e.g. Frankel, 2013; Kurahashi and Irikura, 2011) earthquakes, among others. Such studies have determined that relatively low background rupture velocities of 2.5–2.8 km/s well-represent the character of observed velocity pulses for Tohoku and Maule (e.g. Frankel, 2017); however, Tokachi-Oki was determined to have a rupture velocity at or above the local shear velocity (e.g. Wirth et al., 2017). Depth-dependent models (e.g. Lay et al., 2012) for within-rupture variations of slip and stress parameter, with decreased slip and increased stress parameter at depth, have been supported by simulation studies (e.g. Frankel, 2017) and such models may extend well to between-rupture property variations.

Source studies of slab ruptures have focused on events which were notable either for having large magnitudes or for being particularly deep. Source characterization of shallow slab ruptures by Asano et al. (2003) using stochastic ground-motion simulations identified strong ground-motion generating asperity regions with areas that tended to decrease with depth, but for which stress parameter increased with depth. Wei et al. (2013) analyzed the 2013 $M_w$ 8.3 Sea of Okhotsk earthquake and found that the rupture involved multiple rapid complementary sub-events over a compact fault area, suggesting highly efficient rupture triggering with characteristically short rupture duration. At great depths, the postulated combination of large confining stresses, heat radiating from the core, and macro-stresses of the descending slab has been found to produce ruptures with characteristically high static stress drop and short rupture durations (Chhangte et al., 2021; García et al., 2004; Takeo et al., 1993). Iwata and Asano (2011) determined that for a given magnitude, slab earthquakes have characteristically small rupture and asperity areas with less slip, which implies greater stress parameter. Ichinose et al. (2006) similarly observed very small asperity areas for deep slab earthquakes and greater stress parameter.

Subduction-specific simulation parameters

We developed subduction-specific modifications for stress parameter, rupture velocity, and anelastic attenuation based on synthesis of the observations made and fitting to theoretical considerations supported by the evidence from (1) empirical GMMs, (2) SRCMOD database, and (3) insights from other studies of subduction earthquakes discussed in the prior sections. The specific model coefficients were developed using an iterative approach—roughly equal consideration given to items 1, 2, and 3—with examination of prediction residuals and their dependence on phenomenological parameters, while also attempting to avoid over-fitting limited data with relatively simple models. The adopted models are presented sequentially below and are then compared and contrasted with the available information sources. The adjustments introduced by the subduction-specific models, shown in Figures 5 and 6, and the Electronic Supplement, demonstrate improved alignment of subduction-specific simulation models with global rupture models from the SRCMOD database.

Stress parameter

For stress parameter, $\mathrm{\Delta }\mathrm{\sigma }$ , models with linear depth dependence were determined to well-represent both interface and slab earthquakes.

$$\mathrm{\Delta }{\mathrm{\sigma }}_{\mathrm{if}}=10+1.25D,\mathrm{\Delta }{\mathrm{\sigma }}_{\mathrm{if}}\le 85\mathrm{bar}$$

(6)

$$\mathrm{\Delta }{\mathrm{\sigma }}_{\mathrm{slab}}=50+2D,\mathrm{\Delta }{\mathrm{\sigma }}_{\mathrm{slab}}\le 350\mathrm{bar}$$

(7)

where $D$ is centroid depth in km, and $\mathrm{\Delta }\mathrm{\sigma }$ is in units of bar. These stress parameter models for subduction earthquakes and a histogram of the depths of validation events considered are shown in Figure 7; also included for reference is the crustal-based simulation model (e.g. Lee et al., 2022) which uses $\mathrm{\Delta }{\mathrm{\sigma }}_{\mathrm{crustal}}=50\mathrm{bar}$ . Due to limited validation data for deeper slab events, and observed plateauing of high-frequency pSAs for deep slab events (Figure 3), $\mathrm{\Delta }{\mathrm{\sigma }}_{\mathrm{slab}}$ is limited to $350$ bar, which corresponds to a hypocentre depth of $150$ km. Similarly, $\mathrm{\Delta }{\mathrm{\sigma }}_{\mathrm{if}}$ is limited to $85$ bar which corresponds to $60$ km—deeper than any interface earthquake in the validation dataset.

Figure 7.

(a) Simulation models of stress parameter for crustal, interface, and slab earthquakes, and (b) depth distributions for interface and slab events used for validation shown as a shaded histogram. The models are not applicable for all depths shown, particularly those for crustal and interface, but are extended (where dotted) to facilitate comparison between the three types of earthquakes. (a) Stress parameters from selected studies of historical events: T93 (Takeo et al., 1993), AS97 (Atkinson and Silva, 1997), RB03 (Roumelioti and Beresnev, 2003), G04 (García et al., 2004), M08 (Macias et al., 2008), and L20 (Liu et al., 2020), and from equivalent point-source models: YA15 (Yenier and Atkinson, 2015) and HA21 (Hassani and Atkinson, 2021) are shown for comparison. (a) Simulation models for stress parameter. (b) Depth distribution of events.

When applied to adjust estimated rise times for the SRCMOD rupture models, the subduction-specific models produce a more consistent magnitude-dependence relationship of between-event residuals across tectonic classifications (Figure 5b). The depth dependence of stress parameter for the subduction-specific models (Figure 7a) exhibits similar trends to the ratios of the high-frequency pSAs from empirical GMMs (Figure 3c), the stress parameter models in the equivalent point-source models of Hassani and Atkinson (2021), and matches the high stress parameters which have been widely reported for deep slab ruptures, for which some selected values are shown in Figure 7a (e.g. Iwata and Asano, 2011). While the exact formulation of stress parameter varies within the literature (Atkinson and Beresnev, 1997), the similarity of the subduction-specific models to the range of values from studies of historical ruptures provides confidence in their suitability.

Rupture velocity

Constant values of rupture velocity ratio (relative to shear wave velocity at the source), $V_R/V_S$ , were used for interface and slab earthquakes, these were set at.

$${V_R}_{\mathrm{if}}/V_S=0.75$$

(8)

$${V_R}_{\mathrm{slab}}/V_S=0.9$$

(9)

where $V_S$ is the shear wave velocity at the source—crustal-based simulations use a relative rupture velocity of 0.8. Velocity-strengthening fault zones, which are modeled for crustal earthquakes through reduced rupture velocity in the shallow and deep crust (e.g. Graves and Pitarka, 2016), were not implemented for interface and slab earthquakes.

Figure 6c and d indicates that subduction earthquakes have low rupture velocity ratios compared with crustal earthquakes; comparison of the shear wave velocities in the SRCMOD database with the 1D velocity model used in HF simulations and the empirical relations of Brocher (2005) (see the Electronic Supplement) indicates that shear wave velocities for crustal models from the SRCMOD database are systematically low and explain this observation. Therefore, greater importance was assigned to findings from studies of the phenomenological rupture processes for these events and the observed trends in rupture velocities in the SRCMOD database. To achieve a desired decrease in rupture velocity for interface earthquakes, the ratio to shear wave velocity at the source was reduced slightly from the value used for crustal earthquakes. Similarly, for slab earthquakes, a larger ratio of shear wave velocity at the source is specified to achieve faster rupture velocities. As shown in Figure 6a and b, the subduction-specific models effectively remove most of the differences in between-event residuals between tectonic classifications indicating that the approximate models are reasonable and do not over-correct. Figure 6c and d indicates that the subduction-specific models increase the between-event $V_R/V_S$ residuals for slab ruptures (in contrast to the reductions seen for $V_R$ ); however, this is primarily caused by differences between the shear wave velocity models for the simulations and the SRCMOD database, as examined in the Electronic Supplement. Therefore, the reductions in $V_R$ residuals were more meaningful for assessing subduction-specific models for rupture velocity. The decreased rupture velocity for interface and increased rupture velocity for slab ruptures are consistent with studies of these events which have observed these behaviors (e.g. Chhangte et al., 2021; Frankel, 2017).

Anelastic attenuation

The HF simulation component uses a 1D velocity model which has been developed for NZ. Although this model has been applied successfully for shallow crustal earthquake ground-motion simulation throughout NZ, conceptually it is a poor representation of the deeper and laterally varying crustal structures which are important for subduction earthquakes, such as forearc/backarc differences in anelastic attenuation and the greater anelastic attenuation in the Taupō Volcanic Zone (TVZ). Although such features warrant treatment in the HF simulation component through path-specific 1D velocity models, the computation of seismic moment in the HF component simulation has been calibrated based on the velocity structure in the 1D velocity model and it is not possible to modify this velocity structure without significant changes to the simulation computations. Furthermore, the 1D model $V_S$ values are mainly used for impedance, not attenuation (Baker et al., 2021, Chapter 5), and therefore, a direct adjustment to $Q$ is appropriate in order to isolate attenuation effects of interest here. Several possible approaches were considered to incorporate such effects: using a 3D velocity model (Thomson et al., 2020) for HF, applying depth- and azimuth-based correction factors to the 1D velocity model, or developing 1D velocity models with path-specific velocity and rock quality layers; however, because the HF component is based on simplified physics, the exact 1D velocity structure has limited importance and a pragmatic approach was taken in which the velocity structure is preserved, but the rock quality factors, $Q_P$ and $Q_S$ , are adjusted.

A ray-tracing approach, similar to the cell-based attenuation approach of Dawood and Rodriguez-Marek (2013), was used to adjust the rock quality factors based on the ratio of path attenuation, $t^{*}$ , in the 1D and 3D velocity models—computed along a simplified straight-line source-to-site ray path. For the 1D velocity model, the path attenuation is summed based on the travel time-averaged rock quality across the layers. For the 3D velocity model, the straight-line approximation of the travel path is discretized to ten segments and the travel time-averaged quality factor for each segment is computed based on the properties of the nearest grid point in the 3D velocity model. Separate scale factors are computed for $Q_P$ and $Q_S$ with $t^{*}$ computed as per Eberhart-Phillips et al. (2015):

$$t^{*}={\sum }^{\frac{x}{\mathrm{VQ}}}$$

(10)

where $x$ is the travel distance along each discretized ray path, $V$ is the velocity ( $V_P$ or $V_S$ ), and $Q$ is the quality factor ( $Q_P$ or $Q_S$ ). The quality factor “scaling factor” is then determined from:

$$Q_{\mathrm{scale}}=\frac{t^{*}\left(1D\right)}{t^{*}\left(3D\right)}$$

(11)

Because the velocity is not adjusted, differences in both the velocity and rock quality between the adopted 1D velocity model and the 3D velocity model are mapped to a single adjustment of the rock quality factor in the adjusted 1D velocity models. The adjusted 1D velocity models were found to produce anelastic attenuation which exhibited similarities to isoseismal intensity maps of historical earthquakes collated by Downes (1996); however, as shown in the Electronic Supplement, for most cases the salient features—offset from the epicenter and elliptical shape—are significantly less pronounced. Spatial variation of anelastic attenuation is explicitly included in the structure of the 3D velocity model and 3D wave propagation (Graves, 1996); therefore, corresponding adjustments to the LF component simulation are not needed.

Simulation validation results

The subduction-specific models were implemented in ground-motion simulations and compared with observed ground motions and simulations using crustal-based models. Prediction residuals from Lee et al. (2022) for 5218 ground motions (Table 1) are also included to assess prediction performance relative to crustal ground-motion simulations. Predictions are compared using the previously described mixed-effects regression approach. Predictions for pSAs and additional selected IMs: peak ground acceleration, $\mathrm{PGA}$ ; peak ground velocity, $\mathrm{PGV}$ ; cumulative absolute velocity, $\mathrm{CAV}$ ; and significant durations encompassing 5%–75% and 5%–95% of the Arias intensity, $D_{s575}$ and $D_{s595}$ , respectively, are considered. RotD50 (Boore et al., 2006) IMs were used for observed, simulated, and empirically predicted ground motions. The validation presented in this study is appropriate for Fourier amplitudes at frequencies above approximately $0.33$ Hz or pSAs up to periods of approximately $3$ s beyond which the majority of observed ground-motion records have signal-to-noise ratios below 2. Each observed ground-motion record was used only at periods where the signal-to-noise ratio exceeds 2; therefore, the validation results shown at periods above $3$ s are based on fewer records. Approximately half as many records are used at T = 5 s compared with T = 3 s, as presented in the Electronic Supplement. Additional figures showing the validation results are presented in the Electronic Supplement.

Examination of effects for archetype earthquakes

To elucidate the general effects of the subduction-specific models, the simulated ground motions for well-recorded archetype interface and slab earthquakes from the validation dataset are first presented.

Predictions from simulations conducted with crustal-based and interface-specific simulation models for a single $15$ km deep $M_w$ 4.9 interface earthquake are shown in Figure 8. The observed trends indicate that both short- and long-period pSAs are generally reduced relative to those simulations based on crustal parameters. For large rupture distance, $R_{\mathrm{rup}}$ , the $Q$ adjustment to the 1D velocity model significantly increases the simulated short-period pSAs, and produces larger predictions of short-period pSAs at distances greater than $100$ km relative to the crustal-based models. These effects are consistent with observed characteristics of empirical GMMs for interface earthquakes which exhibit similar trends at short- and long-period pSAs (e.g. Figure 2).

Figure 8.

Observed and predicted pSAs for periods of (a) $0.1$ s and (b) $2.0$ s from simulations with crustal-based and interface-specific models shown by rupture distance. Results for a single $M_w$ 4.9 interface earthquake (ID: 2016p861632) located $15$ km deep at 42.3° S, 173.7°E are shown. Lines of best fit from linear regressions of the results are shown to illustrate the relative differences between the predictions. Observed pSAs outside the simulation velocity model domain are shown with small markers and were not included in the line of best fit. (a) pSA ( $0.1$ s) by rupture distance. (b) pSA ( $2.0$ s) by rupture distance.

Similarly, simulated ground motions using crustal-based and slab-specific simulation models for a single $M_w$ 4.3 slab earthquake are shown in Figure 9. Both short- and long-period pSAs are increased by the implementation of the slab-specific models; however, short-period pSAs are increased much more significantly than long periods, especially at greater rupture distances. Taken together, these observations indicate that the subduction-specific models for stress parameter and rupture velocity have a greater effect for slab ground motions than the adjustments to the 1D velocity model.

Figure 9.

Observed and predicted pSAs for periods of (a) $0.1$ s and (b) $2.0$ s from simulations with crustal-based and slab-specific models shown by rupture distance. Results for a single $M_w$ 4.3 slab earthquake (2020p262796) located $42$ km deep at 39.8° S, 176.2°E are shown. Lines of best fit from linear regressions of the results are shown to illustrate the relative differences between the predictions. Observed pSAs outside the simulation velocity model domain are shown with small markers and were not included in the line of best fit. (a) pSA ( $0.1$ s) by rupture distance. (b) pSA ( $2.0$ s) by rupture distance.

Observed systematic effects

Figures 10 and 11 provide a summary of the mixed-effect regression analysis results for interface and slab earthquakes, respectively. For interface ground motions, Figure 10a illustrates the subduction predictions have comparable magnitude of bias, a, as simulations for crustal ground motions (Lee et al., 2022), and the bias is reduced at short-period pSAs through implementation of the subduction-specific models. Comparing the effects of the different subduction-specific modifications, the bias increases through the reductions of stress parameter and rupture velocity but is then reduced through scaling of quality factor—this supports the suite-based-approach with model-specific adjustments for subduction ground motions, as opposed to individual model refinements. Because the majority of sites are not located within the backarc, the average effect of the path-specific modifications to the 1D velocity models is increased values of rock quality, $Q_P$ and $Q_S$ , which results in greater IM predictions in the HF simulation component; this was also true for slab ground motions which were mainly recorded at a common set of sites.

Figure 10.

Mixed-effect regression results for interface earthquake ground-motion simulations with crustal-based and subduction-specific simulation models. The simulations labeled “Interface $\mathrm{\Delta }\mathrm{\sigma }$ ,”“Interface $\mathrm{\Delta }\mathrm{\sigma }$ and $V_R$ ,” and “All interface parameters” represent the incremental cumulative implementation of subduction-specific models for stress parameter, rupture velocity, and anelastic attenuation, respectively. Crustal-based simulations for small-magnitude crustal earthquakes by Lee et al. (2022) are shown for reference. For (d) depth dependence of $\delta B_e$ , the solid line is the hypocentre depth coefficient from a generalized multivariate linear regression of $\delta B_e$ with hypocentre depth and magnitude considered as independent variables. Gray shading is used to indicate periods where some records have signal-to-noise ratios below 2 and therefore fewer ground motions were used in the mixed-effects regression. (a) Bias, $a$ . (b) Total SD, $\sigma$ . (c) Between-event SD, $\tau$ . (d) Depth dependence of the between-event residual, $\delta {B_e}^{\prime }\left(\mathrm{Depth}\right)$ . (e) Systematic site-to-site SD, ${\varphi }_{S2S}$ . (f) Single-station within-event SD, ${\varphi }_{\mathrm{SS}}$ .

Figure 11.

Mixed-effect regression results for slab earthquake ground-motion simulations with crustal-based and subduction-specific simulation models. The simulations labeled “Slab $\mathrm{\Delta }\mathrm{\sigma }$ ,”“Slab $\mathrm{\Delta }\mathrm{\sigma }$ and $V_R$ ,” and “All slab parameters” represent the incremental cumulative implementation of subduction-specific models for stress parameter, rupture velocity, and anelastic attenuation, respectively. Crustal-based simulations for small-magnitude crustal earthquakes by Lee et al. (2022) are shown for reference. For (d) depth dependence of $\delta B_e$ , the solid line is the hypocentre depth coefficient from a generalized multivariate linear regression of $\delta B_e$ with hypocentre depth and magnitude considered as independent variables. Gray shading is used to indicate periods where some records have signal-to-noise ratios below 2 and therefore fewer ground motions were used in the mixed-effects regression. (a) Bias, $a$ . (b) Total SD, $\sigma$ . (c) Between-event SD, $\tau$ . (d) Depth dependence of the between-event residual, $\delta {B_e}^{\prime }\left(\mathrm{Depth}\right)$ . (e) Systematic site-to-site SD, ${\varphi }_{S2S}$ . (f) Single-station within-event SD, ${\varphi }_{\mathrm{SS}}$ .

For slab ground motions, Figure 11a illustrates the initially large bias of simulations with crustal-based parameter models indicating that the crustal-based models significantly under-predict short-period pSAs. The bias decreases, corresponding to increased predicted pSAs, at all periods due to implementation of each slab-specific simulation model. In particular, the slab-specific model for stress parameter significantly reduces overall bias, especially at short-period pSAs. There is a slight over-prediction at short periods and slight under-prediction at long-period pSAs once the suite of slab-specific models are implemented. Overall, the slab-specific models significantly reduce bias and their adoption is further supported by analysis of the relative change of these models from the crustal-based models which is examined in greater detail in the Electronic Supplement.

For interface and slab ground motions, total SD, $\sigma$ , was slightly reduced at very short periods, that is, less than $0.1$ s, and largely unchanged above periods of about $0.1$ s (Figures 10b and 11b), compared to simulated subduction ground motions with crustal parameters. The SD for both interface and slab ground motions is comparable to that for crustal ground-motion simulations. However, there is a distinct upward trend at periods greater than $3$ s (Figures 10 and 11) which the crustal ground motions (Lee et al., 2022) do not exhibit. This is because the subduction ground motions in this study were only required to have a usable period range up to $3$ s due to reasons of record scarcity, as illustrated in the Electronic Supplement, while Lee et al. (2022) was able to select crustal ground motions with an extended usable period range. For slab ground motions in particular, the SD increases significantly at long periods, that is, greater than $3.0$ s which is a suspected consequence of decreasing signal-to-noise ratio, and therefore a reduction of the number of available records at long periods.

For both interface and slab ground motions, the between-event SD, $\tau$ , reflects the reductions in depth dependence of the between-event residual with a significant associated reduction at short periods (Figures 10c and 11c, respectively). For slab ground motions, the slab-specific parameter models for stress parameter slightly increase the between-event SD at moderate periods. For both interface and slab ground motions, the values are comparable in magnitude to crustal simulation results.

For interface and slab ground motions, an investigation of the depth dependence of the between-event residual, $\delta {B_e}^{\prime }\left(\mathrm{Depth}\right)$ , is presented in Figures 10d and 11d, respectively. With the crustal-based parameter models, both interface and slab between-event residuals exhibit significant depth dependence, especially at short periods, and with stronger dependence for slab earthquakes. Both the interface and slab models for stress parameter reduce the slope of between-event residual with depth, especially at short-period pSAs which are affected by the depth-dependent modifications to stress parameter made for both interface and slab earthquakes. These modifications to stress parameter for slab earthquakes had greater depth dependence and a greater corresponding reduction in depth dependence of the between-event residual for slab ground motions was observed. The results indicate the slab-specific parameterization of stress parameter are effective at reducing depth dependence at all periods below the transition frequency ( $1$ s or $1$ Hz)—depth dependence above this period is not apparent, even with the crustal-based models. Magnitude dependence of the between-event residual was examined for interface and slab ground motions and was not found to be significant for the magnitude range considered in this study ( $M_w$ 3.5–5).

The systematic site-to-site SD, ${\varphi }_{S2S}$ , and single-station within-event SD, ${\varphi }_{\mathrm{SS}}$ , were also compared for interface (Figure 10e and f) and slab (Figure 10e and f) ground motions. The results indicate that the variability in residuals between observation and prediction due to systematic site effects for both interface and slab ground motions was unaffected by the implementation of the new models. For interface ground motions, the systematic site-to-site SD is comparable to the crustal-based simulations for crustal ground motions at short periods, and lower at long periods. For slab ground motions, the systematic site-to-site SD is comparable to the crustal-based simulations for crustal ground motions at moderate periods, and larger at short periods; the systematic site-to-site SD for slab ground motions increases at long periods for reasons discussed earlier for total SD. Single-station within-event SD, which measures the variability of the portion of the residual not partitioned to between- and within-event effects, was also unaffected by the implementation of interface-specific simulation models.

Between-event residuals

Additional analysis was conducted to examine the relationship of between-event residuals, $\delta B_e$ , with earthquake source properties. Figures 12 and 13 show the depth and magnitude dependence of residuals for spectral accelerations at $0.1$ s and $3.0$ s, corresponding to high- and low-frequency responses for interface and slab events, respectively. Results are shown as “ $a+\delta B_e$ ,” that is, including the constant bias term, in order to indicate whether the residuals correspond to over- or under-prediction. Further examination of the depth and magnitude dependence of the between-event residuals for both interface and slab ground-motion simulations is included in the Electronic Supplement.

Figure 12.

Comparison of the (a, c) depth and (b, d) magnitude dependence of the between-event residual, $\delta B_e$ , for (a, b) pSA ( $0.1$ s) and (c, d) pSA ( $3.0$ s) from mixed-effect regression for interface earthquake ground-motion simulations with crustal-based and subduction-specific simulation models. Crustal-based simulations for small-magnitude crustal earthquakes by Lee et al. (2022) are shown for reference. The solid lines are the best fits from generalized multivariate linear regressions of $a+\delta B_e$ . (a) Depth dependence, pSA ( $0.1$ s) (b) Magnitude dependence, pSA ( $0.1$ s) (c) Depth dependence, pSA ( $3.0$ s). (d) Magnitude dependence, pSA ( $3.0$ s).

Figure 13.

Comparison of the (a, c) depth and (b, d) magnitude dependence of the between-event residual, $\delta B_e$ , for (a, b) pSA ( $0.1$ s) and (c, d) pSA ( $3.0$ s) from mixed-effect regression for slab earthquake ground-motion simulations with crustal-based and subduction-specific simulation models.

For both interface and slab ground motions, the implementation of the subduction-specific simulation models reduces $a+\delta B_e$ at all periods, with significant reduction in the depth dependence of $\delta B_e$ , especially for short periods but also for long periods, due to the implementation of the depth-dependent model for stress parameter.

Magnitude dependence of between-event residuals for interface and slab ground motions was also examined, as shown in Figures 12 and 13. Because this study was concerned with small-magnitude events, which are sufficiently represented with point-source rupture models, the subduction-specific models did not include explicit treatment for magnitude-dependent scaling. Therefore, the magnitude dependence of the $\delta B_e$ residuals was not explicitly addressed through magnitude-dependent models for stress parameter or rupture velocity. However, for slab ground motions at HF, the inclusion of a depth-dependent model for stress parameter introduces magnitude dependence in the residuals due to the correlation of depth and magnitude for slab events in the validation dataset. This effect was also observed, although to a smaller degree, for $\delta B_e$ residuals for interface ground motions.

Examination of geospatial trends

Figures 14 and 15 illustrate the geospatial variation of the between-event and systematic site-to-site residuals with the subduction-specific simulation models for pSA ( $0.1$ s) and pSA ( $3.0$ s) for interface and slab ground motions, respectively. For between-event residuals of interface ground motions in Figure 14, systematic and coherent patterns are not apparent which indicates that large-scale regional misfit is not a problem and that the models provide an acceptable representation of seismological factors. This is also generally true for slab ground motions; however, the between-event residuals indicate that there is a slight tendency toward short-period over-prediction for events located in the TVZ. Although this slight short-period over-prediction manifests in the between-event residual, $\delta B_e$ , it could indicate increased path attenuation which is not adequately modeled at short periods. For both interface and slab ground motions, there are not enough events in the Puysegur Subduction Zone to make conclusions regarding regional differences in the residuals observed for the Hikurangi and Puysegur subduction zones.

Figure 14.

Geospatial variation of (a, b) between-event residual, $\delta B_e$ , and (c, d) site-to-site residual, $\delta S2S$ , for (a, c) pSA ( $0.1$ s) and (b, d) pSA ( $3.0$ s) for interface earthquake ground-motion simulations. The extents of the Taup $\overline{o}$ Volcanic Zone (TVZ), and depth contours of the Hikurangi and Puysegur subduction zones are shown for reference. All earthquakes are plotted but some earthquakes are densely clustered and hence are overlapping and not visible. (a) $\delta B_e$ for pSA ( $0.1$ s). (b) $\delta B_e$ for pSA ( $3.0$ s). (c) $\delta S2S$ for pSA ( $0.1$ s). (d) $\delta S2S$ for pSA ( $3.0$ s).

Figure 15.

Geospatial variation of (a, b) between-event residual, $\delta B_e$ , and (c, d) site-to-site residual, $\delta S2S$ , for (a, c) pSA ( $0.1$ s) and (b, d) pSA ( $3.0$ s) for slab earthquake ground-motion simulations. The extents of the Taup $\overline{o}$ Volcanic Zone (TVZ), and depth contours of the Hikurangi and Puysegur subduction zones are shown for reference. All earthquakes are plotted but some earthquakes are densely clustered and hence are overlapping and not visible. (a) $\delta B_e$ for pSA ( $0.1$ s). (b) $\delta B_e$ for pSA ( $3.0$ s). (c) $\delta S2S$ for pSA ( $0.1$ s). (d) $\delta S2S$ for pSA ( $3.0$ s).

For systematic site-to-site residuals of interface ground motions, there are not apparent macro trends—such as toward over-prediction in the backarc region or TVZ, or under-prediction in the forearc region. This generally supports that the adjustments to the 1D velocity models to account for 3D variations in anelastic attenuation are appropriate. Similarly, for slab ground motions, the systematic site-to-site residuals do not strongly indicate over-prediction in the backarc region for either the Hikurangi or Puysegur subduction zones. However, as was the case for between-event residuals, there is over-prediction for sites located within or adjacent to the TVZ, and very large over-prediction for a single site (TOZ), at short periods for travel paths through the backarc which may indicate greater path attenuation than accounted for in this region. Further examination of these geospatial trends in systematic site-to-site residuals focused on assessing the effects of targeted path-specific adjustments to the 1D velocity models are included in the Electronic Supplement.

Systematic site-to-site residuals

A significant area of targeted improvement in this study was reduction of systematic geospatial variation in the residuals resulting from deep and laterally varying features in the velocity structure of NZ. To accomplish this, path-specific adjustments to the 1D velocity model were made for the simulation of each ground motion. The change in site-to-site residuals due to these adjustments are shown in Figure 16 for both interface and slab earthquake ground motions. Although the change in site-to-site residuals has geospatial variation and residuals actually increase for certain sites, the adjustments have the general effect of reducing the residuals. The differences are greatest for short-period pSAs, which are more sensitive to changes in anelastic attenuation. The path-specific adjustments were based on a simplistic straight-line ray path approach. For deep slab earthquakes, deviation of curved ray paths from this straight-line approximation may explain why the reduction of site-to-site residuals is less than anticipated and less than for relatively shallow interface earthquakes, which may be affected less by the straight-line approximation.

Figure 16.

Spatial variation of change to site-to-site residual, $\mathrm{\Delta }\mathrm{ABS}\left(\delta S2S_S\right)$ , of (a, c) pSA ( $0.1$ s) and(b, d) pSA ( $3.0$ s) from ground-motion simulation of (a, b) interface and (c, d) slab earthquakes with the crustal 1D velocity model (VM) and with path-adjusted velocity models. Negative values correspond to a reduction in absolute value of the site-to-site residual. (a) Interface, pSA ( $0.1$ s). (b) Interface, pSA ( $3.0$ s). (c) Slab, pSA ( $0.1$ s). (d) Slab, pSA ( $3.0$ s).

Site-to-site residuals, $\delta S2S_s$ , determined for different tectonic classes of earthquakes, were compared with one another for the same site (refer to the Electronic Supplement) to understand whether similarities or differences exist between tectonic classes. Most sites exhibit very similar behavior for all tectonic classes considered, which indicates that the dominant factors affecting the site-to-site residuals are features related to the site. Given that the majority of sites demonstrate similar residuals for all three tectonic classes, these residuals can be most directly reduced through improved site characterization and treatment of site effects with targeted improvements for poor-performing sites. This finding underscores the importance of site-specific factors relative to regional differences, which were not given explicit treatment in this study (Kuncar et al., 2025).

Discussion

By comparing the ratios of pSA predictions from the crustal-based and subduction-specific models for interface and slab ground motions, as shown in Figure 17, further confidence in the credibility of the adopted subduction-specific models was gained. The relative period-to-period trends for each tectonic type in Figure 17 are consistent with those in Figure 2, as are the relative trends of interface-to-slab ground motions. Although some of the specific trends with rupture distance and magnitude differ between the empirical- and simulation-based ratios, such as for interface pSA at $2.0$ s, this can be explained in part by the differing rupture distances and magnitudes considered by each—empirical models are generally applicable for moderate-to-large magnitudes and rupture distances up to $300$ km.

Figure 17.

Comparison of the smoothed ratios of pSAs predicted for interface and slab earthquake ground motions using the subduction-specific and crustal-based simulation models by (a) rupture distance, $R_{\mathrm{rup}}$ , and (b) magnitude, $M_w$ . For each ground motion, the ratio of the simulation predictions with subduction parameter models to crustal-based models was computed which resulted in $367$ and $975$ points for the interface and slab datasets, respectively. A moving average was then applied to produce a smoothed curve. (a) Predictions relative to crustal models vs $R_{\mathrm{rup}}$ . (b) Predictions relative to crustal models vs $M_w$ .

While the parametric relationships developed for rupture velocity and stress parameter are approximate, this approach was deliberate to avoid over-fitting the limited data available; further work could be done to refine these relationships on an expanded catalogue of global ground motions. Site variability effects were previously investigated for crustal ground-motion simulations (e.g. Lee et al., 2022) for which there is relatively plentiful validation data; findings were inconclusive due to significant scatter and biases from other site effects which obscured any possible trends. Such phenomena were not explicitly considered in this study of subduction earthquake ground motions for which there are insufficient data to conduct a rigorous investigation.

Several areas of study which future research could address were identified. The subduction earthquake ground-motion simulations exhibit partial support of persistent regional effects related to the TVZ which may be reduced through improved knowledge of the deeper velocity structure and refined regional scale velocity models. The adjustments to anelastic attenuation were based on a straight-line ray path assumption. Future work could implement an adjustment with more accurate treatment of travel paths for deep slab earthquakes.

The subduction-specific models developed in this work mainly affect the high-frequency component of the ground-motion simulations above $1$ Hz; therefore, the improvements from crustal-based simulations in the low-frequency range are limited. Future studies should focus on making improvements to the low-frequency component of the hybrid broadband approach. This study considered small-magnitude events which have ample recordings and thus presented a logical domain for validation. The validation should be extended to larger magnitude events, which dominate seismic hazard, and necessitate the use of finite-fault source models for which low-frequency ground motions are sensitive to rupture geometry and slip variability, among other complexities. This should include implementing 3D models of rock quality in the LF simulations, instead of simply using parametric relationships that are a function of velocity, especially as this component will extend to higher frequencies as computation capabilities improve.

This study considered earthquakes along the Hikurangi and Puysegur subduction zones together in a combined regression analysis because too few subduction earthquakes were well recorded from the Puysegur Subduction Zone to facilitate a separate validation. Future work could leverage additional instrumentation and the passage of time to do an independent validation for each subduction zone and possibly extend the validation to other global subduction zones. Similarly, residuals could be partitioned to quantify differences between regions (e.g. Hikurangi and Puysegur) by computing a systematic region-to-region residual in a combined regression analysis.

It is acknowledged that the applicability of the simulation models validated in this study may vary between different subduction systems due to regional differences in tectonic and geologic characteristics. Moreover, the transferability of these modifications may be influenced by factors such as the selection of velocity models, average rupture velocities, and the source generation methodology and thus further validation studies are warranted to evaluate their generalizability across different simulation frameworks.

Uncertainty was not a primary focus of this study; however, the ground-motion dataset considered was limited in size relative to many validation studies of small-magnitude crustal earthquakes (e.g. Lee et al., 2022). Therefore, practitioners should note that although the simulations are able to predict subduction ground motions well, the predictions likely suffer from greater uncertainty than crustal ground-motion simulations.

Conclusion

In this study, subduction-specific models of stress parameter, rupture velocity, and anelastic attenuation were implemented within the hybrid broadband ground-motion simulation approach and validated using small-magnitude ( $M_w$ 3.5–5) subduction earthquake ground-motion records in NZ. Only observed ground-motion records with two high-quality horizontal components and maximum usable vibration periods of T ≥ 3 s were considered. The subduction-specific models were motivated by observed differences in the characteristics of interface and slab finite-fault rupture models and the treatment of subduction earthquakes by empirical GMMs. Focus was on adjustments that have physical basis in addition to being constrained by empirical calibration against validation data.

Relative to crustal-based simulation models, a smaller stress parameter with moderate depth dependence and a slightly slower rupture velocity were used for interface earthquakes. For slab earthquakes, a larger stress parameter with significant depth dependence and a faster rupture velocity were implemented. Ray-path-based adjustments were made to the rock quality factors in the 1D velocity models to account for deep and laterally varying geophysical structures in the 3D velocity model which are important for subduction earthquakes.

The results indicate the depth-dependent parameterization of stress parameter was effective at reducing depth dependence of residuals and that the path-based adjustments to the 1D velocity models slightly reduced geospatial trends of residuals. When applied with the subduction-specific models for rupture velocity, the parameter models generally reduced bias of the prediction residuals compared with crustal-based parameters, especially for slab ground motions at short periods. The validation of simulation predictions done with subduction-specific models achieved performance for subduction earthquakes which is similar to that of the crustal-based models for crustal earthquakes in NZ.

We believe this work has advanced the field by (1) examining small-magnitude ( $M_w$ 3.5–5) subduction event validation and (2) proposing subduction-specific simulation parameter models. While the simulation accuracy and precision are comparable to studies on crustal events, there clearly remains significant scope to further refine the subduction-specific models. This would include consideration of insights from moderate- and large-magnitude events, events beyond NZ, and other theoretical insights.

Supplemental Material

sj-zip-1-eqs-10.1177_87552930251353816 – Supplemental material for Hybrid broadband ground-motion simulation validation of small-magnitude subduction earthquakes in New Zealand

Supplemental material, sj-zip-1-eqs-10.1177_87552930251353816 for Hybrid broadband ground-motion simulation validation of small-magnitude subduction earthquakes in New Zealand by Michael R Dupuis, Robin L Lee and Brendon A Bradley in Earthquake Spectra