Host-region parameters for an adjustable model for crustal earthquakes to facilitate the implementation of the backbone approach to building ground-motion logic trees in probabilistic seismic hazard analysis

Abstract

The backbone approach to constructing a ground-motion logic tree for probabilistic seismic hazard analysis (PSHA) can address shortcomings in the traditional approach of populating the branches with multiple existing, or potentially modified, ground-motion models (GMMs) by rendering more transparent the relationship between branch weights and the resulting distribution of predicted accelerations. To capture epistemic uncertainty in a tractable manner, there are benefits in building the logic tree through the application of successive adjustments for differences in source, path, and site characteristics between the host region of the selected backbone GMM and the target region for which the PSHA is being conducted. The implementation of this approach is facilitated by selecting a backbone GMM that is amenable to such host-to-target adjustments for individual source, path, and site characteristics. The NGA-West2 GMM of Chiou and Youngs (CY14) has been identified as a highly adaptable model for crustal seismicity that is well suited to such adjustments. Rather than using generic source, path, and site characteristics assumed appropriate for the host region, the final suite of adjusted GMMs for the target region will be better constrained if the host-region parameters are defined specifically on the basis of their compatibility with the CY14 backbone GMM. To this end, making use of a recently developed crustal shear-wave velocity profile consistent with CY14, we present an inversion of the model to estimate the key source and path parameters, namely the stress parameter and the anelastic attenuation. With these outputs, the effort in constructing a ground-motion logic tree for any PSHA dealing with crustal seismicity can be focused primarily on the estimation of the target-region characteristics and their associated uncertainties. The inversion procedure can also be adapted for any application in which different constraints might be relevant.

Introduction

A key objective in probabilistic seismic hazard analysis (PSHA) is to ensure both a representation of the best estimate of the site hazard, on the basis of the currently available data and models, and the adequate capture of the associated uncertainty, which arises from the inevitable incompleteness of the data. This objective is often expressed as capturing the center, body, and range of technically defensible interpretations, or CBR of TDI (United States Nuclear Regulatory Commission (USNRC), 2018). The tool ubiquitously deployed in PSHA to capture the CBR of TDI is the logic tree, with a node for each element of the PSHA inputs with weighted branches corresponding to the alternative models or parameter values that could represent each element. The CBR of TDI objective for ground-motion logic trees can be interpreted in terms of the resulting distribution of predicted ground-motion amplitudes. The traditional approach of populating the branches of the logic tree with published ground-motion models (GMMs) is not an effective way of meeting this objective because of the obscure relationship between the weights on the logic-tree branches and the resulting distribution of ground motions (Scherbaum et al., 2005). Consequently, in recent years it has become increasingly common in practice to populate the branches of a ground-motion logic tree with scaled versions of a single GMM (Atkinson et al., 2014). The scaling factors may vary as a function of both magnitude and distance, but in all cases the relationship between the weights on the branches and the resulting distribution—which can be broader than that defined by existing GMM predictions—becomes much more transparent. This general approach of populating a logic tree with scaled versions of a single GMM has been called the backbone approach (Bommer, 2012) and has been widely adopted using this terminology (Akkar et al., 2021; Douglas, 2018; Haendel et al., 2015; Kowsari and Ghasemi, 2021; Weatherill and Cotton, 2020; Weatherill et al., 2020). The approach has also been adapted to avoid similar pitfalls in constructing logic trees for site response analyses (Rodriguez-Marek et al., 2020). Once a backbone GMM is selected, the ideal way to construct the logic tree is through a series of nodes that correspond to adjustments for differences in source, path, and site parameters between the host region of the GMM and the target region of the PSHA. Adding branches for the uncertainty in the estimates of each of the parameters for which adjustments are made provides a physical and tractable basis for capturing epistemic uncertainty. To apply sequential host-to-target adjustments to the backbone GMM, it is desirable to select a model that is amenable to such targeted modifications. Bommer and Stafford (2020) previously identified the three key characteristics of an adaptable GMM to be as follows: (1) reliable host-region characterization; (2) isolated influence of individual source, path, and site characteristics; and (3) a functional form that reflects theoretical scaling of Fourier amplitude spectra. Chiou and Youngs (2014) was identified as the most adaptable of current GMMs for crustal earthquakes, having a functional form that closely mimics the theoretical scaling implicit in stochastic simulations and including terms that can be directly related to source (through the stress parameter, $\mathrm{\Delta }\sigma$ ), path (in terms of anelastic attenuation, as reflected in the quality factor, $Q$ ), and site (through explicit inclusion of $V_{S30}$ ) characteristics. However, while the Chiou and Youngs (2014) GMM—hereafter CY14—satisfies the second and third criteria very well, the challenge remains to obtain a suite of parameters that reflect the host-region source, path, and site characteristics. The focus and purpose of this article is to provide such a suite of host-region parameters for the CY14 model to facilitate implementation of the backbone GMM approach to constructing ground-motion logic trees for PSHA studies involving sources of crustal seismicity.

The article will now move to discuss the key differences between host-region vs model-specific Fourier parameters. Following that, the mathematical framework that is employed to invert the CY14 GMM is presented in detail. This includes the specification of the relevant system of equations, the optimization framework, the specification of the parameter space, explicit definition of the CY14 model to be inverted, and a presentation of the results. Finally, the article closes with a discussion of the results obtained and an interrogation of the model performance.

Host-region vs model-specific parameters

In the original formulation of the hybrid-empirical method (HEM) of Campbell (2003), in which adjustments for all host-to-target region differences are made in a single step, the host-region parameters were to be obtained from studies characterizing source and path features of the region for which the backbone GMM was developed. These source and path parameters could then be combined with a representative shear-wave velocity ( $V_S$ ) profile (Boore, 2016), anchored to a $V_{S30}$ value that is high enough to ensure mainly linear site response but low enough to be well-constrained by recording sites in this range, to complete the host-region characterization. Scherbaum et al. (2006) addressed the fact that many empirical GMMs combine data from several regions and, therefore, there is often no unique geographical region that can be defined as the host region. Using stochastic simulations (Boore, 2003) based on hypocentral distance, Scherbaum et al. (2006) applied a genetic algorithm to determine the combination of source, path, and site parameters that are most compatible with several GMMs, as defined by a minimum misfit. For GMMs developed through the use of the ergodic assumption (Anderson and Brune, 1999), the host region is really a pseudo-region in the sense that it does not relate to a single geographical region.

The CY14 model was developed using the ergodic assumption, pooling strong-motion data from around the world (particularly for moderate-to-large magnitude scenarios). However, differences in site response (Japan) and anelastic attenuation (Japan, Italy, and the Wenchuan sequence) were identified for certain regions, and bespoke coefficients for these “non-California” regions were presented. The implication is that the default CY14 model is developed with the intention that it be representative of ground motions in California. That said, Figure 1 shows the magnitude–distance distribution of the data used by CY14. For magnitudes above $\mathit{M}=6$ , a significant portion of the data is not from California and it is more appropriate to consider the CY14 model representing some pseudo-region with crustal characteristics similar to those of California, rather than California itself. That is, identifying host-region Fourier parameters strictly appropriate for California would not necessarily mean that these were compatible with CY14.

Figure 1.

Magnitude–distance distribution of data used to develop the Chiou and Youngs (2014) GMM.

We believe that it is, therefore, still more appropriate to define a suite of source, path, and site parameters that are consistent with the backbone GMM—which may, in some sense, correspond to a pseudo-region—rather than bring in additional uncertainty related to any misfit between this pseudo-region and average conditions in the assumed host region of California. The objective is to adjust the backbone model to match the target-region characteristics and for this purpose the actual provenance of the backbone GMM is not important, provided that the host-region parameters collectively provide simulations that are a very good approximation to the predictions from the empirical model. Another advantage of determining a suite of inverted parameters that are compatible with the backbone GMM is that it allows constraints to be applied that facilitate the application of the subsequent adjustments. For example, the geometric spreading can be modeled to mimic the spreading implicit in the backbone GMM, whereas published studies reporting source, path, and site parameters for the host region may have assumed different spreading functions, which would then need to be accounted for in applying an adjustment for anelastic attenuation. In addition, the inversion can focus on matching the backbone GMM in the magnitude and distance ranges of most relevance to the PSHA integrations for the target site. Moreover, for this specific application, we are fortunate to have access to a $V_S$ profile that has been defined on the basis of compatibility with the CY14 GMM (Al Atik and Abrahamson, 2021), which is adopted as a constraint for our inversions thus facilitating the definition of a truly model-consistent suite of parameters.

Inverting the CY14 GMM

Mathematical framework

Within the random vibration theory (RVT) framework, response spectral ordinates for a given natural period, $T_n$ (or natural frequency, $f_n=1/T_n$ ), and damping ratio ${\zeta }_n$ , are computed as:

$$Sa\left(f_n,{\zeta }_n\right)=\psi \left(f_n,{\zeta }_n\right){\omega }_n^2{u}_{rms}\left(f_n,{\zeta }_n\right)=\psi \left(f_n,{\zeta }_n\right)\sqrt{\frac{m_0\left(f_n,{\zeta }_n\right)}{D_{rms}\left(f_n,{\zeta }_n\right)}}$$

(1)

