Methods to account for shallow site effects in hybrid broadband ground-motion simulations

Overview

Figure 1a presents five general approaches (Approaches A–E) that can be used to model site effects in a broad range of applications (e.g. empirical or physics-based ground-motion modeling), and Figure 1b illustrates their relationship with the five methods (Methods 1–5) examined in this article, which are specifically tailored to adjusting hybrid broadband ground-motion simulations. Approaches A–E are discussed below, while the site-adjustment Methods 1–5 are described in the “Methods for performing site adjustments to hybrid broadband ground-motion simulation” section.

OPEN IN VIEWER

Figure 1.

(a) Time-domain and frequency-domain approaches to model site effects, and their respective input data. (b) The relationship between these approaches and the investigated methods to perform site adjustments to hybrid broadband ground-motion simulations. L: linear; NL: nonlinear; EQL: equivalent-linear.

As illustrated in Figure 1a, going from Approach A to E implies the use of incrementally less site information and the introduction of some simplifications in the description of the site conditions and wave-propagation phenomena. Table 1 provides a comparison between these five approaches in terms of their relative advantages and disadvantages, and further discussion is provided below.

OPEN IN VIEWER

Table 1.

Comparison of different approaches to model site effects

Approach	Pros	Cons
A. 3D/2D time-domainsite-response analysis	3D/2D wave-propagation effects and soil nonlinearity can be explicitly modeled	Requires more site-characterization data than is usually availableSignificant computational cost and analyst competency
B. 1D time-domain site-responseanalysis	1D wave-propagation effects and soil nonlinearity can be explicitly modeledRequired site-characterization data can be collected at a reasonable cost compared to Approach A	Performance diminishes when the 1D SH wave propagation (SH1D) assumptions do not holdAdvanced constitutive models require the calibration of multiple parameters
C. 1D frequency-domain site-responseanalysis	Linear 1D wave-propagation effects are explicitly modeledRequired site-characterization data can be collected at areasonable cost compared to Approach A	Performance diminishes when the SH1D assumptions do not holdLimitations to model soil nonlinearity
D. Square-root-impedance method	Site-specific impedance and attenuation effects are modeledRequired site-characterization data can be collected at areasonable cost compared to Approach AComputed siteamplification is not sensitive to the details in the shear-wavevelocity ( V S ) profile (which is generally characterized by considerable levels of uncertainty)	It can underpredict site response when resonance effects are relevantIt is not able to model soil nonlinearity
E. Site response component of ground-motion models	3D wave-propagation effects (consistent with the sites considered in the GMM development) are implicitly capturedRequires considerably less site-characterization data andanalyst competency than Approaches A–D	Characterization via proxy parameters (e.g., V S 30 ) may ignore site-specific featuresBiased predictions if appropriate host-to-target adjustment is not considered

Approach A: 3D/2D time-domain site-response analysis

In this approach, the 3D or 2D wave equation is solved numerically (e.g. using finite element, finite difference, or spectral element methods) at different time steps, and soil nonlinearity is generally simulated using plasticity models (e.g. Yang et al., 2003). This method is the only one with the capability of explicitly simulating complex 3D/2D wave-propagation features such as basin effects (e.g. Ayoubi et al., 2021) or wave scattering due to soil heterogeneities (e.g. de la Torre et al., 2022a; Thompson et al., 2009), offering the highest potential to realistically simulate site effects among the approaches presented in Figure 1a. However, this would involve significant challenges in the model characterization, such as data collection within a great lateral extension (e.g. Hallal and Cox, 2023) or the definition of random-field (e.g. de la Torre et al., 2022b) and spatially distributed nonlinear parameters. Furthermore, defining an input motion on the model boundaries that adequately reflects the spatial variability of the incident wavefield can be challenging. This approach is not considered in the methods examined in this study to adjust ground-motion simulations.

Approach B: 1D time-domain site-response analysis

The numerical simulation in this approach relies on the 1D wave-propagation assumptions: the medium is represented by laterally continuous and homogeneous layers overlying a half-space; wavefronts are considered to be planar; and only horizontally polarized shear waves propagating in the vertical direction are modeled (the so-called SH1D assumptions). In addition to plasticity models, simpler cyclic stress-strain relationships (e.g. Groholski et al., 2016; Shi and Asimaki, 2017) can be used to capture soil nonlinearity. Adopting the SH1D assumptions has two implications: (1) it allows the problem to be represented with a relatively simple site-response model, and using site-characterization data (e.g. geophysical and cone penetration tests [CPTs]) that are typically collected within a relatively small area, corresponding to the site of interest, and (2) sites with conditions that violate these assumptions can result in inaccurate ground-motion predictions, as shown by several validation studies (e.g. Afshari and Stewart, 2019; Kaklamanos et al., 2013; Pilz and Cotton, 2019; Thompson et al., 2009, 2012; Zhu et al., 2022). As illustrated in Figure 1b, this approach is used in Method 5.

Approach C: 1D frequency-domain site-response analysis

This approach also relies on the SH1D assumptions. However, in this case, the problem is solved in the frequency domain through the direct application of a complex-valued transfer function, T F SH 1 D , to the Fourier transform of the input motion (ground motion at a reference condition, that needs to be propagated through the soil column). Alternative methods are available to compute T F SH 1 D for a given configuration of horizontal soil layers over an elastic half-space (e.g. Haskell, 1953; Kramer, 1996; Thomson, 1950). This complex transfer function (TF) accounts for both the modification of the amplitude and the phase of the input motion. The equivalent-linear (EQL) method (Seed and Idriss, 1969) is typically used to consider soil nonlinearity if necessary. This method has been widely adopted in engineering practice due to its simplicity but presents some limitations such as the overdamping (i.e. underprediction) of ground motion at high frequencies when the soil experiences relatively large shear strains (e.g. Kaklamanos et al., 2013; Kim et al., 2016). As shown in Figure 1b, this approach is used in Methods 4 and 5.

Approach D: square-root-impedance method

This seismological approach combines the square-root-impedance (SRI) method (Boore, 2003, 2013; Joyner et al., 1981) to model site amplification, I ( f ) , with a diminution operator, D ( f ) , to capture site attenuation (Anderson and Hough, 1984). The resulting “outcrop” TF, T F SRI ( f ) , can be applied in the same way as in Approach C:

T F SRI ( f ) = I ( f ) D ( f ) = ρ R V S , R ρ ¯ V S ¯ exp ( − π f κ 0 )

(1)

where ρ R and V S , R are the density ( ρ ) and shear-wave velocity ( V S ), respectively, at the reference condition; ρ ¯ and V S ¯ are the corresponding average ρ and travel-time weighted average V S , respectively, computed at each frequency, f , considering a depth equivalent to the quarter-wavelength; and κ 0 is the site attenuation parameter (Anderson and Hough, 1984; Ktenidou et al., 2014) estimated at the surface. The diminution operator can also account for a portion of κ 0 , corresponding to the contribution of the near-surface materials down to a reference horizon (e.g. Cabas and Rodriguez-Marek, 2017), in which case, κ 0 is replaced by Δ κ 0 in Equation 1, where Δ κ 0 = κ 0 − κ 0 , R (being κ 0 , R the value estimated at the reference horizon). Approach D does not account for soil nonlinearity.

Several differences exist between T F SRI (Approach D) and T F SH 1 D (Approach C): (1) T F SRI does not consider any phase adjustment, which is accounted for in T F SH 1 D ; (2) the amplification term in T F SRI represents only impedance effects and does not account for resonance effects (Thompson et al., 2011; Zhu et al., 2022), which are explicitly modeled in T F SH 1 D ; (3) the attenuation model is conceptually different, resulting in lower levels of damping for T F SRI in the case of a one-layer system when equivalent attenuation parameters are used for T F SH 1 D and T F SRI (Boore, 2013); and (4) T F SRI is significantly less sensitive to details in the V S profile than T F SH 1 D (Boore, 2013). Due to aspect (2), T F SH 1 D produces peaks and troughs in the TF, with significant amplification over narrow frequency bands, which are absent in T F SRI . The existence of these resonant features strongly depends on the suitability of the SH1D assumptions and the considered horizontal stratigraphy, as well as stiffness and damping properties. Soil heterogeneity and other non-1D wave phenomena can prevent, or significantly reduce, their development (e.g. de la Torre et al., 2022a; Thompson et al., 2011). Due to several sources of uncertainty that exist in the estimation of the V S profile (Rodriguez-Marek et al., 2021), aspect (4) means that important levels of uncertainty can be also expected in the computation of T F SH 1 D (Thompson et al., 2011). Boore and Abrahamson (2023) proposed a modification to the SRI approach that gives results closer to T F SH 1 D for gradient profiles, but it is not explored in this article. As illustrated in Figure 1b, this Approach D is used in Methods 2, 3, 4, and 5.