where $\mathrm{Sa}(f_n,{\zeta }_n)$ is the pseudo-spectral acceleration, ${\omega }_n=2\pi f_n$ is the natural angular frequency of the oscillator, $u_{\mathrm{rms}}$ is the root-mean-square (RMS) displacement response of the oscillator, $\psi$ is the peak factor (the expected value of the ratio between the maximum and RMS displacement response: $\psi \equiv u_{\max }/u_{\mathrm{rms}}$ ), $m_0$ is the zeroth spectral moment of the oscillator response, and $D_{\mathrm{rms}}$ is the root-mean-square duration. All of these components depend on both $f_n$ and ${\zeta }_n$ .

The zeroth moment is computed using the Fourier amplitude spectrum (FAS) of the oscillator’s displacement response $\left|U(f;f_n,{\zeta }_n)\right|$ , and the general expression for the kth moment:

$$\begin{array}{c}\multicolumn{1}{c}{m_k(f_n,{\zeta }_n)=2{\omega }_n^4\underset{0}{\overset{\infty }{\int }}{\left(2\pi f\right)}^k{\left|U(f;f_n,{\zeta }_n)\right|}^2\mathrm{df}}\\ \multicolumn{1}{c}{=2\underset{0}{\overset{\infty }{\int }}{\left(2\pi f\right)}^k{\left|H(f;f_n,{\zeta }_n)\right|}^2{\left|A\left(f\right)\right|}^2\mathrm{df}}\end{array}$$

(2)

with $k=0$ , $\left|H(f;f_n,{\zeta }_n)\right|$ being the modulus of the frequency response function of the oscillator, and $\left|A\left(f\right)\right|$ being the FAS of the ground acceleration. The modulus of the oscillator frequency response function is defined as:

$$\left|H(f,f_n,{\zeta }_n)\right|=\frac{1}{\sqrt{{[1-{(f/f_n)}^2]}^2+{(2{\zeta }_{n}f/f_n)}^2}}\equiv \frac{{\omega }_n^2\left|U(f;f_n,{\zeta }_n)\right|}{|A(f)|}$$

(3)

Note that this frequency response function differs slightly from that in the study by Boore (2003) as it embeds the influence of ${\omega }_n^2$ . The FAS of acceleration is specified in terms of a source spectrum, $E\left(f\right)$ , and path, $P\left(f\right)$ , and site, $S\left(f\right)$ , modifiers as:

$$\left|A\left(f\right)\right|={\left(2\pi f\right)}^{2}E\left(f\right)P\left(f\right)S\left(f\right)$$

(4)

The overview article of Boore (2003) provides detailed descriptions of these components, including a number of alternative representations that can be used for the source spectrum.

In this study, the primary objective is to obtain a set of FAS parameters that can be used for making comparisons against equivalent parameter sets inverted from empirical data in some target region. Such data are overwhelmingly dominated by recordings of rupture scenarios for which a simple point-source representation of the earthquake source is appropriate. Comparisons between host-region and target-region FAS parameter sets are greatly facilitated if the same parametric formulation is adopted. In addition, with the consideration of an appropriate distance metric, a simple point-source spectrum can also represent motions from larger events that will often be considered within hazard calculations (Boore, 2009). Therefore, while the CY14 GMM is applicable for very large events where double-corner frequency source spectra (Atkinson and Silva, 2000; Boore et al., 2014) can perform well, the source spectrum used in this study is a relatively simple single-corner frequency ${\omega }^2$ spectrum (Aki, 1967) defined by:

$$E\left(f\right)=\frac{\mathcal{C}{M}_0}{1+{(f/f_c)}^2}$$

(5)

The corner frequency, $f_c$ , in Equation 5 is defined in terms of the seismic moment, $M_0$ (Nm), and stress parameter, $\mathrm{\Delta }\sigma$ (MPa), following Brune (1970, 1971)

$$f_c=4.9058\times {10}^4{\beta }_s{\left(\frac{\mathrm{\Delta }\sigma }{M_0}\right)}^{1/3}$$

(6)

The term $\mathcal{C}$ is defined as:

$$\mathcal{C}=\frac{{\mathcal{R}}_{\theta \varphi }\mathrm{VF}}{4\pi {\rho }_s{\beta }_s^3{r}_0}\times {10}^{-13}$$

(7)

with ${\mathcal{R}}_{\theta \varphi }$ being an average radiation pattern (taken to be 0.55, as is common), $V=1/\sqrt{2}$ is used to partition energy among horizontal components, $F=2$ represents the free-surface effect, ${\rho }_s=2.75$ (mg/m³) and ${\beta }_s=3.5$ (km/s) are the density and S-wave velocity at the source (Al Atik and Abrahamson, 2021), and $r_0=1$ km is the reference distance to which the source amplitude is anchored. When the seismic moment is specified in units of Nm, and other terms adopt the units previously indicated, the 1×10⁻¹³ factor will define spectral amplitudes in Equation 4 with units of cm/s. It is important that the reference distance is specified to be consistent with the geometric spreading function, within $P\left(f\right)$ , and the source density and velocity should also be consistent with the impedance effects encapsulated within the site response $S\left(f\right)$ . That is, distance scaling and impedance effects must be modeled relative to the reference distance and source impedance.

The path scaling is comprised of geometric spreading, $g\left(r\right)$ , and anelastic attenuation, $q(f,r)$ , effects.

$$P\left(f\right)=g\left(r\right)q(f,r)=g\left(r\right)\exp [-\frac{\pi \mathrm{fr}}{Q_0{f}^{\eta }{c}_Q}]$$

(8)

The geometric spreading functions considered will be discussed subsequently, as will the choice of distance metric $r$ . In Equation 8, $c_Q$ is the velocity (km/s) used to determine the quality factor $Q_0$ . These two terms clearly trade-off directly with each other, so $c_Q\equiv {\beta }_s$ is typically assumed.

Finally, the site response $S\left(f\right)$ is comprised of both impedance and damping effects. Impedance effects are typically computed using the quarter wavelength approximation (Boore, 2003, 2013), which requires the specification of a crustal profile (velocity and density variation with depth). Damping effects are modeled using a “kappa” filter of the form $\exp (-\pi {\kappa }_{0}f)$ , where the ${\kappa }_0$ term should be consistent with the crustal profile used for the impedance effects and the path attenuation, $q(f,r)$ . The CY14 GMM has an elaborate site response formulation that is capable of providing estimates of spectral amplitudes for a broad range of $V_{S30}$ values. In principle, this range of $V_{S30}$ values corresponds to a range of crustal profiles. However, for this study, we focus upon a single value of $V_{S30}=760$ m/s, and therefore require the specification of a single profile, or a single $S\left(f\right)$ , that is consistent with the target $V_{S30}$ . Al Atik and Abrahamson (2021) have recently developed a method that can be used to determine a representative profile, and hence $S\left(f\right)$ , for a given GMM, and demonstrated the method using models including CY14. We, therefore, adopt their findings for $V_{S30}=760$ m/s directly and constrain $S\left(f\right)$ within the inversion. The impedance function and overall site response from Al Atik and Abrahamson (2021) is shown in Figure 2. The Al Atik and Abrahamson (2021) model is a combination of both impedance effects and damping, with ${\kappa }_0=0.039$ s, and these components should be used together.

Figure 2.

Impedance effects and overall site response. The model of Al Atik and Abrahamson (2021) is used to constrain $S\left(f\right)$ , while the Boore (2016) model is shown for comparative purposes. Curves with damping effects included correspond to ${\kappa }_0=0.039$ s in both cases.

As noted in the previous section, we ideally want model-specific rather than host region–specific parameters, and Figure 2 demonstrates the potential implications of this point. In the absence of the Al Atik and Abrahamson (2021) model, an obvious choice to constrain the impedance effects within $S\left(f\right)$ would be to adopt the generic velocity profile of Boore (2016) for the same target condition of $V_{S30}=760$ m/s. Although the Boore (2016) impedance function can reasonably be regarded as representative of the host region associated with CY14, Figure 2 shows some important differences over a broad range of frequencies. If the Boore (2016) impedance function were adopted, we would need to invert for a value of ${\kappa }_0$ , and would also need the other FAS parameters to accommodate the differences seen in Figure 2.

The remaining terms of Equation 1 yet to be defined are $D_{\mathrm{rms}}$ and $\psi$ . The RMS duration, $D_{\mathrm{rms}}$ , is computed from the excitation duration, $D_{\mathrm{ex}}$ , as:

$$D_{\mathrm{rms}}=\mathrm{\Gamma }(f_n,{\zeta }_n)D_{\mathrm{ex}}$$

(9)

where $\mathrm{\Gamma }$ is itself a function of the excitation duration, in addition to the natural frequency of the oscillator and the damping ratio. This term $\mathrm{\Gamma }$ is computed as the adjustment to the excitation duration that is required for the RVT-based estimates of peak oscillator response to match those obtained from time-domain simulations (Boore and Thompson, 2012, 2015) In this study, the $\mathrm{\Gamma }$ factors for active crustal regions from the study by Boore and Thompson (2015) are assumed appropriate and are adopted. The overall excitation duration is represented as:

$$D_{\mathrm{ex}}=D_{\mathrm{src}}+D_{\mathrm{path}}$$

(10)

where $D_{\mathrm{src}}$ is the source duration and $D_{\mathrm{path}}$ is the path duration. The study by Boore and Thompson (2014) presented a model for path duration that was exclusively a function of the rupture distance, $R_{\mathrm{RUP}}$ . This model was revisited by Boore and Thompson (2015), and was deemed appropriate for continued use with the caveat that $R_{\mathrm{RUP}}$ be replaced by the equivalent point-source distance, $R_{\mathrm{PS}}$ . The Boore and Thompson (2015) path duration model is used within this study, but the particular definition of $R_{\mathrm{PS}}$ can vary, as discussed subsequently. The source duration is represented as $D_{\mathrm{src}}=1/f_c$ and is, therefore, a function of the seismic moment and stress parameter. As a result, in this study, the only component of the duration model that is inverted for is $\mathrm{\Delta }\sigma$ —which controls the source duration within $D_{\mathrm{ex}}$ through its relationship with the corner frequency.

Finally, the peak factor $\psi$ is computed following the recommendations of Boore and Thompson (2015). Specifically, this entails the use of the Der Kiureghian (1980) modification to the peak factor distribution by Vanmarcke (1975).

All terms of Equation 1 have now been defined, and this set of equations collectively define the procedure for computing a response spectral ordinate via the RVT procedure. As can be seen, a number of components are constrained within the framework. Specifically, the path duration and $\mathrm{\Gamma }$ duration correction factors are fixed, and the site response $S\left(f\right)$ is fixed. As previously noted, the stress parameter $\mathrm{\Delta }\sigma$ is a free parameter and influences the source duration. A number of other parameters, related to geometric spreading, source depth scaling, and anelastic attenuation, influence $\left|A\left(f\right)\right|$ . As all spectral moments depend on $\left|A\left(f\right)\right|$ , both the zeroth moment and the peak factor in Equation 1 depend on these FAS parameters. We will collectively define all free FAS parameters by the vector $\mathit{\theta }$ .

In addition, while not explicitly noted in the equations presented previously, all terms of Equation 1 depend on variables that characterize the rupture scenario, such as magnitude, distance, and depth. These independent variables, along with $T_n$ and ${\zeta }_n$ , are collectively represented by ${\mathit{x}}_{\mathrm{ijkl}}$ , where the indexing is $i$ for depth, $j$ for distance, $k$ for magnitude, and $l$ for period. Therefore, once a particular rupture scenario and single-degree-of-freedom (SDOF) oscillator are defined, that is, ${\mathit{x}}_{\mathrm{ijkl}}$ is formed, the spectral acceleration associated with these properties is represented as ${\mathrm{Sa}}_{\mathrm{ijkl}}^{\mathrm{RVT}}({\mathit{x}}_{\mathrm{ijkl}};\mathit{\theta })$ . For the same set of independent variables, the CY14 model can be evaluated to compute ${\mathrm{Sa}}_{\mathrm{ijkl}}^{\mathrm{CY}}\left({\mathit{x}}_{\mathrm{ijkl}}\right)$ . The computation of ${\mathrm{Sa}}_{\mathrm{ijkl}}^{\mathrm{CY}}\left({\mathit{x}}_{\mathrm{ijkl}}\right)$ is not quite a direct evaluation of the CY14 GMM, the reason being that the constrained $S\left(f\right)$ corresponds to purely linear site response. Values of ${\mathrm{Sa}}_{\mathrm{ijkl}}^{\mathrm{CY}}\left({\mathit{x}}_{\mathrm{ijkl}}\right)$ are, therefore, computed by enforcing pure linear site response within the CY14 GMM. The host-to-target site adjustments required for any given application can then define linear and nonlinear site corrections relative to linear site response within the host model.

To identify the optimal FAS parameters, we solve for the $\mathit{\theta }$ that minimizes the function in Equation 11:

$$\ln \mathcal{L}(\mathit{\theta };\mathit{x})=\sum _{l=1}^{n_T}\sum _{k=1}^{n_M}\sum _{j=1}^{n_R}\sum _{i=1}^{n_Z}{w}_i{[\ln {\mathrm{Sa}}_{\mathrm{ijkl}}^{\mathrm{CY}}({\mathit{x}}_{\mathrm{ijkl}})-\ln {\mathrm{Sa}}_{\mathrm{ijkl}}^{\mathrm{RVT}}({\mathit{x}}_{\mathrm{ijkl}};\mathit{\theta })]}^2$$

(11)

The summation over the indices $i$ , $j$ , $k$ , and $l$ correspond to the consideration of all rupture scenarios ( $n_Z$ depths, $n_R$ distances, and $n_M$ magnitudes) and $n_T$ response spectral ordinates used within the inversion. The term $w_i$ represents weights of a discrete depth distribution, to be defined in detail later. All numerical analyses are implemented in Julia (Bezanson et al., 2017), with all RVT operations performed using the bespoke StochasticGroundMotionSimulation.jl package (Stafford, 2021b). Optimal parameter estimates are obtained using gradient-based optimization algorithms implemented in the NLopt library (Johnson, 2021b) via the NLopt.jl Julia interface (Johnson, 2021a). While multiple algorithms were tested, the sequential least squares quadratic programming algorithm, SLSQP, of Kraft (1988, 1994) was found to be relatively fast and stable for this particular problem. Gradients required for the optimization process were computed exactly using automatic differentiation (Molkenthin et al., 2014) via the ForwardDiff.jl package (Revels et al., 2016).

The weighted sum-of-squares expression in Equation 11 is essentially proportional to the log-likelihood function. As such, estimates of the standard errors and covariances among the FAS parameters can be computed using the Cramér–Rao bounds and the Fisher information matrix (Stafford, 2019). The Cramér–Rao bounds define the covariance matrix of the parameters $\mathit{\theta }$ as:

$${\mathbf{\sum }}_{\widehat{\mathit{\theta }}}\ge \mathit{I}{(\widehat{\mathit{\theta }})}^{-1}$$

(12)

where $\mathit{I}(\widehat{\mathit{\theta }})$ is the Fisher information matrix evaluated at the optimal parameter estimates $\hat{\mathit{\theta }}$ . When we have $p$ parameters, the information matrix will be a $p\times p$ matrix with individual elements defined by:

$${\mathit{I}}_{ij}(\widehat{\mathit{\theta }})=-\mathbb{E}\left[\frac{{\partial }^2\ln ℒ(\widehat{\mathit{\theta }};\mathit{x})}{\partial {\theta }_i\partial {\theta }_j}\right]$$

(13)

That is, the information matrix is the expectation of the Hessian matrix evaluated at the optimal parameter estimates. The computation of these mixed second-order partial derivatives for the mathematical framework presented throughout this section appears to be a formidable undertaking. However, like the gradients required for the optimization procedure, the Hessian matrix is readily evaluated using exact automatic differentiation via ForwardDiff.jl (Revels et al., 2016), as StochasticGroundMotionSimulation.jl was developed with this functionality in mind.

Parameter space

It is well known that significant parameter correlations can be encountered when performing Fourier spectral inversions (Boore, 2012; Boore et al., 1992). For this reason, inversions of empirical Fourier spectra are generally conducted in a sequence of steps in the hope of decoupling potential dependencies, as discussed in detail by Abercrombie (2021). When working with empirical data, the parameter trade-offs that are observed arise from multiple sources. Although each of the Fourier spectral parameters introduced in the previous section plays a distinct role in influencing either the amplitude or shape of spectra, the combination of the variability in Fourier spectral ordinates and the bandwidth limitations of recording instruments dictates that their unique role can rarely be isolated. Compounding this issue is the fact that the rupture scenarios represented within empirical databases impose an implicit weighting upon certain parameters over others. For example, many empirical datasets are dominated by recordings of small-magnitude events at large distances (Figure 1). Such scenarios produce recordings with a limited usable bandwidth at relatively high frequencies where the influence of $\mathrm{\Delta }\sigma$ , $Q$ , and ${\kappa }_0$ , all play a role.

When inverting a response spectral GMM, there are different challenges to address to limit parameter correlations. The two primary issues relate to how the parameter space is sampled, and how the influence of FAS parameters manifests for response spectral ordinates (Bora et al., 2016; Stafford et al., 2017). To ensure that constraint can be imposed upon each FAS parameter of interest, it is important to sample the parameter space of independent variables to include rupture scenarios for which the loss function of Equation 11 exhibits the greatest sensitivity to each parameter.

The value of the loss function in Equation 11 depends on the particular rupture scenarios and periods that are considered. For any given set of $n$ RVT calculations, a particular set of optimal parameters $\hat{\mathit{\theta }}$ will be obtained. Therefore, an objective of designing the sampling of the parameter space is to specify the ${\mathit{x}}_{\mathrm{ijkl}}$ scenarios in such a manner that $\hat{\mathit{\theta }}$ is relatively stable. That is, adding or removing scenarios from ${\mathit{x}}_{\mathrm{ijkl}}$ would not lead to significant changes to $\hat{\mathit{\theta }}$ . Ensuring a broad coverage of the magnitude, distance, depth, and period space helps to provide this stability. However, each scenario that is considered also adds to the computational demands of the inversion. We, therefore, look to strike a balance between the coverage of the parameter space and the robustness of the $\hat{\mathit{\theta }}$ estimates.

To explore and aid the design of the parameter space, formal sensitivity analyses were conducted. Starting with some initial estimates of parameter values, ${\mathit{\theta }}_0\equiv \{\mathrm{\Delta }\sigma =100,{\gamma }_1=1.15,Q_0=200,\eta =0.65\}$ , which map into $y\left({\mathit{\theta }}_0\right)\equiv \ln {\mathrm{Sa}}_{\mathrm{ijkl}}^{\mathrm{RVT}}({\mathit{x}}_{\mathrm{ijkl}};{\mathit{\theta }}_0)$ , the relative sensitivities can be computed as (Molkenthin et al., 2014):