Approach E: site response component of a GMM

Semi-empirical GMMs (Baker et al., 2021, Ch. 4) provide a prediction of the mean value of a given intensity measure (IM) in natural log scale, for an earthquake rupture, rup , and the site considered, site , with the general form,

μ ln IM ( rup , site ) = f E + f P + f S

(2)

where f E , f P , and f S are source, path, and site response scaling factors, respectively. In the context of this article, two IMs will be primarily considered: pseudo-spectral acceleration (SA) and Fourier amplitude spectra (FAS). The site response scaling factor ( f S ) can be generally expressed as the sum of three terms (e.g. Bayless and Abrahamson, 2019; Campbell and Bozorgnia, 2014):

f S = f L + f NL + f sed

(3)

where f L and f NL represent the scaling of the linear and nonlinear site response, respectively, with one or more site parameters (e.g. V S 30 ). In addition, the term f NL is a function of some IM at a reference condition ( I M rock ). For example, the CB14 model (Campbell and Bozorgnia, 2014) uses the estimated median value of peak ground acceleration (PGA) on rock with V_S30 = 1100 m/s. The term f sed represents the scaling with the thickness of the site column, with V S lower than a certain threshold value. This is typically characterized by the parameter Z 1 . 0 , corresponding to the depth to a Vs of 1.0 km/s (e.g. Bayless and Abrahamson, 2019; Boore et al., 2014), although some GMMs (e.g. Campbell and Bozorgnia, 2014) use Z 2 . 5 (depth to 2.5 km/s) instead. The f sed term can also depend on V S 30 when the scaling factor is based on the relationship between Z 1 . 0 and the average Z 1 . 0 expected for a given region computed using a V S 30 -based correlation (e.g. Bayless and Abrahamson, 2019; Boore et al., 2014).

The terms f L , f NL , and f sed are based on observational ground-motion data, theoretical models, and numerical simulations. In particular, due to the lack of recordings with strong soil nonlinearity, several models have been developed for f NL using a large number of 1D site-response analyses (e.g. Hashash et al., 2018; Kamai et al., 2014; Shi, 2019; Walling et al., 2008) or combining these analyses with observations (e.g. Seyhan and Stewart, 2014), and they are implemented in different GMMs. Most of these models represent only an adjustment to the amplitude of the ground motion in the Fourier or response spectral domain, but a phase adjustment can also be included (e.g. Shi, 2019: Ch.4). As shown in Figure 1b, this approach is used in Methods 1, 2, 3, and 4.

Overview

Currently, 3D numerical ground-motion simulations of the earthquake source and (3D) wave propagation are utilized in engineering applications only up to a limited maximum frequency of f  = 1–3 Hz (Graves and Pitarka, 2016; Paolucci et al., 2021). This is primarily due to knowledge limitations regarding small-scale features of the fault rupture and the Earth’s crust (Graves and Pitarka, 2016; Mai et al., 2010) and is compounded by the large computational demands of simulating small length scales. As a result, hybrid simulation methods have been applied (e.g. Graves and Pitarka, 2010; Mai et al., 2010), which utilize this comprehensive 3D physics-based approach for the LF simulation and a simplified physics-based approach for the HF simulation. The resulting LF and HF waveforms are then combined using a matched filter around a transition frequency, f t . The transition frequency f t often corresponds to 1 Hz but generally is increasing over time as computation, theory, and velocity models improve.

The method developed by Graves and Pitarka (2010, 2015, 2016) (herein referred to as the GP method for brevity) has been applied and validated in different regions, including California, Europe, and New Zealand (e.g. Dreger et al., 2015; Lee et al., 2022; Razafindrakoto et al., 2021) and is considered in this article to illustrate the challenges in the modeling of shallow site effects in hybrid broadband simulations. In this method, the LF component is simulated using a kinematic description of the earthquake fault rupture and a visco-elastic finite-difference algorithm to solve the 3D wave propagation equation. The HF component is simulated using a stochastic representation of source radiation and a simplified Green’s function determined for a given 1D velocity structure. In particular, a frequency-dependent quality factor ( Q ) is used to incorporate path anelastic attenuation, and the SRI method along with the multiplicative operator P ( f ) = exp ( − π f κ 0 , sim ) (Approach D) is used to represent impedance effects and capture HF spectral decay, where κ 0 , sim is the site attenuation parameter considered in the simulation. Further details on the implementation of this method in this study are provided in Razafindrakoto et al. (2018) and Lee et al. (2022).

Considerations for the modeling of shallow site effects

The modeling of shallow site effects in hybrid broadband ground-motion simulations is usually performed using the “uncoupled approach” (e.g. de la Torre et al., 2020; Graves and Pitarka, 2010; Pilz et al., 2021; Razafindrakoto et al., 2021), and therefore, it takes the form of a posterior adjustment to the ground-motion time series produced by the regional-scale simulation. The four main aspects that need to be considered when performing this adjustment are discussed below.

First, the hybrid nature of the simulation method requires different considerations for the adjustment of the LF and HF components. For example, in the case of the GP method, the 3D velocity structure considered in the LF approach allows for the incorporation of any level of spatial variability within the simulation domain, whereas the HF approach is usually performed using a unique 1D velocity profile representative of the whole region of interest (e.g. Graves and Pitarka, 2010; Lee et al., 2022). Consequently, the 1D velocity profile implicit in the LF and HF components of the simulated ground motions at a given site can be significantly different, including the minimum V S ( V S , min ).

Second, due to knowledge and computational limitations, the minimum V S considered in the regional-scale simulation typically corresponds to a stiff soil or rock condition (e.g. V S , min  = 500 m/s in Lee et al., 2022), and the spatial variations of V S are represented with a relatively coarse resolution (e.g. a constant grid spacing of 100 m was considered for the LF simulation in Lee et al., 2022). In reality, near-surface materials typically reach lower values and exhibit significant spatial variation in V S within the top 100 m. Thus, site effects adjustments should be able to capture these features not modeled in the regional simulation.

Third, although near-surface sediments are generally not properly modeled in the regional-scale simulation (e.g. due to the high V S , min value and coarse grid spacing in Lee et al., 2022), local site effects may be captured to some extent. For instance, “basin effects” (Ayoubi et al., 2021) may be dominated by the deep velocity structure of the site, which can be relatively well represented in the LF component of the simulation (e.g. a relatively well-constrained 3D crustal velocity model of the Canterbury Region is considered in the LF simulation conducted by Lee et al., 2022), even if the shallower sedimentary layers are coarsely represented. Another example is related to the HF site attenuation, which is modeled in the HF component of the simulation through the parameter κ 0 (e.g. κ 0 , sim  = 0.045 s was used in Lee et al., 2022). Both “basin effects” and HF attenuation are already incorporated in the regional simulation to some degree, which must be taken into account when applying the site response adjustment, to avoid double-counting them.

Finally, regional-scale simulations are usually performed considering linear viscoelastic materials (e.g. Graves and Pitarka, 2010; Lee et al., 2022). Hence, the adjustment has to model soil nonlinearity, which is particularly relevant for maximum shear strains γ max  >  0 . 01 % – 0 . 1 % (e.g. Kaklamanos et al., 2013). While several 3D ground-motion simulations have explicitly included non-linear near-surface response (e.g. Fu et al., 2017; Taborda et al., 2012), such simulation results are still limited by the minimum spatial scale challenges alluded to above, as well as the ability to represent site-specific near-surface stratigraphy a priori in the regional simulation (as opposed to a posterior adjustment via an uncoupled approach which allows for site-specific data collection and site-response analysis for regional simulations that already have been performed).

General description

A case study is considered to examine the different methods to adjust ground-motion simulations and compare and contrast their features. It includes four strong motion station (SMS) sites in Christchurch, New Zealand, for which site characterization data (e.g. Teague et al., 2018; Wotherspoon et al., 2015) and simulation results of historical events (Lee et al., 2022; Razafindrakoto et al., 2018) are available. Table 2 shows that the four sites represent a wide range of site conditions. V S 30 , actual and Z 1 . 0 , actual are the measured or inferred values from the preferred V S profile available at each location. In the case of PRPC and CACS, the measured V S profile was relatively shallow, and it did not reach 1000 m/s, so the New Zealand Velocity Model (NZVM) (Thomson et al., 2020), version 2.07, was used to estimate Z 1 . 0 , actual . Z G is the depth to the first gravel formation, which in the central and eastern part of Christchurch corresponds to the Riccarton Gravel Formation and represents a significant velocity contrast relative to the softer sediments above (Teague et al., 2018; Wotherspoon et al., 2015). Based on these parameters, PRPC and CMHS represent relatively soft sites; with CMHS a shallower soil deposit (i.e. smaller Z G and Z 1 . 0 , actual values), consistent with its location near the basin edge (see Supplemental Figure A.1 in Electronic Supplement A). RHSC and CACS are stiffer sites, with CACS an example where the gravel formation starts at the surface.