$${\mathrm{\Phi }}_i=\frac{{\mathit{\theta }}_{0,i}}{y\left({\mathit{\theta }}_0\right)}{\frac{\partial y\left(\mathit{\theta }\right)}{\partial {\theta }_i}|}_{{\mathit{\theta }}_0}\approx \left(\frac{\mathrm{\Delta }y\left(\mathit{\theta }\right)}{y\left({\mathit{\theta }}_0\right)}\right)\left(\frac{{\mathit{\theta }}_{0,i}}{\mathrm{\Delta }{\theta }_i}\right)$$

(14)

This expression determines how a percentage change in ${\theta }_i$ maps into a percentage change in $y\left(\mathit{\theta }\right)\equiv \ln S{a}^{\mathrm{RVT}}(\mathit{x};\mathit{\theta })$ . As the loss function of Equation 11 depends directly on $\ln S{a}^{\mathrm{RVT}}(\mathit{x};\mathit{\theta })$ , it is important to sample the parameter space to ensure that scenarios that lead to strong values (positive or negative) of ${\mathrm{\Phi }}_i$ are obtained. Note that as $S{a}^{\mathrm{RVT}}(\mathit{x};\mathit{\theta })<1$ over the parameter space considered, and all elements of $\mathit{\theta }$ are positive, negative values of ${\mathrm{\Phi }}_i$ correspond to situations where $\partial \ln \mathrm{Sa}/\partial {\theta }_i>0$ , and vice versa.

Figure 3 shows the relative sensitivities to four key FAS parameters for short ( $T=0.05$ s), intermediate ( $T=0.5$ s) and long ( $T=5.0$ s) periods. The base FAS model uses the source spectrum of Equation 5, that depends on $\mathrm{\Delta }\sigma$ , the anelastic attenuation formulation within Equation 8, depending on both $Q_0$ and $\eta$ , and a geometric spreading function that smoothly transitions from a near-field geometric spreading rate of ${\gamma }_1$ to a fixed far-field spreading rate of ${\gamma }_f=0.5$ . This spreading function is shown in Equation 15, where $r_t=50$ km and $r_0=1$ km:

$$\ln g\left(R_{\mathrm{PS}}\right)=-{\gamma }_1\ln \left(R_{\mathrm{PS}}\right)+\frac{({\gamma }_1-{\gamma }_f)}{2}\ln \left(\frac{R_{\mathrm{PS}}^2+r_t^2}{r_0^2+r_t^2}\right)$$

(15)

Figure 3.

Relative sensitivities, ${\Phi }_i$ (Equation 14) representative of periods in the short, intermediate, and long range. Each column corresponds to the period shown in the column title, while each row relates to the parameter specified in the row annotation. Common color ranges are adopted for all panels in each row.

In Equation 15, the equivalent point-source distance $R_{\mathrm{PS}}$ , to be defined later, is computed using the model of Boore and Thompson (2015).

Consideration of figures like those shown in Figure 3 ensures that the ${\mathit{x}}_{\mathrm{ijkl}}$ combinations ultimately selected should have a good chance of isolating the effects of individual parameters. For example, the left and central columns of Figure 3 both show that the sensitivity to $Q_0$ and $\eta$ follow very similar patterns in magnitude–distance space (although the shorter periods have a stronger relative sensitivity to $\eta$ ). This similarity in sensitivity leads to strong correlations among $Q_0$ and $\eta$ if only these short-to-intermediate periods are considered—which is often necessarily the case for empirical data. However, the right column of Figure 3 shows that the sensitivity to these parameters diverges for large-magnitude, long-distance, scenarios. The key to decoupling the correlation between $Q_0$ and $\eta$ is to ensure that these “upper right” regions of the magnitude–distance space are sampled for long periods.

Similar considerations can be used to isolate the role of other parameters. For example, the left column of Figure 3 shows that ${\gamma }_1$ plays a significant role for large-magnitude, short-distance scenarios—as it interacts strongly with the near-source saturation constraints. Furthermore, in the same column, Figure 3 also shows a strong local sensitivity to $\mathrm{\Delta }\sigma$ for intermediate magnitudes at short distances.

In addition to the above considerations, it is also important to ensure that the inverted parameter set leads to RVT predictions that mimic the CY14 GMM over a broad range of scenarios for which this GMM is applicable. If good performance over the full range of applicability is not obtained, then complications can arise when developing host-to-target adjustments as a result of extrapolation issues.

After considering the above constraints, a very comprehensive set of 7,140 magnitude–distance depth scenarios was defined. For each of these scenarios, a set of 20 periods were considered. As a result, for every evaluation of Equation 11, 142,800 RVT calculations were performed (as well as the simultaneous computation of partial derivatives with respect to each considered FAS parameter). Naturally, within the optimization process, Equation 11 is evaluated a large number of times as we move from initial parameter estimates toward the final converged solution. Formally, the optimization is terminated when the relative tolerance in estimates of $\hat{\mathit{\theta }}$ drop below 1×10⁻⁵. The procedure is computationally intensive, but the package implementing the RVT procedure (Stafford, 2021b) takes a few tens of microseconds on a standard laptop for an individual RVT calculation, and the overall loss functions can be evaluated in several seconds (depending on the model). Performing the full optimization can take multiple hours, but this depends on the model, convergence tolerances, numbers of scenarios, and parameters.

The full specification of ${\mathit{x}}_{\mathrm{ijkl}}$ can be obtained from:

• $n_T=20$ periods, approximately logarithmically spaced, with specific values of 0.01, 0.02, 0.03, 0.04, 0.05, 0.075, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5, 0.75, 1.0, 1.5, 2.0, 3.0, 5.0, 7.5, and 10.0 s.

• $n_M=28$ moment magnitude, $\mathit{M}$ , values, linearly spaced with increments of $\mathrm{\Delta }\mathit{M}=0.2$ from 3.0 to 8.4, inclusive. This covers the full range of applicability of the CY14 GMM, with an additional extension to smaller magnitudes given their relevance for deriving model adjustments.

• $n_R=35$ distance values (defined as Joyner–Boore, $R_{\mathrm{JB}}$ , distances), with values of 0, 1, 2, 3, 5, 7.5, 10, 12.5, 15, 17.5, 20, 25, 30, 35, 40, 45, 50, 60, 70, 80, 90, 100, 120, 140, 150, 160, 170, 180, 190, 200, 220, 240, 260, 280, and 300 km. Tight spacing is adopted at short distances to help constrain the near-source saturation. This spacing coarsens with increasing distance, but is linear within certain distance ranges to help constrain the anelastic attenuation parameters. That is, we balance the implicit weight assigned to geometric spreading and anelastic attenuation parameters. As discussed in the following section, all considered rupture scenarios correspond to vertical strike–slip events, so the rupture distance required by CY14 is computed as $R_{\mathrm{RUP}}\left(\mathit{M}\right)=\sqrt{R_{\mathrm{JB}}^2+Z_{\mathrm{TOR}}^2\left(\mathit{M}\right)}$ , with $Z_{\mathrm{TOR}}\left(\mathit{M}\right)$ being the magnitude-dependent depth to the top of rupture.

• $n_Z=9$ depth values, defined as magnitude-dependent fractiles from the depth distribution shown in Figure 4. The depth distribution is derived from the model for the expected depth from CY14, and the empirical data underpinning that model. Specifically, the expected depth model (derived in $\sqrt{Z_{\mathrm{TOR}}}-\mathit{M}$ space) is adopted directly, and a model for the magnitude-dependent standard deviation of $\sqrt{Z_{\mathrm{TOR}}}$ is derived from the empirical data. These parameters are assumed to define a normal distribution of $\sqrt{Z_{\mathrm{TOR}}}$ that is truncated to apply over the range $\sqrt{Z_{\mathrm{TOR}}}\in [0,20]$ . For each considered magnitude, this truncated normal distribution is then discretised into nine pairs of nodes, $z_{\mathrm{tor},i}$ , and weights, $w_i$ , using the approach of Miller and Rice (1983). These discrete depths are shown in Figure 4, with line thicknesses and color intensities representing the weights, $w_i$ . Note that these discrete depth weights correspond to the weights in Equation 11 and that $\sum _i{w}_i=1$ .

Figure 4.

Discrete depth distribution used in the definition of the parameter space. The empirical data used for the model development are shown using circular markers. The dark solid line shows the expected depth model of CY14, while light shaded lines show nodes of the discrete depth distribution. The square markers indicate the weighted mean of these discrete depths, confirming their consistency with CY14.

Note that the full factorial combination of magnitude, distance, and depth values leads to 8,820 rupture scenarios. The lower number (7,140) of scenarios actually used arises because the depth distribution collapses to a singular value for the largest events, that is, the depth to the top of rupture for very large events is assumed to be zero (Figure 4).

Host CY14 GMM for inversions

Because the source within the FAS model is treated as an equivalent point source, it is not possible to explicitly model effects associated with style of faulting and fault-finiteness (such as hanging-wall effects and directivity). For this reason, the rupture scenarios that are considered all correspond to vertical strike–slip ruptures, and directivity effects are ignored. This modeling decision effectively simplifies the CY14 model by rendering many of the functional terms of the model redundant. For vertical strike–slip ruptures, ignoring directivity effects (setting $\mathrm{\Delta }\mathrm{DPP}=0$ ), constraining site response to be purely linear, fixing $V_{S30}=760$ m/s and using the expected depth to the 1.0 km/s velocity horizon (setting $\mathrm{\Delta }{Z}_{1.0}=0$ ), the functional form of the CY14 GMM reduces to:

$$\begin{array}{c}\multicolumn{1}{c}{\ln \mathrm{Sa}=c_1\prime +\underset{\mathrm{depth}\mathrm{to}\mathrm{top}-\mathrm{of}-\mathrm{rupture}\mathrm{scaling}}{\underbrace{(c_7+\frac{c_{7b}}{\cosh [2max(\mathit{M}-4.5,0)]})\mathrm{\Delta }{Z}_{\mathrm{TOR}}}}}\\ \multicolumn{1}{c}{+\underset{\mathrm{magnitude}\mathrm{scaling}}{\underbrace{c_2(\mathit{M}-6)+\frac{c_2-c_3}{c_n}\ln [1+e^{c_n(c_M-\mathit{M})}]}}}\\ \multicolumn{1}{c}{+\underset{\mathrm{near}-\mathrm{source}\mathrm{saturation}\mathrm{and}\mathrm{geometric}\mathrm{spreading}}{\underbrace{c_4\ln (R_{\mathrm{RUP}}+c_5\cosh [c_{6}max(\mathit{M}-c_{\mathrm{HM}},0)])+(c_{4a}-c_4)\ln \left(\sqrt{R_{\mathrm{RUP}}^2+c_{\mathrm{RB}}^2}\right)}}}\\ \multicolumn{1}{c}{+\underset{\mathrm{anelastic}\mathrm{attenuation}\mathrm{scaling}}{\underbrace{(c_{{\gamma }_1}+\frac{c_{{\gamma }_2}}{\cosh [max(\mathit{M}-c_{{\gamma }_3},0)]})R_{\mathrm{RUP}}}}}\end{array}$$

(16)

where the $c_1\prime$ coefficient encapsulates the ${\varphi }_1\ln (760/1130)$ associated with linear site response. All other coefficients correspond directly to those of the CY14 GMM.

In Equation 16, the distinct functional terms that relate to $Z_{\mathrm{TOR}}$ scaling, magnitude scaling, geometric spreading (including near-source saturation), and anelastic attenuation are indicated. While we should not expect a one-to-one mapping between the functional terms of the FAS and those of Equation 16 (Bora et al., 2016), it does help to consider which terms depend on magnitude, and where breaks in scaling arise within the CY14 GMM. For example, the depth to top-of-rupture scaling is clearly magnitude dependent only for $\mathit{M}>4.5$ , so it is worth keeping that in mind when developing functional expressions within the FAS inversions. As another example, the anelastic attenuation term is also magnitude dependent, which suggests that one of the parameters $Q_0$ or $\eta$ may need to include magnitude dependence to replicate the scaling.

We can also see that the geometric spreading functional form shares a lot in common with the function shown previously in Equation 15. The equivalent point-source distance, $R_{\mathrm{PS}}$ , considered by Boore and Thompson (2015) can be generally represented as:

$$R_{\mathrm{PS}}={(R_{\mathrm{RUP}}^n+h{\left(\mathit{M}\right)}^n)}^{1/n}$$

(17)

where Boore and Thompson (2015) use $n=2$ , and $h\left(\mathit{M}\right)$ is their finite-fault factor model. In the CY14 GMM we have:

$$h\left(\mathit{M}\right)\equiv c_5\cosh [c_{6}max(\mathit{M}-c_{\mathrm{HM}},0)]$$

(18)

with $c_5$ , $c_6$ , and $c_{\mathrm{HM}}$ being period dependent, and $n=1$ in Equation 17. However, the second term of the geometric spreading function that controls the geometric spreading at large distances (where the interaction with anelastic attenuation is important), scales directly in terms of $R_{\mathrm{RUP}}$ rather than $R_{\mathrm{PS}}$ . These considerations of individual terms within the overall functional form of the CY14 GMM inform the parameterization of the FAS model used within the inversions.

Inversion results

A challenge often faced when developing model adjustments is the reconciliation of “apparent parametric” and “real physical” differences that arise from considering different magnitude ranges in the host and target regions. In particular, target region FAS parameter sets are often inferred from recordings of small magnitude events, while potential host regions are generally considered because they have empirical constraint for larger rupture scenarios of greater relevance to hazard calculations. It is, therefore, important to ensure that the inversion of the CY14 GMM is robust for small magnitude events, and that any magnitude dependence in the FAS parameters is adequately modeled. When this is done successfully, differences in host and target region parameter sets obtained over the small magnitude range can be appropriately scaled to corresponding differences at larger magnitudes (assuming the physical FAS model is appropriate).

As a result, the first step of the inversion procedure is to investigate any magnitude dependence of the FAS parameters by performing inversions for specific magnitude values. This first step is discussed in the following section, and informs the parameterization of the full inversions presented subsequently.

Initial magnitude-dependent scaling

For a fixed value of magnitude, it can be appreciated that many elements of the functional complexity within Equation 16 disappear. To mimic the functional form of Equation 16, a single-corner source spectrum was adopted, with the stress parameter made a function of the relative depth to the top-of-rupture, $\mathrm{\Delta }{Z}_{\mathrm{TOR}}$ . The stress parameter scaling is therefore:

$$\ln \mathrm{\Delta }\sigma =\ln \mathrm{\Delta }{\sigma }_{\mathit{M}}+{\delta }_{\ln \mathrm{\Delta }\sigma ,z}\mathrm{\Delta }{Z}_{\mathrm{TOR}}$$

(19)

where $\ln \mathrm{\Delta }{\sigma }_{\mathit{M}}$ is the expected value of $\ln \mathrm{\Delta }\sigma$ for magnitude $\mathit{M}$ , and ${\delta }_{\ln \mathrm{\Delta }\sigma ,z}$ determines how strongly changes to the depth of rupture influence the logarithmic stress parameter. The geometric spreading function is defined as:

$$\ln g(R_{\mathrm{PS}},R_{\mathrm{RUP}})=-{\gamma }_1\ln \left(R_{\mathrm{PS}}\right)+\frac{({\gamma }_1-{\gamma }_f)}{2}\ln \left(\frac{R_{\mathrm{RUP}}^2+r_t^2}{r_0^2+r_t^2}\right)$$

(20)

with $r_0=1$ km, $r_t=50$ km, ${\gamma }_f=0.5$ , and the equivalent point-source distance defined as $R_{\mathrm{PS}}=R_{\mathrm{RUP}}+h_{\mathit{M}}$ . In this formulation, $h_{\mathit{M}}$ is the best estimate of the saturation distance for magnitude $\mathit{M}$ . The anelastic attenuation is parameterized according to Equation 8, with $R_{\mathrm{RUP}}$ used as the distance metric. Site effects associated with impedance and damping are constrained to equal the Al Atik and Abrahamson (2021) model. Therefore, with this model formulation, there are six free FAS parameters for each magnitude level considered: $\mathrm{\Delta }{\sigma }_{\mathit{M}}$ , ${\delta }_{\ln \mathrm{\Delta }\sigma ,z}$ , ${\gamma }_1$ , $h_{\mathit{M}}$ , $Q_0$ , and $\eta$ . However, for the largest magnitudes, only five of these parameters are relevant as $Z_{\mathrm{TOR}}=0$ km for all scenarios and ${\delta }_{\ln \mathrm{\Delta }\sigma ,z}=0$ as a result.

Inversions are performed for the periods, distances, and depths previously presented in the Parameter Space section, and for individual magnitudes from 3.0 to 8.4 in 0.2 unit increments. The results of these initial inversions are presented in Figure 5. In all cases, both the parameter estimates and the associated standard errors are shown. Figure 5a shows that stress parameter estimates initially grow exponentially with magnitude, up to $\mathit{M}=5$ , before reaching a plateau over the approximate range $\mathit{M}\in (5,7)$ . Thereafter, the stress parameter increases at an exponential rate, and reaches extremely high values that would not ordinarily be regarded as physically meaningful. However, Figure 5c also shows that the geometric spreading rate is approximately constant for small to moderate magnitudes before increasing significantly for the largest magnitudes. For the largest magnitudes, these parameters have a strong positive correlation, with high stress parameter values trading off with high geometric spreading rates. This correlation can be observed in Figure 6, along with all other parameter correlations. In cases like this, the apparent magnitude dependence needs to be regarded with suspicion. However, previous studies are consistent with the result that stress parameter values increase over the small magnitude range before reaching an approximate plateau for larger magnitudes (Rietbrock et al., 2013; Yenier and Atkinson, 2015b). This prior intuition is used to inform the parametric form used for the full inversion, and to help decouple the correlation between stress parameter and geometric spreading rates, to some extent.

Figure 5.

Variation of FAS parameters with magnitude, from magnitude-specific inversions. (a) Logarithmic stress parameter estimates are shown, along with annotated $\mathrm{\Delta }\sigma$ (MPa) reference values. (b) The stress parameter increment associated with $\mathrm{\Delta }{Z}_{\mathrm{TOR}}$ scaling is shown, with “CY sech” showing the functional scaling of $\mathrm{\Delta }{Z}_{\mathrm{TOR}}$ within the CY14 GMM. (c) The geometric spreading rate is shown. (d) The near-source saturation distances obtained from the inversions are shown, along with other models. The CY14 ranges indicate the variation of $h(\mathit{M})$ with period, “BT15” is the Boore and Thompson (2015) model, and “YA15” is the Yenier and Atkinson (2015a) model. (e) and (f) The anelastic parameters $Q_0$ and $\eta$ , respectively, are shown.

Figure 6.

Parameter correlations from magnitude-specific inversions. (a) The variation of parameter correlations with magnitude are shown. Vertical gray dashed lines correspond to the magnitudes for which the correlation matrices are plotted in (b) and (c).

Figure 5b shows the influence of $\mathrm{\Delta }{Z}_{\mathrm{TOR}}$ on $\mathrm{\Delta }\sigma$ via the parameter ${\delta }_{\ln \mathrm{\Delta }\sigma ,z}$ from Equation 19. For small magnitudes, estimates of ${\delta }_{\ln \mathrm{\Delta }\sigma ,z}$ are well constrained, as the discrete depth nodes considered within the inversion span a broad range of depths (Figure 4). For moderate to large magnitudes, the range of $\mathrm{\Delta }{Z}_{\mathrm{TOR}}$ values reduces significantly and translates into very uncertain estimates of ${\delta }_{\ln \mathrm{\Delta }\sigma ,z}$ . However, Figure 5b shows a superimposed function, “CY sech,” that has the same functional form as the $\mathrm{\Delta }{Z}_{\mathrm{TOR}}$ scaling in Equation 16. While there are clear deviations for the smallest magnitudes, the best estimates of ${\delta }_{\ln \mathrm{\Delta }\sigma ,z}$ match the scaling embedded within the CY14 GMM very well. This finding is used to inform the parameterization of ${\delta }_{\ln \mathrm{\Delta }\sigma ,z}$ in the full inversion.

Figure 5d compares multiple models for the near-source saturation distance, $h\left(\mathit{M}\right)$ . In particular, two published models of Boore and Thompson (2015) and Yenier and Atkinson (2015a) are shown, along with the range of $h\left(\mathit{M}\right)$ predictions associated with Equation 18 over the full period range of the model. As can be appreciated, the estimates of $h_{\mathit{M}}$ from the inversion are very close to those of Boore and Thompson (2015) for very small magnitudes, and all models become remarkably consistent for the largest magnitudes—albeit with very different slopes. However, it is important to note that the effective $h\left(\mathit{M}\right)$ model of CY14 gives significantly longer saturation distances for the smallest events. This is important because a key reason why $\mathrm{\Delta }\sigma$ values decrease for the smallest magnitude events is that the $h\left(\mathit{M}\right)$ values over this range are larger than would be expected from physical arguments (Atkinson et al., 2016; Boore, 2012; Boore and Thompson, 2015).

Finally, Figure 5e suggests that $Q_0$ is largely insensitive to magnitude, while $\eta$ indicates a relatively strong magnitude dependence. Note from Equation 8 that values of $\eta \to 1$ correspond to a situation where $\ln q(f,r)\approx -\frac{\pi }{Q_0{c}_Q}r$ , that is, the anelastic attenuation becomes frequency independent and filters all Fourier ordinates by a similar amount. Large values of $\eta$ are, therefore, similar to using a different functional form for geometric spreading.

The correlation values in Figure 6 are encouraging in that we do not observe parameter combinations where the correlation stays very strong (positive or negative) over the full magnitude range. The most extreme correlations correspond to the combinations of $\{\ln \mathrm{\Delta }{\sigma }_{\mathit{M}},{\gamma }_1\}$ and $\{Q_0,\eta \}$ . The first of these was previously discussed, and can be addressed, in part, through the imposition of a functional form for $\mathrm{\Delta }\sigma (\mathit{M})$ in the full inversion. However, the correlation between $Q_0$ and $\eta$ is harder to deal with, as it arises largely from the fact that high-frequency response spectral ordinates do not correlate strongly with high-frequency Fourier spectral ordinates, where $\eta$ is most influential (Bora et al., 2016).

Optimal parameterization

Following the magnitude-specific inversions discussed in the previous section, a number of alternative parameterizations for the full inversion were considered. These alternatives were judged on the basis of their “scores,” as defined by Equation 11, and upon their general ability to mimic the scaling within the CY14 GMM in a parsimonious way. The latter characteristic was judged from visual inspection of diagnostic plots, but is also clearly related to the formal metric provided through Equation 11.

Note that while physical reasoning influences the selected parameterization in part, the ultimate goal is to derive a set of FAS parameters that reproduce the scaling embedded within the CY14 GMM. It must also be recalled that the specification of the FAS model is only one part of the overall inversion, as the site response, duration models, and peak factor model are all constrained. The specific scaling embedded within these RVT components also influences how the FAS model should behave. For example, Equation 1 shows that spectral amplitudes scale with both $m_0$ , which depends on the FAS, and $D_{\mathrm{rvt}}$ , which includes a fixed distance dependence within the path-scaling component. Therefore, to match the complex geometric spreading within the CY14 GMM, we need to find the FAS parameterization that combines with the fixed distance dependence within the duration model in the most appropriate manner. This point is easy to overlook as it is common to make direct comparisons between the scaling of response and Fourier spectral amplitudes.

For the full inversion, the stress parameter was defined to include explicit magnitude dependence for small magnitudes, and to also embed the effects of the depth to top of rupture. The magnitude dependence largely offsets the near-source saturation distance of the CY14 GMM, as previously discussed. The scaling with $\mathrm{\Delta }{Z}_{\mathrm{TOR}}$ should not necessarily relate to the stress parameter that primarily controls high-frequency Fourier amplitudes. For example, factors such as the source rigidity, ${\mu }_s={\rho }_s{\beta }_s^2$ , appearing in Equation 7, will reduce at shallow depths and should effect all FAS ordinates across the entire frequency band. However, the $\mathrm{\Delta }{Z}_{\mathrm{TOR}}$ scaling within CY14 is strongest at short periods, and weakens for longer periods. As such, embedding the $\mathrm{\Delta }{Z}_{\mathrm{TOR}}$ scaling within the stress parameter allows the source amplitudes to still be controlled by a single effective parameter and mimics the scaling of CY14 well. The functional form for the stress parameter is, therefore, defined as:

$$\begin{array}{c}\multicolumn{1}{c}{\ln \mathrm{\Delta }\sigma (\mathit{M},\mathrm{\Delta }{Z}_{\mathrm{TOR}})=s_{\alpha }+s_{\beta }min(\mathit{M}-5,0)}\\ \multicolumn{1}{c}{+(s_{\gamma }+s_{\delta }\mathrm{sech}[2max(\mathit{M}-4.5,0)])\mathrm{\Delta }{Z}_{\mathrm{TOR}}}\end{array}$$

(21)

where $s_{\alpha }$ , $s_{\beta }$ , $s_{\gamma }$ , and $s_{\delta }$ are all parameters solved for during the inversion. Note that the functional form for the $\mathrm{\Delta }{Z}_{\mathrm{TOR}}$ scaling is the same as that employed by CY14 as $1/\cosh \left(x\right)\equiv \mathrm{sech}\left(x\right)$ . While Figure 5b suggested that scaling with $\mathrm{\Delta }{Z}_{\mathrm{TOR}}$ may deviate from this form for very small magnitudes, we retain the CY14 parameterization.

Geometric spreading is defined in Equation 20, with the equivalent point-source distance defined as:

$$R_{\mathrm{PS}}=R_{\mathrm{RUP}}+h\left(\mathit{M}\right)$$

(22)

and the near-source saturation distance coming from:

$$\ln h\left(\mathit{M}\right)=h_{\alpha }+h_{\beta }\mathit{M}+\frac{h_{\beta }-h_{\gamma }}{h_{\delta }}\ln (1+e^{-h_{\delta }(\mathit{M}-h_{\epsilon })})$$

(23)

where $h_{\alpha }$ , $h_{\beta }$ , $h_{\delta }$ , and $h_{\epsilon }$ are free parameters to be solved within the inversion. The parameter $h_{\gamma }$ is constrained to equal $h_{\gamma }=\ln \left(10\right)/2\approx 1.15$ so that $h\left(\mathit{M}\right)$ scales in proportion to the expected source dimensions of the rupture for small events. Although the parameterization of Equation 23 is certainly more elaborate than the model of Yenier and Atkinson (2015a), shown in Figure 5d, it requires the specification of fewer parameters than the model of Boore and Thompson (2015) while aiming to achieve similar goals. That is, Equation 23 transitions from essentially exponential scaling with a rate of $h_{\gamma }$ at small magnitudes, to exponential scaling at a reduced rate of $h_{\beta }$ for large magnitudes. The parameter $h_{\delta }$ controls how sharply this transition occurs, while $h_{\epsilon }$ controls where the transition is centered. The functional form is the same as that used in the magnitude scaling term within Equation 16, with the significance of the various parameters clearly illustrated in the study by Chiou et al. (2010).

Within the inversion, a nonlinear constraint is imposed upon the combination of $h_{\beta }$ and ${\gamma }_1$ to avoid possible over-saturation of response spectral amplitudes for very large magnitudes. Initial inversions performed without imposing this constraint led to over-saturation—which does not exist within the CY14 GMM. Details of this nonlinear constraint are provided in the study by Stafford (2021a).

Finally, the quality exponent within the anelastic attenuation filter is modeled as magnitude dependent, based on the results presented in Figure 5f. The distance metric used within the anelastic attenuation filter is $R_{\mathrm{RUP}}$ . The quality factor, $Q_0$ , is a free parameter in the inversion, but $\eta$ is modeled as:

$$\eta \left(\mathit{M}\right)={\eta }_{\alpha }+{\eta }_{\beta }\tanh (\mathit{M}-{\eta }_{\gamma })$$

(24)