OPEN IN VIEWER

Table 2.

Parameters representative of the actual and simulation site conditions for the four sites considered

Site	Actual site conditions			Simulation site conditions
	V S 30 , actual (m/s)	Z G (m)	Z 1 . 0 , actual (m)	V S 30 , sim (LF and HF) (m/s)	Z 1 . 0 , simLF (m)	Z 1 . 0 , simHF (m)
PRPC	196	28	473	500	500	600
CMHS	203	14	57	500	100	600
RHSC	294	6.5	555	500	400	600
CACS	435	0.0	557	500	600 ( * )	600

Z G is the depth to the first gravel formation. The Z 1 . 0 values of the LF and HF simulations differ due to the differences in the 3D and 1D velocity structures used.

( * ) The V S value actually reached at the 600-m depth is 0.98 km/s, which remains constant down to 1200-m depth, where V S increases to 1.15 km/s. Because of this particular stratigraphy, 600 m was chosen as a more meaningful value for Z 1 . 0 , simLF .

A low- and a high-intensity ground motion are considered for each site in order to examine the different methods in the near-linear and nonlinear response regimes. The associated events are described in Electronic Supplement A. Only one horizontal component of these ground motions is used.

Site-characterization data

Figure 2 presents the V S profile (via surface-wave testing) and CPT data used to characterize the four sites considered, down to 30 m. As observed in the soil behavior type index ( I c ) plot, in the top 28 m, the site PRPC is mainly characterized by clean and silty sands, with some localized layers of more cohesive materials; whereas in the top 13.8 m, the site CMHS has a more variable stratigraphy, with a greater presence of cohesive soils. For RHSC and CACS, CPT data were not available, but boreholes and standard penetration test (SPT) data were considered. Further details and site data available at these four locations are provided in Wotherspoon et al. (2015), Teague et al. (2018), and Electronic Supplement A.

OPEN IN VIEWER

Figure 2.

V S profiles and CPT data used for the four sites considered.

Ground-motion simulations considered

The events considered in this case study were simulated using the GP method. The simulations of the low-intensity events were performed by Lee et al. (2022), and those of the high-intensity event by Razafindrakoto et al. (2018). For the LF component, the NZVM (Thomson et al., 2020) v2.02 with a grid spacing of 100 m and a minimum V S  = 500 m/s was utilized. For the HF component, the nationwide 1D velocity model described in Lee et al. (2022) was used, which also considers a minimum V S  = 500 m/s. The transition frequency between LF and HF components was set to f t  = 1 Hz.

Figure 3a shows the 1D V S profile extracted at the site PRPC from the 3D velocity model used in the LF simulation, and the corresponding HF V S profile (herein, the LF and HF simulation profiles), illustrating that the two profiles have different deep velocity structures. The LF profile captures a localized volcanic formation that produces strong V S contrasts at depths of 500 m and 1000 m. These velocity contrasts do not appear in the HF profile, which is a regional average representation. Figure 3b illustrates the common case where the regional-scale LF and HF simulations are not able to accurately capture the V S of the surficial soil, which is typically characterized by V S  < 500 m/s with significant variation at a scale of meters. Table 2 provides the site parameters V S 30 and Z 1 . 0 associated with the simulation conditions for the LF and HF components. Similar figures for all the sites are provided in Electronic Supplement A.

OPEN IN VIEWER

Figure 3.

Site-specific, LF simulation, HF simulation, and NZVM (v2.07) V S profiles at PRPC. The horizontal dashed lines represent the grid spacing considered in the LF simulation: (a) 0–5000 m deep. (b) 0–220 m deep.

Overview

In this section, five methods to perform the site adjustment are presented, and the case study sites are used to illustrate their application and compare the methods under different soil conditions. As shown in Figure 1b, these methods are based on a range of different approaches discussed in the “Approaches for modeling local site effects” section; in particular, Approaches B–E. Approach A was not considered in this study because its use in engineering applications is limited. For example, even though the sites considered in the case study are well-characterized near the SMS, proper implementation of Approach A would require spatially distributed data over a considerably larger area (Hallal and Cox, 2023), and these data are significantly sparser. Also, defining an input motion that is compatible with the regional-scale ground-motion simulation is challenging. While methods exist that allow determining the equivalent input forces derived from 3D ground-motion simulations (e.g. Bielak, 2003), their implementation requires significant analyst competency and dealing with the hybrid nature of the simulation method.

The site adjustment can be performed in the frequency domain (Methods 1–4) or in the time domain (Method 5). The “frequency-domain adjustment” involves five steps:

Since the simulated acceleration time series are different for the LF and HF components, the adjustment is performed separately for each component, and then the adjusted LF and HF time series are combined in the time domain. This procedure is described in the “Application of the site factor” section of Electronic Supplement B. Given that the simulation profiles are also different for each component, this may result in different LF and HF SFs. In the following sections describing Methods 1–4, the LF and HF SFs are only plotted in their associated frequency ranges ( f < 1 Hz and f  > 1 Hz, respectively), but it should be noted that when they are applied to the simulated ground motions, their influence slightly exceed these frequency ranges in the vicinity of f  = 1 Hz, due to the low- and high-pass filters applied before merging, which ensure a smooth transition. Also, due to the differences alluded to earlier, the LF and HF SFs may not perfectly match at f  = 1 Hz, which manifests as a discontinuity in some of the figures presenting the SFs.

The “time-domain adjustment” involves the use of time-domain site response analysis. Method 5 corresponds to a specific implementation of this approach, based on 1D nonlinear inelastic site-response analysis, and it is further explained in the “Method 5 - 1D time-domain wave propagation analysis” section.

The site PRPC is first used to explain the methods, but results for all four case study sites are provided in Electronic Supplement C and in the “Comparison of the site amplification obtained with the five methods” section.

Method 1— V S 30 -based amplification

This method uses the site response scaling factor ( f S ) from existing semi-empirical GMMs (Approach E) and corresponds to the method used in Graves and Pitarka (2010) to capture shallow site effects. In principle, the computation of the SF ( S F 1 ) is based on the ratio between the median prediction of the FAS for the actual site conditions ( FA S actual , and the corresponding prediction for the simulation site conditions ( FA S sim ):

S F 1 ( f ) = exp [ μ ln FA S actual ( f ) ] exp [ μ ln FA S sim ( f ) ]

(4)

Since semi-empirical models for FAS have only recently been developed (e.g. Bayless and Abrahamson, 2019; Bora et al., 2019), several previous studies have used GMMs for SA to approximately compute S F 1 instead (e.g. Graves and Pitarka, 2010; Lee et al., 2022). Due to the differences between SA and FAS (Bora et al., 2016), S F 1 estimates based on SA models have been truncated to a value of 1 at high frequencies to mitigate the inconsistencies between the two IMs (e.g. as discussed in de la Torre et al., 2020).

Since source and path effects are common to the numerator and denominator in Equation 4, it reduces to the ratio of site terms. Specifically:

S F 1 ( f ) = exp [ f S , actual ( f ) ] exp [ f S , sim ( f ) ] = exp ( f L , actual + f NL , actual + f sed , actual ) exp ( f L , sim + f NL , sim + f sed , sim )

(5)

where all terms have been previously defined in the “Approach E: Site response component of a ground-motion model” section. The terms f sed , actual and f sed , sim may or may not be included in Equation 5. For example, Graves and Pitarka (2010), who used the CB08 model (Campbell and Bozorgnia, 2008), did not include them, arguing that the effect of deep-basin amplification was already included in their LF simulation. If the regional ground-motion simulation is performed adopting a linear viscoelastic constitutive model, which is usually the case, the term f NL , sim should in principle be removed from Equation 5. This has been considered in some studies (e.g. Rodgers et al., 2020), but others (e.g. Lee et al., 2022) have not removed this term. Considering that usually f L ~ f ( V S 30 ) , f NL ~ f ( V S 30 , I M rock ) , and f sed ~ f ( V S 30 , Z 1 . 0 ) , Equation 5 can be rewritten as,

S F 1 ( f ) = exp [ f L ( f , V S 30 , actual ) + f NL ( f , V S 30 , actual , I M rock ) + f sed ( f , V S 30 , actual , Z 1 . 0 , actual ) ] exp [ f L ( f , V S 30 , sim ) + f sed ( f , V S 30 , sim , Z 1 . 0 , sim ) ]

(6)

This results in a SF that can be expressed as,

S F 1 ( f ) = S F 1 , L ( f ) · S F 1 , NL ( f )

(7)

where S F 1 , L = exp ( f L , actual + f sed , actual − f L , sim − f sed , sim ) ~ f ( V S 30 , actual , Z 1 . 0 , actual , V S 30 , sim , Z 1 . 0 , sim ) is the linear component of S F 1 , and S F 1 , NL = exp ( f NL , actual ) ~ f ( V S 30 , actual , I M rock ) is the nonlinear component. Several studies have truncated this SF to a value of S F 1  = 1 at low frequencies to avoid double-counting LF amplification already captured in the LF simulation component, based on a 3D velocity model (e.g. Graves and Pitarka, 2010; Lee et al., 2020, 2022). The I M rock parameter generally corresponds to the PGA at a certain reference V S 30 that varies between GMMs. Here, the PGA obtained from the HF simulation (HF PGA) at V S 30 , sim is used and scaled to the model-dependent reference V S 30 , applying a linear SF based on the corresponding GMM. This procedure is explained in more detail in Electronic Supplement B.

Table 3 summarizes some of the most recent studies on ground-motion simulations where Method 1 has been used to capture site effects. The table illustrates that there is no consistent approach to date to implement Method 1, and subjective decisions are involved in the selection of the GMM, the application of LF and HF truncations, and the inclusion of the f sed terms.

OPEN IN VIEWER

Table 3.

Recent ground-motion simulation studies that have used Method 1

Study	Simulation method	GMM used for computing SF(f)	Implementation
Shi (2019), Ch.4 Southern California, US	Hybrid (Graves and Pitarka, 2010) f t = 1 Hz , V s , min as implemented in the SCEC Broadband Platform (BBP) v17.3.0	Combination of BSSA14 (Boore et al., 2014) and CB14 (Campbell and Bozorgnia, 2014) IM: SASAG19 (Shi, 2019: Ch.4) IM: FAS	As implemented in the SCEC BBP v17.3.0Different from Equation 5
Lee et al. (2020) Canterbury, New Zealand	Hybrid (Graves and Pitarka, 2010) f t = 1 Hz , V S , min = 500 m / s	CB14 (Campbell and Bozorgnia, 2014) IM: SA	Not including f sed Truncation at LF and HF
Rodgers et al. (2020) Northern California, US	3D numerical simulation f max = 10 Hz , V S , min = 500 m / s	BA18 (Bayless and Abrahamson, 2019) IM: FAS	Including f sed Without truncations
Razafindrakoto et al. (2021) Upper Rhine Graben, Europe	Hybrid (Graves and Pitarka, 2010) f t = 1 Hz , V S , min = 800 m / s	CB14 (Campbell and Bozorgnia, 2014) IM: SABCS19 (Bora et al., 2019) IM: FAS	Not including f sed Truncation at LFand HFWithout truncations
Lee et al. (2022) New Zealand	Hybrid (Graves and Pitarka, 2010) f t = 1 Hz , V S , min = 500 m / s	CB14 (Campbell and Bozorgnia, 2014) IM: SA	Not including f sed Only applied to the HF component, with HF truncation

In order to illustrate the implications of alternative combinations of models, Figure 4 presents the resulting SFs ( S F 1 ) for the site PRPC. This comparison comprises two models developed in the response spectral domain (CB14, BSSA14) and one model developed in the Fourier spectral domain (BA18), which have been used in previous studies (Table 3). In the case of the BSSA14 and BA18 models, the results are presented with and without the f sed ( f , V S 30 , Z 1 . 0 ) terms ( f sed , actual and f sed , sim ), to illustrate the effect that this model component has on the SF. Further details of the models considered and additional results including the BCS19 and SAG19 models are provided in Electronic Supplements B and C, respectively.

OPEN IN VIEWER

Figure 4.

Site factors obtained for the site PRPC using Method 1 and different models: (a) and (b) Show the linear and nonlinear components for the low-amplitude and high-amplitude motion, respectively. (c) and (d) Present the resulting site factors without any truncations. (e) and (f) Show the site factors with LF and HF truncations.

Figures 4a and b display the linear ( S F 1 , L ) and nonlinear ( S F 1 , NL ) components of the SFs for the two events considered; Figures 4c and d present the resulting SFs ( S F 1 ) without any truncations; and Figures 4e and f show the SFs with two truncations applied: (1) for all the models, the LF SF is truncated to a value of S F 1  = 1 for f  < 0.2 Hz with a taper from f  = 0.5–0.2 Hz, as considered in some previous studies (e.g. de la Torre et al., 2020; Lee et al., 2020; Razafindrakoto et al., 2021); and (2) in the case of the SA-based models (CB14, BSSA14), the HF SF is truncated to a value of S F 1  = 1 for f  > 15 Hz, with a taper from f  = 10–15 Hz, as also considered in some previous studies (e.g. de la Torre et al., 2020; Lee et al., 2020, 2022; Razafindrakoto et al., 2021).

Significant model-to-model variability is observed in the linear and nonlinear components of the SF in Figures 4a and b, but in the latter case, this variability only manifests under the high-amplitude ground motion (Figure 4b). For the low-amplitude motion, the nonlinear component of the three GMMs is close to unity over the entire frequency range, whereas for the high-amplitude motion, this component substantially departs from 1 at frequencies greater than 0.5 Hz, resulting in de-amplification. These figures also show that the BA18 model accounts for the HF decay in the site response, a feature that cannot be properly captured in the response spectral domain (Bora et al., 2016; Cabas and Rodriguez-Marek, 2017).

A common feature observed in all the SFs, clearly illustrated in Figures 4c and d, is the relatively high amplification at low frequencies, particularly around f  = 0.3–1.0 Hz. This amplification is reduced when the f sed terms are included in the BSSA14 and BA18 models. This is mainly driven by an increase in the exp ( f sed , sim ) term at low frequencies, as explained in the “Effect of f sed on S F 1 . L ” section of Electronic Supplement C. Several validation studies conducted in New Zealand (de la Torre et al., 2020; Lee et al., 2020, 2022), for example, have shown that this high amplification may lead to significant overprediction in the LF range, even when an LF truncation is applied. Considering this, Lee et al. (2022) showed that the application of the SF to the HF component only (i.e. f  > 1 Hz) can significantly reduce the model bias at long vibration periods. One of the downsides of this approach is that the proper truncation may be site- and region-dependent.

Host-to-target adjustment issue

To investigate the root cause of the apparent overamplification at relatively low frequencies, Figure 5a illustrates a V S profile representing the actual site conditions (actual profile) at a hypothetical soil site, with an associated V S 30 , actual value. It also shows the corresponding simulation profile, which is assumed to be equivalent to the actual profile in the deeper portion (i.e. the deep velocity structure is well represented by the regional-scale simulation) but differs in the shallow portion because the regional-scale simulation is conducted truncating the V S profile to a V S , min value, due to computational and knowledge constraints (e.g. V S , min  = 500 m/s in the case study considered in this article). This results in V S 30 , sim  = V S , min . Under these conditions, the purpose of the SF should be to capture the site response generated by the actual profile relative to the simulation profile (due to the near-surface soil layers characterized by lower velocities in the actual profile) and accounting for soil nonlinearity. As discussed before, the determination of the site response in Method 1 is based on a GMM and one or more proxy parameters, such as V S 30 . For a given GMM and V S 30 value, there is a corresponding theoretical V S profile which reflects the implicit site conditions considered in the development of the GMM, which is refered to as the “host” V S profile (Al Atik and Abrahamson, 2021; Baker et al., 2021: Ch. 8). Figure 5a includes the host profiles associated with V S 30 , actual and V S 30 , sim for a hypothetical GMM.

OPEN IN VIEWER

Figure 5.

(a) Illustration of the host-to-target conversion issue present in Method 1. (b) Computation of the reference correction (RC) factor introduced in Method 2.

Figure 5a illustrates that a significant inconsistency may exist between the “target” condition to be modeled (i.e. relative site response between the actual and simulation profile) and the “host” condition implicit in the SF S F 1 (i.e. relative site response between the two host profiles). In the context of empirical ground-motion modeling and site-specific seismic hazard analysis, this problem is known as the host-to-target adjustment issue (e.g. Bard et al., 2020; Campbell, 2003; NRC, 2021; Williams and Abrahamson, 2021). In the specific case of Method 1, there are three main related issues:

Method 2— V S 30 -based amplification with host-to-target adjustment

In Method 2, a host-to-target adjustment is proposed to partially overcome the inconsistency involved in Method 1. Particularly, this adjustment addresses the first issue explained in the previous section (“Host-to-target adjustment issue”), related to the difference between the simulation profile and the host profile for V S 30 , sim . This method is inspired by previous studies that have suggested similar adjustments but in different contexts (e.g. Cabas and Rodriguez-Marek, 2017; NRC, 2021; Williams and Abrahamson, 2021). Method 2 is a V S correction approach based on the SRI method (Approach D), which allows for adjusting the Method 1 SF ( S F 1 ) to minimize the incompatibility between the simulation profile and site condition implicit in the GMM for V S 30 , sim . More specifically, the Method 2 SF ( S F 2 ) can be expressed as,

S F 2 ( f ) = exp [ f S , actual ( f ) ] exp [ f S , sim ( f ) ] · I host ( V S 30 , sim ) ( f ) I sim ( f ) = S F 1 ( f ) · RC ( f )

(8)

where I host ( V S 30 , sim ) and I sim are the SRI-based site amplifications for the host V S profile associated with V S 30 , sim and for the simulation profile, respectively; RC is the reference correction factor. The computation of this factor is illustrated in Figure 5b and is conducted using the following equation:

RC ( f ) = I host ( V S 30 , sim ) ( f ) I sim ( f ) = ρ R V S , R ρ ¯ host ( V S 30 , sim ) V ¯ S , host ( V S 30 , sim ) ρ R V S , R ρ ¯ sim V ¯ S , sim

(9)

where the subscripts host ( V S 30 , sim ) and sim indicate properties of the host profile for V S 30 , sim and of the simulation profile, respectively; and ρ R and V S , R are the density and V S , respectively, at a reference depth common to the two profiles, considered as the elastic half-space. Figure 5b indicates the depth at which the two profiles merge. The half-space can be located at this location or deeper. This reference correction factor intends to ensure consistency with the host profile for V S 30 , sim , but it does not address the potential inconsistency with the host profile for V S 30 , actual , which is a common limitation of any ergodic method based on V S 30 (for Methods 1 and 2, it is assumed that the only site information available is V S 30 , actual ).

The main challenge for the application of Method 2 is deriving the host profile implicit in the GMM(s) considered for V S 30 , sim . As discussed by Williams and Abrahamson (2021), there are two approaches to estimate this profile if the GMM developers do not provide it: (1) assuming that available measurements of V S within the region for which the GMM was developed are representative of the host profile (e.g. Kamai et al., 2016), and (2) deriving the profile from the GMM by solving the inverse problem (e.g. Al Atik and Abrahamson, 2021). In this study, the host profiles derived by the Al Atik and Abrahamson (2021) method for the CB14 and BSSA14 models are considered. Al Atik and Abrahamson (2021) showed that the profiles derived with their method (for four NGA-West2 GMMs) provide a better match with the site response implicit in the GMMs for very stiff site conditions (e.g. V S 30  = 620 m/s or greater), which are less influenced by 2D/3D effects. They showed that for lower V S 30 values, the method was less effective. The inversion was originally based on the SRI method, but the use of a modified version of the SRI method proposed by Boore and Abrahamson (2023) improves the results for V S 30  = 490 m/s in the shallow portion of the host profiles (L. Al Atik, personal communication, April 6, 2023). The updated profiles for V S 30  = 490 m/s (which represent a condition very close to V s 30 , sim  = 500 m/s) were provided by L. Al Atik and are used here to illustrate the application of Method 2.

Figure 6 presents the computation of the reference correction factor for the site PRPC, considering the CB14 model (results for the BSSA14 model are provided in Electronic Supplement C). As shown in Figure 6a, the different LF and HF simulation profiles result in slightly different merging depths at each frequency range. Given that the accuracy of the host profile estimation reduces with depth, this profile was modified to be equivalent to the simulation profile below the selected merging depth. The reference depth for the calculating RC was defined at 5000 m.

OPEN IN VIEWER

Figure 6.

Reference adjustment in Method 2 for the site PRPC: (a) LF and HF profiles, from 0 to 5000 m; (b) HF profiles from 0 to 220 m; and (c) resulting reference correction factors.

Figure 6c shows that the resulting reference correction factor is less than 1 for f  < 5 Hz, reducing the amplification (i.e. mitigating the overprediction issue) at relatively low frequencies in Method 1. Interestingly, Figure 6c reveals that the overamplification issue in Method 1 is not restricted to the LF component but also affects part of the HF range. As illustrated in Figures 6b and c, this is due to the difference between the HF simulation profile and the host profile for V S 30 , sim in the top ≈ 100 m. Although both profiles are characterized by V S 30 , sim , the HF simulation profile displays a constant value of V S  = 500 m/s down to a depth of 100 m, whereas the host V S increases with depth to a value close to 700 m/s, generating increasing levels of overamplification in Method 1 from f ≈ 5 Hz to f ≈ 1.2 Hz, as can been seen in Figure 6c. This overamplification is removed by Method 2, through the application of RC . For f > 5 Hz, RC is positive, which is due to the host V S values < 500 m/s at the very shallow depths. However, this has limited effect on S F 2 because a HF truncation is considered for the CB14 and BSSA14 models (developed in the response spectral domain) in the same fashion as in Method 1. The resulting S F 2 and its comparison with S F 1 are investigated in the “Comparison of the site amplification obtained with the five methods” section and in Electronic Supplement C.

Method 3—SRI-based amplification

Methods 1 and 2 only require proxy parameters (e.g. V S 30 , Z 1 . 0 ) to characterize the site conditions, which on many occasions is the only information available (e.g. in regional applications). If more site-specific site-characterization data exist, alternative approaches can be used. Method 3 uses a V S profile and the SRI method (Approach D) to represent the linear component of the SF, as follows:

S F 3 , L ( f ) = T F SRI , actual T F SRI , sim = ρ R V S , R ρ ¯ actual V ¯ S , actual ρ R V S , R ρ ¯ sim V ¯ S , sim · exp [ − π f κ 0 , actual ] exp [ − π f κ 0 , sim ]

(10)

where the subscripts actual , sim , and R indicate properties of the actual profile, simulation profile, and reference condition (i.e. elastic half-space considered for the amplification computation), respectively; and in particular, κ 0 , actual and κ 0 , sim are the HF site attenuation parameters estimated at the top of the actual and simulation profiles, respectively.

The objective of the denominator in Equation 10 is to remove the impedance- and attenuation-based site response introduced by the regional ground-motion simulation within the depth over which the site adjustment is considered. The use of the SRI method in this equation is fully consistent with the treatment of the linear amplification and HF attenuation adopted in the HF ground-motion simulation method (as described in the “Hybrid broadband ground-motion simulation methodology” section). Although the LF simulation is based on a different approach (3D time-domain wave propagation), the frequencies affected by T F SRI , sim are mainly in the HF range (i.e. f > 1 Hz) for shallow site adjustments of concern here. Nonetheless, as the transition frequency of hybrid ground motion simulations increases over time (e.g. f t > > 1 Hz), a more consistent approach for the LF range could be adopted.

Equation 10 can be rewritten as

S F 3 , L ( f ) = ρ ¯ sim V ¯ S , sim ρ ¯ actual V ¯ S , actual · exp [ − π f ( κ 0 , actual − κ 0 , sim ) ] = S F 3 , L , I ( f ) · S F 3 , L , D ( f )

(11)

where S F 3 , L , I and S F 3 , L , D are the amplification and attenuation operators for the S F 3 , L factor. S F 3 , L , I represents the relative linear (impedance-based) site amplification between the actual and simulation profiles with respect to a common elastic half-space, and S F 3 , L , D represents the relative site attenuation between the actual and simulation site conditions. Since κ 0 reflects the cumulative effect of damping through the full site profile (Afshari and Stewart, 2019; Campbell, 2009), the term ( κ 0 , actual − κ 0 , sim ) in Equation 11 accommodates any difference in the amount of damping represented by κ 0 , actual and κ 0 , sim , for a theoretical material column that can extend to a depth significantly greater than the near-surface (e.g. deeper than 100 m). In contrast, the amplification operator, S F 3 , L , I , is generally limited to capturing impedance effects in the near-surface (e.g. within the top 100 m), where site-specific V S data are most likely to be available.

Since Equation 11 (i.e. Approach D) is limited to the estimation of the linear site response, the nonlinear operator derived in Method 1 is proposed for capturing soil nonlinearity in Method 3. In this way, the total SF in Method 3 ( S F 3 ) is computed as

S F 3 ( f ) = S F 3 , L ( f ) · S F 1 , NL ( f )

(12)

Figure 7 illustrates the computation of the linear component of Method 3 for the site PRPC. Figure 7a shows the derivation of the actual V S profile and the selection of the elastic half-space. The definition of the half-space in this method requires that the actual and simulation profiles merge at a certain depth; however, the available site-specific V S profile for PRPC has a depth of 30 m and does not reach the simulation profile. To achieve this, alternative sources of information can be utilized to extend the site-specific profile to a proper depth. The V S profile extracted from the NZVM (without imposing a minimum V S ) is used, and the half-space is selected at a depth of 100 m. Given that for Method 3, it is assumed that the only site information available is the V S , actual profile, the actual soil density, ρ actual , is estimated using the V S -based correlation for unit weight provided in Rix et al. (2019).

OPEN IN VIEWER

Figure 7.

(a) V S profiles used in Methods 3 (and subsequently, in Method 4), for the site PRPC. (b) Resulting linear site factor.

Since the HF ground-motion simulation methodology uses the SRI method in the same fashion as Method 3, κ 0 , sim is a known input parameter. In particular, κ 0 , sim  = 0.045 s was used in all the regional simulations considered in this article. When this value is used in Equation 10, the generic HF spectral decay considered in the regional HF simulation is fully removed and replaced by a site-specific estimate based on κ 0 , actual . κ 0 , actual can be estimated using a number of different methods (e.g. Ktenidou et al., 2014), including its direct derivation from recordings and the use of correlations with V S 30 and other site parameters. In this study, κ 0 , actual is estimated using two different V S 30 -based correlations, proposed by Bayless and Abrahamson (2019) and Xu et al. (2020), which give values of κ 0 , actual ( BA 19 )  = 0.052 s and κ 0 , actual ( X 20 )  = 0.060 s, respectively, for PRPC. It must be acknowledged that κ 0 - V S 30 correlations generally present significant scatter (e.g. Ktenidou et al., 2014), with a standard deviation in natural log units in the order of 0.30 (e.g. Xu et al., 2020). Figure 7b presents the linear SFs from Method 3 (Equation 11) for the two estimates of κ 0 , actual , along with the linear site response associated with the actual ( T F SRI , actual ) and simulation ( T F SRI , sim ) profiles. The figure shows that the selection of κ 0 , actual has a greater impact at higher frequencies, especially for f > 4 Hz. In particular, the frequency at which S F 3 , L starts displaying deamplification vary significantly, with f ~ 13 Hz for κ 0 , actual ( X 20 ) and f ~ 30 Hz for κ 0 , actual ( BA 19 ) .

Figure 8 presents the SFs ( S F 3 ), along with the underlying operators, S F 3 , L , I , S F 3 , L , D , and S F 1 , NL , for the low-amplitude and high-amplitude motions, considering three alternative GMMs for the computation of S F 1 , NL . Figure 8b shows that changes in the slope of the linear amplification operator ( S F 3 , L , I ) are associated with the presence of additional velocity contrasts (shown in the inset within Figure 8d) in its computation, illustrating the direct relationship between depth and frequency in the SRI method. Particularly, the linear amplification at frequencies greater than 3 Hz is controlled by the top 12 m of the soil profile. Figures 8a and c show that in the case of the low-amplitude motion, the selection of the κ 0 model generally has a greater impact than the choice of the nonlinear model ( S F 1 , NL ). However, Figures 10b and d indicate that for the high-amplitude motion, the selection of the nonlinear model has a greater influence in the frequency range 0.5 < f < 2.5 Hz, and the sensitivity to the estimation of κ 0 , actual starts to dominate at f > 10 Hz.

OPEN IN VIEWER

Figure 8.

(a) Site-response operators of Method 3 for the site PRPC considering the low-amplitude motion, (b) the high-amplitude motion, (c) resulting site factors for the low-amplitude motion, and (d) high-amplitude motion.

Method 4—1D TF-based amplification

Method 4 is similar to Method 3 but uses the theoretical 1D TF (Approach C) to account for linear site effects. The linear component of the SF can be expressed as follows:

S F 4 , L ( f ) = | T F SH 1 D , actual ( V S , actual , ρ actual , D min , actual ) | ρ R V S , R ρ ¯ sim V ¯ S , sim · exp [ − π f Δ κ 0 , sim ]

(13)

where the numerator represents the linear site response for the actual profile, computed through the modulus of the outcrop 1D TF ( T F SH 1 D , actual ), and the denominator represents the linear site response for the simulation profile ( T F SRI , sim ), with respect to a common elastic half-space (more details regarding this are provided in the “Definition of the reference condition in Method 4” section of Electronic Supplement B). Since the objective of the denominator is to remove the site response introduced by the regional ground-motion simulation, the SRI method is used in the denominator of Equation 13 for consistency with the HF ground-motion simulation (similar to Method 3). In Method 4, this is implemented in the same fashion as the V S − κ 0 deconvolution approach proposed by Cabas and Rodriguez-Marek (2017). Δ κ 0 , sim is defined as Δ κ 0 , sim = κ 0 , sim − κ 0 , R , where κ 0 , sim and κ 0 , R are the HF site attenuation parameters estimated at the top of the simulation profile and at the top of the elastic half-space, respectively. The elastic half-space, and the actual and simulation V S profiles, are established in the same way as in Method 3 (see Figure 7a).

In the same manner as in Method 3, nonlinearity is modeled in Method 4 using the nonlinear operator from Method 1. Hence, the total SF in Method 4 ( S F 4 ) is computed as

S F 4 ( f ) = S F 4 , L ( f ) · S F 1 , NL ( f )

(14)

Other studies have also utilized theoretical 1D TFs to adjust hybrid broadband ground-motion simulations (e.g. Ojeda et al., 2021; Pilz et al., 2021), but either they have not accounted for soil nonlinearity (e.g. Ojeda et al., 2021) or have used a different approach to model it, such as the equivalent linear method (e.g. Pilz et al., 2021). The method proposed here to treat nonlinearity along with the linear amplification from the 1D TF is simpler and only requires V S 30 for the nonlinear component.

Several approaches have been proposed to estimate the small-strain damping D min (e.g. Afshari and Stewart, 2019; Cabas et al., 2017; Darendeli, 2001; Pretell et al., 2023; Xu et al., 2020), which can be used in Method 4 to compute D min , actual . Here, three formulations are investigated: (1) a laboratory-based damping, D m i n , a c t u a l l a b , computed using Darendeli (2001); (2) D min , actual lab with a multiplicative factor of 3 (Pretell et al., 2023); and (3) a V S -based damping, D min , actual V S , computed via D min = 1 / ( 2 Q ) (Kramer, 1996), where the quality factor, Q , is estimated using Q = 7 . 17 + 0 . 0276 V S (Campbell, 2009, Model 1), with V S in m/s. The multiplicative factor considered in (2) is based on the recognition that the intrinsic material damping measured in the laboratory ( D min lab ) does not fully capture the actual damping observed in the field due to additional mechanisms such as wave scattering (Afshari and Stewart, 2019; Cabas et al., 2017; de la Torre et al., 2022a), which results in the overprediction of the theoretical fundamental mode when D min lab is used (Pretell et al., 2023). The multiplicative factor should be ideally determined on a site-specific basis, which is possible, for example, in the case of borehole array sites (e.g. Afshari and Stewart, 2019; Xu et al., 2020). A generic value of 3 is used in this study based on Pretell et al. (2023), who found that this factor improves the overall site response prediction at multiple borehole arrays in Japan and California. The 100-m depth considered in this article to perform the site adjustments is within the depth range considered by Pretell et al. (2023) to derive this multiplicative factor (sites with depths mostly less than 230 m); caution should be used for significantly deeper profiles. The Campbell (2009) Q − V S relationship adopted in (3) was selected based on its previous utilization in several site response studies (e.g. Afshari and Stewart, 2019; Cabas et al., 2017; Xu et al., 2020; Zhu et al., 2022). Some of them (e.g. Cabas et al., 2017; Xu et al., 2020) have found that this relationship better captures the field attenuation inferred from downhole array sites than the Darendeli (2001) model ( D min lab ). Given that for Method 4, it is assumed that the only site information available is the V S , actual profile, it is considered that the plasticity index (PI) is 0 for all layers, and the V S -based correlation provided in Rix et al. (2019) is used for estimating ρ actual .

In the case of Method 3, the parameter κ 0 , sim was used in the denominator of Equation 10 to remove the HF attenuation introduced by the full simulation site profile. In Method 4, Δ κ 0 , sim is used in the denominator of Equation 13 instead, which only accounts for the portion of the site attenuation produced above the elastic half-space considered. Unlike κ 0 , sim in Method 3, Δ κ 0 , sim is not a direct input parameter of the HF simulation, and hence, Δ κ 0 , sim has to be estimated in Method 4. Under the assumption that this parameter represents additive attenuation effects over the simulation soil layers above the elastic half-space, Δ κ 0 , sim can be expressed as (Cabas and Rodriguez-Marek, 2017):

Δ κ 0 , sim = ∑ i = 1 n H sim , i Q sim , i V S , sim , i

(15)

where H sim , i , Q sim , i , and V sim , i are the thickness, quality factor, and shear-wave velocity of each layer i in the simulation profile, respectively, and the sum is performed over the n layers above the elastic half-space. Using Q sim = 1 / ( 2 D min , sim ) , Q sim can be estimated via any of the three formulations considered for D min , actual . For internal consistency, the same approach that is used for obtaining D min , actual is also used for estimating D min , sim . In all the cases investigated in this study, Δ κ 0 , sim is computed on the basis of a 100-m deep layer with V S , sim  = 500 m/s, which results in the following Δ κ 0 , sim values for the damping formulations based on D min lab , 3 · D min lab , and D min , actual V S , respectively: 0.00192 s, 0.00576 s, and 0.00953 s.

Figure 9 illustrates the computation of the linear component of Method 4 for the site PRPC. Figure 9a shows the D min , actual and D min , sim profiles, and Figure 9b presents the resulting linear SFs ( S F 4 , L ), along with the linear site response associated with the actual and simulation profiles. Figure 9a shows that the V S -based damping formulation ( D min , actual V S ) generally results in a larger damping than 3 · D min , actual lab for the site PRPC. Figure 9b illustrates that the selection of the damping formulation has a greater impact at the peaks of the actual linear site response, and at very high frequencies ( f > 10 Hz).

OPEN IN VIEWER

Figure 9.

(a) D min profiles used in Method 4 for the site PRPC. (b) Resulting linear site factor from Method 4 for the site PPRC.

Figure 10 shows the SFs ( S F 4 ), along with the underlying operators, S F 4 , L and S F 1 , NL , for the low-amplitude and high-amplitude motions, considering two alternative formulations for the computation of D min , and three alternative GMMs for the computation of S F 1 , NL . D min lab was excluded considering that the evidence suggests that 3 · D min lab or D min V S provide more realistic values (e.g. Cabas et al., 2017; Pretell et al., 2023; Xu et al., 2020). Figures 10a and c illustrate that in the case of the low-amplitude motion, the selection of the damping formulation generally introduce a greater variability than the choice of the nonlinear model ( S F 1 , NL ). On the other hand, Figures 10b and d show that for the high-amplitude motion, this sensitivity is frequency-dependent. In the frequency range 0.5 < f < 2 Hz, the selection of the nonlinear model is more relevant, whereas for f > 25 Hz, the estimation of damping plays a more significant role.

OPEN IN VIEWER

Figure 10.

(a) Site-response operators of Method 4 for the site PRPC considering the low-amplitude motion, (b) the high-amplitude motion, (c) resulting site factors for the low-amplitude motion, and (d) high-amplitude motion.

Method 5—1D time-domain wave propagation analysis

Method 5 is a two-step procedure that generally requires more site-characterization data than Methods 3 and 4. Similar to Method 4, it relies on the SH1D assumptions but uses 1D time-domain wave propagation analysis (Approach B) to capture the near-surface nonlinear site response. In this method, the actual profile is divided into two portions: the “actual 1” and “actual 2” profiles. As illustrated in Figure 11a, for the site PRPC, the actual 2 profile corresponds to the shallower portion, which is bounded by a significant impedance contrasts at the bottom (corresponding to the top of the Riccarton Gravel Formation, characterized by V S  = 400 m/s). As previously shown in Figure 2, the actual 2 profile is well characterized, with site-specific measures of V S and CPT. On the other hand, the actual 1 profile is the deeper portion of the actual profile, for which soil nonlinearity can be assumed negligible. In the case of the site PRPC, location-specific site-characterization data are not available for this depth range.

OPEN IN VIEWER

Figure 11.

(a) V S profiles used in Methods 5 for the site PRPC. (b) Resulting linear site factor.

Figure 12 illustrates Method 5. In Step 1, a linear SF, S F 5 , L , is applied to the simulated ground motion (separately for the LF and HF components, as in the previous frequency-domain methods). This linear SF is obtained using the same equation used in Method 4 to compute S F 4 , L (Equation 13), but with the “actual 1” properties in the numerator (instead of the “actual” properties). S F 5 , L is computed relative to the half-space 1, shown in Figure 11a for the site PRPC. Figure 11b presents the resulting linear SF, for two different damping formulations, illustrating relatively small variability due to damping selection. The output of Step 1 is an outcrop broadband (BB) acceleration time series, a input ( t ) , which is then used as input motion for the 1D time-domain nonlinear inelastic site-response analysis performed in Step 2. This analysis is conducted using the actual 2 profile and the half-space 2, illustrated in Figure 11a for the site PRPC.

OPEN IN VIEWER

Figure 12.

Illustration of Method 5. In Step 1, a frequency-domain adjustment is performed to the LF and HF simulated ground motions ( a sim ( t ) ), which results in a BB ground motion ( a input ( t ) ) used in Sept 2; in Step 2, a 1D time-domain nonlinear inelastic site-response analysis is performed, which results in the adjusted acceleration time series, a adj ( t ) , which accounts for shallow site effects.

The main reason for computing the site-response of the actual profile in two steps (i.e. using the actual 1 and 2 profiles) is to deal with different half-space 1 properties for the LF and HF components. In the case of the site PRPC, the difference between LF and HF V S at half-space 1 is minor, but other sites (e.g. CMHS) may display significant deviations. 1D Time-domain site-response analysis requires a broadband input ground motion, and Step 1 allows for obtaining it at a common reference condition for LF and HF (half-space 2). This procedure also has the benefit of reducing the size of the site-response model for time-domain nonlinear analysis, limiting it to the depth range where soil nonlinearity is expected and where more site-characterization data are available to estimate site-response model parameters. However, decoupling the site-response analysis via this two-step procedure may introduce some errors in the site-response computation (Cabas et al., 2022), and hence, the selection of half-space 2 must be done with caution. Cabas et al. (2022) showed, for example, that a half-space located at a shallow impedance contrast mitigates this issue.

For the 1D time domain site-response analysis, different codes (Hashash et al., 2020; Jsh9, 2023; McKenna, 2011) and constitutive models (Groholski et al., 2016; Shi and Asimaki, 2017; Yang et al., 2003) were used to investigate the impact of alternative modeling decisions on the resulting site adjustment. The results of these analyses, summarized in the “Method 5” section of Electronic Supplement C, show that when the comparison is performed in a consistent manner (i.e. using the same modulus reduction (MR) curves), the variability in the site adjustment is significantly reduced. Therefore, in the following section, the results for Method 5 are only provided for the case in which the program OpenSees (McKenna, 2011) is used, along with the constitutive models PressureDependMultiYield02 and PressureIndependMultiYield (Yang et al., 2003, 2008), and user-defined MR curves. The MR curves were defined based on the Darendeli (2001) model, with a large-strain adjustment (Yee et al., 2013) for consistency with the soil shear strength. The estimation of the soil properties (e.g. actual soil density, PI, friction angle, and undrained shear strength) was based on the V S profile, and the CPT, borehole, and SPT data available at each site, including the use of CPT-based correlations provided in Robertson and Cabal (2022). This resulted in the characterization of some soil properties (e.g. ρ actual ) with different values compared to those considered in Methods 3 and 4, where only the V S profile was used as site information.

Comparison of the site amplification obtained with the five methods

Prior figures for each method focused on intra-method differences due to input parameter and model alternatives, using the site PRPC as an example. This section provides a comparison between the five methods for all four case study sites. Figure 13 presents a comparison of the linear SFs ( S F L ) obtained with Methods 1–4. Given that the treatment of the nonlinear component of the SFs ( S F NL ) is the same for these four methods, only the linear component is compared in this figure, but the effect of soil nonlinearity in the resulting AFs is examined subsequently. Two alternative GMMs (CB14 and BA18) are considered for Method 1 ( V S 30 -based amplification), and the CB14 model is considered for Method 2 ( V S 30 -based amplification with host-to-target adjustment). The CB14 is considered because it has been used in several previous studies (e.g. Lee et al., 2022; Razafindrakoto et al., 2018), and the BA18 model is included to illustrate the differences between SA- and FAS-based models. Also, two different estimates of κ 0 , actual and D min , actual are examined for Method 3 (SRI-based amplification) and Method 4 (1D TF-based amplification), respectively. For each site, the predominant frequency of the near-surface soil profile, f site , is plotted, which is estimated based on the lowest frequency peak in S F 4 , L .

OPEN IN VIEWER

Figure 13.

Comparison of the linear factors obtained with Methods 1–4 for the sites: (a) PRPC, (b) CMHS, (c) RHSC, and (d) CACS. Dashed lines in Methods 1 and 2 indicate that the LF component is not used in the adjustment.

Figure 13 shows that in the four sites investigated, the linear SFs obtained with Methods 3 and 4 mainly affect the HF range and tend to a value of 1 at low frequencies, reflecting the fact that the adjustment is only accounting for shallow site effects (within the top 100 m); whereas Methods 1 and 2 generally produce considerable amplification in the LF range. Following this, and considering that the compatibility with the regional-scale simulation is explicitly controlled in Methods 3 and 4, in the case of Methods 1 and 2, the SF is only applied to the HF simulation component (such as in Lee et al., 2022). The implementation of this procedure is described in Electronic Supplement B. Figure 13 suggests that the frequency range of application of these SFs should be defined in a site-specific manner, but for Methods 1 and 2, it is assumed that a site-specific V S profile is not available.

In the case of the four sites considered, Method 1 produces higher amplification than Methods 3 and 4 in the vicinity of f  = 1 Hz, reflecting the previously discussed limitations of this approach. Method 2, which better accounts for the compatibility with the regional ground-motion simulation, reduces this relative overamplification, particularly in the frequency range 1 <  f  < 2 Hz, but generally does not eliminate it.

The ergodic assumption implicit in Method 1 is illustrated when comparing the linear SFs for the sites PRPC and CMHS. These two sites are characterized by very similar V S 30 values, resulting in almost identical S F L produced by Methods 1 with the CB14 model (which only depends on V S 30 ). Methods 3 and 4, on the other hand, produce significantly different linear site amplification at these two sites, reflecting the differences in the site-specific V S profiles (see Figure 2), which are accounted for by these methods.

Method 4 produces significant variations in the SFs over narrow frequency bands, whereas Method 3 results in SFs varying smoothly with frequency. In particular, Method 4 generates significantly greater amplification than Method 3 around f site in all four sites considered and, in some cases, also around higher frequencies, which is driven by resonance effects not modeled by Method 3.

Figure 13 illustrates that the definition of κ 0 , actual in Method 3 can have a considerable impact on the resulting linear SFs at high frequencies. This is particularly true in the case of the sites RHSC and CACS, which are characterized by relatively high V S 30 , actual values, resulting in lower estimates of κ 0 , actual . The V S 30 − κ 0 correlation of Bayless and Abrahamson (2019) in these cases results in a κ 0 , actual (i.e. κ 0 , actual ( BA 19 ) ) that is lower than the constant value κ 0 , sim  = 0.045 s used in the regional ground-motion simulations. This negative Δ κ 0 = ( κ 0 , actual − κ 0 , sim ) value (−0.001 s for RHSC and −0.007 s for CACS) produces HF amplification instead of attenuation, which explains the big differences observed at very high frequencies when comparing the S F L obtained using κ 0 , actual ( BA 19 ) and κ 0 , actual ( X 20 ) . The reason why the impact of κ 0 , actual is significant in Method 3, compared to the role that D min , actual has on Method 4 (as clearly observed in Figure 13d), is that the parameter κ 0 , actual accounts for cumulative effects of damping over the full site profile. In this way, Δ κ 0 can be accounting for differences produced by layers much deeper than 100 m, whereas D min , actual only acts in the shallow portion of the site profile. An alternative implementation of Method 3 that accounts for attenuation in a manner analogous to Method 4 is possible, but the current implementation was considered because of the convenience of using κ 0 , sim and κ 0 , actual as input parameters: κ 0 , sim is a direct input parameter (i.e. known value) of the regional ground-motion simulation, and there is a growing interest in the development of methods and models for estimating κ 0 at a given site (i.e. κ 0 , actual ) (e.g. Ktenidou et al., 2014; Xu et al., 2020).

Figures 14 and 15 present the AFs, which now include the influence of nonlinearity, that result from applying Methods 1–5 to the four sites investigated, for the low- and high-amplitude ground motion, respectively. For a given IM, AF is defined as

AF = I M adj I M sim

(16)

where I M adj is the IM for the adjusted ground motion, and I M sim is the IM for the simulated (unadjusted) ground motion. The IMs considered are SA for vibration periods between T  = 0.01–10 s, PGA, peak ground velocity (PGV), cumulative absolute velocity (CAV), Arias intensity (AI), and the 5%–75% and 5%–95% significant duration ( D s 575 and D s 595 , respectively). In the case of Methods 1–4, the same models and inputs used in Figure 13 are considered. For Methods 3 and 4, the nonlinear operator derived from the BA18 model is used, and for Method 5, only the OpenSees analysis using user-defined MR curves is examined. Figures 14 and 15 include the predominant period of the near-surface soil profile, estimated as T site = 1 / f site

OPEN IN VIEWER

Figure 14.

Amplification factors of different IMs resulting from applying Methods 1–5 to the sites: (a) PRPC, (b) CMHS, (c) RHSC, and (d) CACS, for the low-amplitude ground motion. For visual completeness, the AI AFs for the sites CMHS and RHSC associated with Methods 4 and 5 are plotted at the vertical axis limit (3.0). The actual values for CMHS are 3.45 (Method 4 using 3 · D min lab ), 3.25 (Method 4 using D min , actual V S ), and 3.49 (Method 5); and for RHSC, the actual values are 3.64 (Method 4 using 3 · D min lab ), 3.64 (Method 4 using D min , actual V S ), and 3.53 (Method 5).

OPEN IN VIEWER

Figure 15.

Amplification factors of different IMs resulting from applying Methods 1–5 to the sites: (a) PRPC, (b) CMHS, (c) RHSC, and (d) CACS, for the high-amplitude ground motion.

In the case of the low-amplitude motion (Figure 14), which produces relatively low levels of soil nonlinearity, the response spectral amplification generally follows a shape and amplitude similar to the linear SFs plotted in Figure 13 in the intermediate period range (approximately 0.1 < T < 1 s). However, at shorter and longer vibration periods, this is not the case. For example, there are several cases in Figure 13 where the linear SF produces de-amplification at high frequencies, but this de-amplification is not observed in AF (Figure 14) at short vibration periods, which is due to the difference between FAS and SA (Bora et al., 2016). At low frequencies, all the linear SFs tend to S F L  = 1 in Figure 13, whereas in Figure 14, AF > 1 values are observed at very long vibration periods. However, these periods are associated with relatively low SA values, and therefore, these AFs > 1 do not imply a noticeable increment in SA in absolute terms. These trends can be observed more clearly in the “Comparison between methods” section of Electronic Supplement C, where AF and SA are shown superimposed, for the two ground-motion intensity levels considered. Furthermore, Figure 14 illustrates that (1) the reduction in amplification produced by Method 2 relative to Method 1 in the SF domain (which mainly occurs in the intermediate frequency range) is also produced in the AF domain (for SA); and (2) the significant differences observed in the SFs at high frequencies between the CB14 and BA18 models in Method 1, and between the two estimates of κ 0 , actual in Method 3, are reduced considerably in the case of AF.

Figure 14 also shows that Method 4 and 5 generally display similar shapes and amplitudes in the response spectral AF for the low-amplitude motions, which is expected given that both methods rely on 1D site-response analysis. However, at some vibration periods, considerable differences are observed due to the different treatment of nonlinear site-response, the use of more site-specific data, and the decoupling of the site-response analysis in Method 5. For example, for the site CMHS, at very short vibration periods, Method 5 results in considerably lower values of amplification. This is mainly due to the more significant soil nonlinearity predicted by Method 5, which is modeled in a site-specific fashion, accounting for the relatively low values of cone resistance at this site (see Figure 2).

AI displays particularly high levels of amplification and method-to-method variability in Figure 14, especially in the case of the softer sites. On the other hand, the significant durations D s 575 and D s 595 are almost not influenced by the shallow site-response adjustments under the low-amplitude motions.

Figure 15 shows that for the high-amplitude motion considered, soil nonlinearity has a strong influence on AF, especially for the softer sites, PRPC, and CMHS. At these sites, severe de-amplification is observed at short-period SAs, PGA, and AI, and amplification is generally produced in the significant durations D s 575 and D s 595 . Figure 15 also illustrates the difference between an ergodic treatment of soil nonlinearity based on V S 30 (Methods 1–4) and a site-specific treatment as in Method 5. In the case of PRPC and CMHS, Method 5 generally resulted in stronger de-amplification at short vibration periods, and it was able to capture the softening of the soil deposits, manifested in the shift of the maximum response spectral AFs toward longer vibration periods, whereas the ergodic approach simply results in general reduction and smoothing of the linear amplification. This produces considerable differences in AF over certain period ranges between Methods 4 and 5 (e.g. between 1 and 2 s at CMHS). Also, Method 5 produces more significant amplification of D s 575 and D s 595 in these sites.