with ${\eta }_{\alpha }$ , ${\eta }_{\beta }$ , and ${\eta }_{\gamma }$ free parameters within the inversion. At very small magnitudes, the function tends to $\eta \left(\mathit{M}\right)\to {\eta }_{\alpha }-{\eta }_{\beta }$ ; while at very large magnitudes, we have $\eta \left(\mathit{M}\right)\to {\eta }_{\alpha }+{\eta }_{\beta }$ . There is no physical basis to this functional form, it is simply chosen to mimic the scaling observed in Figure 5f. However, aside from providing an improved fit over a magnitude-independent $\eta$ , that is, Equation 11 results in smaller values for $\eta \left(\mathit{M}\right)$ , there are additional reasons to include this magnitude dependence. Specifically, inversions of empirical data dominated by recordings of small magnitude events in California typically suggest relatively small values of $\eta$ (Lin and Jordan, 2018; Raoof et al., 1999), broadly consistent with the results shown in Figure 5d. If a magnitude-independent $\eta$ were adopted, this parameter would reflect the average over the $\eta (\mathit{M})$ model. Comparisons made with recordings of small magnitude events could then suggest apparent parametric differences in $\eta$ in the host and target regions, potentially leading to unnecessary model adjustments. Therefore, while magnitude-dependent $\eta$ is atypical for FAS parameter sets, there are good reasons to include it within this application.

In total, the optimal parameterization requires solving for 13 free parameters. The resulting parameter estimates and associated standard errors are presented in Table 1.

Table 1.

Estimates and standard errors for the optimal FAS parameters

Parameter	Estimate	Standard error	Parameter	Estimate	Standard error
$s_{\alpha }$	2.296	0.031	$h_{\gamma }$	1.1513	0.0
$s_{\beta }$	0.4624	0.0311	$h_{\delta }$	5.0948	0.725
$s_{\gamma }$	0.0453	0.0136	$h_{\epsilon }$	7.2725	0.0566
$s_{\delta }$	0.109	0.0166	$Q_0$	205.4	5.53
${\gamma }_1$	1.1611	0.00601	${\eta }_{\alpha }$	0.6884	0.0131
$h_{\alpha }$	−0.8712	0.373	${\eta }_{\beta }$	0.1354	0.00654
$h_{\beta }$	0.4451	0.0474	${\eta }_{\gamma }$	5.1278	0.0794

Note that $h_{\gamma }$ is fixed in the inversion, and that the $s_{\alpha }$ parameter leads to stress parameter values in units of MPa.

Convenience parameterization

While the optimal parameterization of the previous section is simple to implement in StochasticGround vMotionSimulation.jl, existing workflows based on alternative software may not be able to accommodate this parameterization. In particular, the mixing of distance metrics in Equation 20 and the magnitude dependence of $\eta$ in Equation 24 may be problematic. For this reason, a “convenience parameterization” was also considered.

This alternative parameterization closely follows the geometric spreading formulation by Boore (2003), and makes use of $R_{\mathrm{PS}}$ throughout both the geometric spreading and anelastic attenuation functions (in contrast to the optimal parameterization that used $R_{\mathrm{RUP}}$ as the distance metric for anelastic attenuation). However, $R_{\mathrm{PS}}$ remains calculated using Equation 17 with $n=1$ , that is, the form in Equation 22. In addition, only a magnitude-independent $\eta$ was considered.

The form of the geometric spreading adopted is trilinear in $\ln A(f)-\ln R_{\mathrm{PS}}$ space:

$$g\left(R_{\mathrm{PS}}\right)=\{\begin{array}{c}1\mathrm{for}{R}_{\mathrm{PS}}\le r_0\hfill \\ \multicolumn{1}{c}{{\left(\frac{r_0}{R_{\mathrm{PS}}}\right)}^{{\gamma }_1}\mathrm{for}{r}_0<R_{\mathrm{PS}}\le r_1}\\ \multicolumn{1}{c}{{\left(\frac{r_0}{r_1}\right)}^{{\gamma }_1}{\left(\frac{r_1}{R_{\mathrm{PS}}}\right)}^{{\gamma }_2}\mathrm{for}{r}_1<R_{\mathrm{PS}}\le r_2}\\ \multicolumn{1}{c}{{\left(\frac{r_0}{r_1}\right)}^{{\gamma }_1}{\left(\frac{r_1}{r_2}\right)}^{{\gamma }_2}{\left(\frac{r_2}{R_{\mathrm{PS}}}\right)}^{{\gamma }_3}\mathrm{for}{R}_{\mathrm{PS}}>r_2}\end{array}$$

(25)

with $r_1$ and $r_2$ being the hinge distances where the geometric spreading rates ${\gamma }_1$ , ${\gamma }_2$ , and ${\gamma }_3$ change. To identify appropriate hinge distances, an inversion was performed that sought to identify the parameters of Equation 25 that led to spectral ordinates as close as possible to those obtained from Equation 20—with all other FAS parameters fixed. The fixed FAS parameters corresponded to a preliminary model that made use of magnitude-independent $\eta$ , but otherwise had the same parameterization as the optimal model. For the optimal model, the overall geometric spreading is defined by the free parameter ${\gamma }_1$ , and the fixed parameters $r_t=50$ km and ${\gamma }_f=0.5.$ For the alternative geometric spreading function of Equation 25, the far-field spreading rate ${\gamma }_3$ was fixed equal to ${\gamma }_3={\gamma }_f=0.5$ for consistency with the assumptions of Equation 20 and CY14. That left the hinge distances, $r_1$ and $r_2$ , and the spreading rates, ${\gamma }_1$ and ${\gamma }_2$ , to be solved to fully specify Equation 25. The other key difference for this intermediate inversion was that the optimal model used $R_{\mathrm{RUP}}$ within the anelastic filter, while the alternative parameterization moved to use $R_{\mathrm{PS}}$ throughout the entire path-scaling model. The inversion used Equation 11 and sought to minimize the differences in predictions arising from the use of Equation 25 rather than Equation 20. The results of this inversion are shown in Figure 7 and resulted in the identification of hinge distances of $r_1=25$ km and $r_2=85$ km.

Figure 7.

Comparison of the geometric spreading functions of Equations 20 and 25 (“CY14mod target” shows the scaling of Equation 20, while “Piecewise result” corresponds to Equation 25). (a) The spreading function for various magnitudes against the equivalent point-source distance is shown, while (b) shows the same function plotted against the rupture distance. Vertical lines in both panels correspond to either hinge distances, $r_1$ and $r_2$ , or the transition distance $r_t$ . These can be magnitude dependent, depending on whether $R_{\mathrm{PS}}$ or $R_{\mathrm{RUP}}$ is used on the abscissa.

It is not possible to obtain perfect agreement between the two path-scaling variants. A key challenge is the mixing of distance metrics, with hinge distances being defined in terms of $R_{\mathrm{PS}}$ for Equation 25, and the transition distance $r_t$ relating to $R_{\mathrm{RUP}}$ in Equation 20. The magnitude dependence of $R_{\mathrm{PS}}$ dictates that differences observed in Figure 7 depend on magnitude. However, in addition to the geometric spreading parameterization differing, the anelastic attenuation function operates upon $R_{\mathrm{RUP}}$ under the optimal model and $R_{\mathrm{PS}}$ for the convenience model. The apparently subtle change in anelastic attenuation actually has a non-trivial influence for the large-magnitude scenarios, and we see this impact in Figure 7b, in particular. For large magnitudes, even short $R_{\mathrm{RUP}}$ distances map to non-trivial $R_{\mathrm{PS}}$ values, and lead to anelastic effects being relevant.

With the hinge distances defined at $r_1=25$ and $r_2=85$ km, a full inversion is then conducted with ${\gamma }_1$ and ${\gamma }_2$ free to adopt new values, and all other FAS parameters free to vary. All the practical details (the parameter space, the numerical algorithms, etc) associated with this inversion are the same as for the inversion to find the optimal parameterization. The results of this inversion are presented in Table 2.

Table 2.

Parameter estimates and standard errors for the convenience parameterization (geometric spreading defined according to Equation 25, magnitude-independent $\eta$ , and $R_{\mathrm{PS}}$ used throughout)

Parameter	Estimate	Standard error	Parameter	Estimate	Standard error
$s_{\alpha }$	2.767	0.0297	$h_{\beta }$	0.4768	0.0668
$s_{\beta }$	0.6451	0.0317	$h_{\gamma }$	1.1513	0.0
$s_{\gamma }$	0.4077	0.0142	$h_{\delta }$	3.418	0.779
$s_{\delta }$	0.117	0.0173	$h_{\epsilon }$	7.088	0.118
${\gamma }_1$	1.1680	0.00644	$Q_0$	183.7	4.69
${\gamma }_2$	0.9293	0.0214	$\eta$	0.7077	0.0107
$h_{\alpha }$	−0.7771	0.541

Note that $h_{\gamma }$ is again fixed in the inversion, and that the $s_{\alpha }$ parameter leads to stress parameter values in units of MPa.

Model performance

The performance of both the optimal and convenience models are demonstrated through visual comparison with CY14 in Figures 8 to 10. In all cases, the dashed lines represent the target CY14 model predictions, the solid lines correspond to the optimal FAS model, and the dotted lines are for the convenience model.

Figure 8.

Comparison of magnitude scaling in the CY14 GMM and the two RVT-based models defined by the parameters in Tables 1 (“Optimal RVT”) and 2 (“Convenience RVT”). (a)–(f) Each show the scaling for a particular period. The legends in (a) and (f) apply to all panels.

Figure 9.

Comparison of distance scaling in the CY14 GMM and the two RVT-based models defined by the parameters in Tables 1 (“Optimal RVT”) and 2 (“Convenience RVT”). (a)–(f) Each show the scaling for a particular period. The legends in (a) and (f) apply to all panels.

Figure 10.

Comparison of spectral scaling in the CY14 GMM and the two RVT-based models defined by the parameters in Tables 1 (“Optimal RVT”) and 2 (“Convenience RVT”). (a)–(f) each show the scaling for a particular magnitude. The legends in (a) and (f) apply to all panels.

Figure 8 shows the magnitude scaling of the models for six periods that span the range of typical engineering interest. The most obvious feature of interest in this figure is the difference in magnitude scaling at the very short distance of $R_{\mathrm{JB}}=1$ km and for magnitudes in the approximate interval (5,7). This behavior is primarily a consequence of differences in near-source saturation models over this magnitude range. The nonlinear constraint imposed within the inversions to prevent oversaturation works for the larger magnitudes (Stafford, 2021a), but differences in $\partial h\left(\mathit{M}\right)/\partial \mathit{M}$ (that can be appreciated from Figure 5d) can lead to the observed scaling at $R_{\mathrm{JB}}=1$ km. However, it can also be appreciated that this feature disappears very quickly, and is hard to perceive at $R_{\mathrm{JB}}=10$ km.

For all other periods and distances, we generally observe a close agreement between all three sets of model predictions. In judging the performance of the RVT-based predictions, it is important to keep in mind that the primary use of the RVT parameterization is to develop host-to-target model adjustments. In this context, absolute differences in amplitude are arguably less important than achieving a good similarity of the “shapes” of the curves. Furthermore, “local” issues like the magnitude scaling at $R_{\mathrm{JB}}=1$ km are very rarely of practical interest when deriving model adjustments (e.g. those scenarios need not be considered when developing the adjustment factors). The ability of the RVT models to largely replicate the scaling of CY14 is a direct consequence of that model basing its magnitude scaling on Fourier spectral theory. This was one of the reasons cited by Bommer and Stafford (2020) for recommending CY14 as a robust backbone GMM.

It is also noteworthy that the loss in performance associated with using the convenience parameterization over the optimal parameterization is not visually obvious from Figures 8 to 10. Figure 9, focusing upon distance scaling, should be expected to highlight differences between these two RVT-based models most clearly. However, there are many rupture scenarios (often at relatively short distances) where the convenience parameterization outperforms the optimal model. Where the optimal model tends to show that superior performance is at large distances and relatively short periods. For these scenarios, the use of magnitude-dependent $\eta \left(\mathit{M}\right)$ helps the optimal model, and the change of the geometric spreading formulation is less significant. Although these large distance scenarios are rarely important for hazard studies, they frequently coincide with scenarios for which empirical data are available for the target region. Consequently, improving the performance for these scenarios can have important practical consequences.

Part of the reason why the change in geometric spreading formulation is not so crucial is that despite Equation 20 attempting to mimic the functional form of CY14, the overall path scaling of $\mathrm{Sa}$ also depends upon the adopted scaling within the Boore and Thompson (2014, 2015) path duration model. Some of the deviations in spreading rates observed in Figure 9 for either RVT-based model or the CY14 model can be attributed to the scaling of the duration model, or at least its interaction with the geometric spreading functions in the FAS formulations.

Figure 9 shows predictions obtained using $\mathrm{\Delta }{Z}_{\mathrm{TOR}}=0$ in all cases. The shortest rupture depths for the smallest magnitudes are, therefore, influenced by the expected $Z_{\mathrm{TOR}}$ model shown in Figure 4. If we allow very shallow depths to be considered, differences between the RVT and CY14 predictions arise as a result of $h\left(\mathit{M}\right)$ effects. Specifically, the FAS parameters in Tables 1 and 2 both indicate very short saturation distances at small magnitudes, while CY14 has quite large $h\left(\mathit{M}\right)$ , as seen in Figure 5d. The relative decrease in $\mathrm{\Delta }\sigma$ associated with shallower-than-expected ruptures, does not fully compensate for the differences in $R_{\mathrm{PS}}$ between the RVT models and CY14. However, for the purposes of developing host-to-target model adjustments, these scenarios are rarely of importance.

The comparison of spectral scaling in Figure 10 shows that we achieve very good agreement in both spectral shape and absolute amplitude for the small-to-moderate magnitude scenarios. This is largely to be expected given that the point-source approximation should work most effectively for these magnitudes. In Figure 10e and f, we see that the RVT-based models have difficulty replicating the spectral shape of CY14. In particular, at long periods, the rate of spectral decay with period is stronger for CY14 than in the RVT models. This portion of the response spectrum is primary controlled by low frequencies in the Fourier domain where our FAS model has relatively few degrees of freedom. In particular, once a single-corner source spectrum is adopted, we have very limited control over these low-frequency FAS ordinates. At short periods and large distances, we also have challenges replicating the spectral shape for large magnitudes. Here, the role of the magnitude-dependent $\eta \left(\mathit{M}\right)$ is most clearly appreciated. For the large magnitudes, it is clear that the magnitude-independent $\eta$ within the convenience parameterization systematically under-predicts the amplitudes at short periods and over-predicts at long periods. While this is largely true of the optimal model as well, it can be seen that the optimal model is able to mimic the target shape slightly better, that is, the spectral peak is slightly better aligned with the target.

The degree of agreement between the CY14 target and the RVT-based predictions shown in Figures 8 to 10 can be slightly misleading as a result of plotting predictions over multiple decades. For this reason, Figures 11 to 13 show the relative performance of the optimal and convenience models through ratios with respect to the CY14 target. In all cases, values of unity correspond to a perfect match with the CY14 prediction. While these figures demonstrate that the RVT-based models are generally within 50% of the CY14 target, it is clear that there are local variations in the degree of agreement with respect to magnitude. To put this into some perspective, Al Atik and Youngs (2014) showed that the model-to-model standard deviation among the NGA-West2 GMMs is on the order of 0.2 natural log units—corresponding to a factor of ≈20% (with the model-to-model range being larger). For the distance ratios in Figure 12, the ratios are largely flat to beyond 100 km, before the differences in anelastic attenuation play a stronger role. For the spectral ratios, Figure 13 shows that for most distances, the ratios are also relatively flat at short periods (indicating similar spectral shape to CY14), but that departures increase at longer periods.

Figure 11.

Ratios of the RVT-based predictions to the CY14 target with respect to magnitude, for a number of distances and periods.

Figure 12.

Ratios of the RVT-based predictions to the CY14 target with respect to distance, for a number of magnitudes and periods.

Figure 13.

Ratios of the RVT-based predictions to the CY14 target with respect to period, for a number of magnitudes and distances.

In general, consideration of Figures 11 to 13 indicate that both the optimal and convenience parameterizations face similar challenges in replicating the scaling of the CY14 model. Comparisons of corresponding pairs of curves in these figures tends to show that the optimal model (solid lines) is slightly closer to unity than the convenience model (dotted lines), on average. However, it is also clear from the complete set of Figures 8 to 13 that the optimal parameterization is not drastically superior to the more familiar convenience parameterization.

As noted previously, the primary purpose of these inversions is to facilitate host-to-target model adjustments, and it is important to keep in mind where, in magnitude–distance space, empirical constraint for the target region tends to come from. Generally speaking, for the small-to-moderate magnitudes and intermediate-to-long distances of greatest practical importance, both the FAS parameterizations do a good job of replicating the absolute amplitude and scaling of the CY14 GMM.

Discussion and conclusion

In this article, we have presented host-region parameters that are broadly consistent with the CY14 GMM, with the intention of facilitating the application of a backbone GMM approach to constructing a logic tree for use in PSHA in a region of crustal seismicity. The CY14 accounts for scaling features that cannot be captured through standard parameterizations of the Fourier spectrum based on point-source seismological theory, and it is never the intention to perfectly match the CY14 predictions using the inverted FAS parameter set. However, it is important to achieve as close a match as possible so that model adjustments derived within the HEM framework do not distort the general scaling of the CY14 model in undesirable ways when it is adapted for application in the target region. In addition, it is important to capture magnitude dependencies of the FAS parameters so that parametric comparisons between the host model and the target region relate to similar rupture scenarios.

If the approach and these parameter suites presented in the article are accepted, then the task of developing a ground-motion logic tree can focus exclusively on determination of the target-region parameters. This is likely to be achieved through a combination of in situ site measurements and analyses of ground-motion recordings from the target region, including estimation of the uncertainty associated with each of the target-region parameters. The source and path differences can be accommodated through direct adjustments to the appropriate coefficients of CY14, or by developing adjustment factors using the HEM framework. The site adjustment can be made in one of two ways, either through the application of $V_S-{\kappa }_0$ adjustments (e.g. Al Atik et al., 2014) or through applying site response analyses to determine amplification factors for both the host and target profiles (with ${\kappa }_0$ distributed as damping in the profiles) and then applying the ratios of these amplification factors to the backbone GMM. If the resulting range of uncertainty is judged to be insufficient, for example, by comparisons with the ranges implicit in suites of high-quality GMMs such as those developed for NGA-West2 (Gregor et al., 2014) or NGA-East (Goulet et al., 2018), then additional nodes or branches can be incorporated within the logic tree. The comparison should not be made with the range of predicted amplitudes but rather with the range of uncertainty (i.e. model-to-model variability), the philosophy being that in a data-poor target region, the epistemic uncertainty should be at least as large as that defined for ground-motion logic trees in the data-rich host regions such as California (Al Atik and Youngs, 2014).