Small-strain site response of soft soils in the Sacramento-San Joaquin Delta region of California conditioned on V_S30 and mHVSR

The ground motion models (GMMs) used in California are from the Next Generation Attenuation (NGA) West2 project (Bozorgnia et al., 2014) and have site response components that apply either globally for active tectonic regions (e.g. Boore et al., 2014; hereafter BSSA14) or for broad regions (e.g. California vs Japan; Campbell and Bozorgnia, 2014). However, differences have been observed when investigating site effects at more local scales (e.g. Landwehr et al., 2016; Nweke et al., 2022). The application of ergodic site response models to regions with soils having unusual characteristics, such as the Sacramento-San Joaquin Delta region of California (hereafter Delta), carries large epistemic uncertainty. The peaty-organic soils in the Delta produce typical time-averaged shear wave velocities in the upper 30 m (V_S30) lower than 150 m/s (Figure 1c), which is softer than the lower limit for NGA-West2 site response models. Furthermore, the soft soils in the Delta typically overlie relatively firm, non-organic soils, which can give rise to more pronounced impedance and resonance effects than would be typical at non-Delta sites. These features of site response are not captured by common V_S30-based models, although they could be potentially captured with the addition of site response terms conditioned on peak frequency (f_p) as a site parameter (e.g. Harmon et al., 2019; Hassani and Atkinson, 2018a, 2018b; Kwak et al., 2017; Wang et al., 2022a).

Accordingly, the ergodic site terms in current GMMs used in California are expected to have bias and large uncertainty when applied to peaty-organic sites, such as those encountered in the Delta. In this article, we utilize a subset of an expanded ground motion database (relative to NGA-West2) to develop a local linear site response model based on features of site response derived from non-ergodic analyses. We condition our model on V_S30 as the primary site parameter, and present optional additive terms conditioned on parameters derived from microtremor-based horizontal-to-vertical spectral ratios (mHVSRs). The models presented here represent a refinement to the ergodic site response model in BSSA14 (i.e. Seyhan and Stewart, 2014; hereafter SS14). The results are expected to be useful for hazard and risk analyses in the Delta region. Moreover, the methodology developed for this study could be useful for future studies seeking to couple V_S30 with mHVSR-based parameters for regional- or local-scale site response modeling.

Following this introduction, we describe the ground motion data set and site parameters upon which the site response models are developed. We then present data analyses that include non-ergodic methods using residual calculations. The results are used to develop the V_S30-scaling model and the additive terms conditioned on mHVSR-based site parameters. We conclude by comparing the performance of the proposed models with SS14 and provide a summary and recommendations.

Since 1969, the Delta has been instrumented by 54 distinct seismic stations operated at different times by the Berkeley Digital Seismograph Network (network code: BK, number of stations: 2), California Department of Water Resources (WR, 11), California Strong Motion Instrumentation Program (CE, 15), and the United States Geological Survey (NC, 1; NP, 14; and YU, 11). This section presents the data compiled for analysis of non-ergodic site response at these stations.

We use a subset of the California ground motion database described in Buckreis et al. (2023a), which is a significantly expanded version of the NGA-West2 database (Ancheta et al., 2014) for California sites. We defer to that publication for details, but present here the selection criteria implemented as part of this work. We consider data from moment magnitude (M) ≥ 4.0 events originating in California or western Nevada that are well recorded with at least one record at any of the 54 distinct seismic stations located within the Delta. The limiting magnitude of 4.0 was used to minimize potential effects of spectral shape on data-derived site response, as discussed by Stafford et al. (2017). M-distance criteria defined by Boore et al. (2014) are enforced to minimize potential bias from poor sampling of distant sources, and we only use data within their usable bandwidth (PGA to a limiting oscillator frequency of 1.25 times the high-pass corner frequency, f_c,hp, used during signal processing). Horizontal components are combined to median-component (RotD50) intensity measures (IMs) as defined by Boore (2010), and the lowest longest usable period of the two individual components is used. Only data from sites which have recorded at least four earthquakes are used during model development. Figure 1 shows the locations of 69 events, 36 “Delta” stations, and 45 “non-Delta” stations that satisfy our selection criteria and were subsequently used in this study.

As shown in Figure 2, most of the data are associated with small magnitude events (M < 4.75) originating relatively far from the Delta (Joyner-Boore distance, R_JB > 100 km). These motions are relatively weak as shown by the distribution of median-component peak ground accelerations (PGA) in Figure 2b. These observations provide evidence that nonlinearity is not expected to be significant in the data set; therefore, linear site response is anticipated, which is useful for developing the linear component of a local site amplification model.

The Delta is the largest freshwater tidal estuary on the west coast of the United States, which provides geologic conditions well suited for the deposition of peaty-organic soils (Drexler, 2011). Peat makes up about 18% of soils encountered within depths of 0–10 m across much of the Delta with notable exceptions on Sherman, Twitchell, and several other neighboring islands where deposits can be up to 15 m thick (Department of Water Resources (DWR), 2013, 2015). Underlying the peat are interbedded layers of clays, silts, and liquefiable sands. It is important to point out that not all Delta sites contain peat; therefore, the site effects associated with the peaty-deposits will not be observed across all Delta sites. This subsection presents the geotechnical data used to assign site parameters for the purpose of modeling local site response within the Delta. We have not considered sediment depth effects because the available basin depth model (USGS SFCVMv21.1; Hirakawa and Aagaard, 2021) provides uniform depth estimates (e.g. depth to 1000 m/s shear wave iso-surface, z_1.0) across the study region.

The most commonly used site parameter in site response models is V_S30, which is computed as the ratio of 30 m to the shear wave travel time in the upper 30 m of the site (e.g. Borcherdt, 1994). Typical V_S30 values in the Delta computed from measured V_S profiles uploaded to the community shear wave velocity database (www.vspdb.org; Kwak et al., 2021) range from approximately 60 to 350 m/s (Figure 1c), where stations with V_S30 > 275 m/s are unlikely to have peat. V_S30 values at seismic stations lacking measured V_S data are assigned using a locally calibrated peat-thickness proxy-based model (Buckreis, 2022) where peat-thickness is taken from nearby borings or cone penetration tests (CPTs), or is estimated using the approach of Deverel and Leighton (2010) in which kriging interpolation of peat thicknesses from over 1100 borings is applied. The range of V_S30 values for Delta stations used in this study is from 100 to 390 m/s, as shown in Figure 1(c) and summarized in Table 1.

Peaks in amplitude–frequency mHVSR plots can identify resonance features in the site response (Bonilla et al., 1997, 2002; Cadet et al., 2012; Field and Jacob, 1993, 1995; Nakamura, 1989; Nogoshi and Igarashi, 1970, 1971; Satoh et al., 2001; Theodulidis et al., 1996). These peaks and other mHVSR features are evaluated in this study for their potential to improve site response predictions and reduce uncertainties relative to those from ergodic V_S30-based models. We derived mHVSR at 41 stations using data from temporarily deployed seismometers (Buckreis et al., 2021), and ambient noise data for 11 permanently deployed seismometers queried from the Incorporated Research Institutions for Seismology (IRIS; Trabant et al., 2012) using the “ObsPy” package in Python (Krischer et al., 2015). The mHVSR were processed using procedures described in Wang et al. (2022b), which are revised from earlier protocols (Site EffectS assessment using AMbient Excitations, 2004).

Each mHVSR curve is independently assessed for peak features using the automated algorithm presented by Wang et al. (2023). The peak detection algorithm implements a regression tree to represent the mean mHVSR curve as a step-wise function, from which relative amplitudes and lengths of adjacent steps are used to identify peak features, as illustrated in Figure 3. Sites identified to possess peak features are fit using a Gaussian function (Ghofrani and Atkinson, 2014):

where f_p is the fitted peak frequency for the ith HVSR peak (i.e. the location of the HVSR peak, regardless of whether it actually corresponds to a site resonant frequency), c₁ is peak amplitude, w_p is peak width, c₀ is a frequency-independent constant, and f is frequency in Hz. The absolute amplitude of the peak (denoted a_p) is calculated as the sum of c₀ and c₁. The fit operation estimates f_p, c₀, c₁, and w_p for peak i at each site with a significant peak.

Note that the lowest peak frequency (f_p1) is not necessarily equal to the fundamental site frequency (f₀) for various reasons, including (1) f₀ could be out of the mHVSR usable frequency range for sites in deep basins, and (2) horizontal amplitudes near f₀ may not exceed the instrument noise threshold. Table 1 presents a summary of the peak identification and fitting results (lowest frequency peak only) for the Delta stations used in model development.

Approximately 71% of all Delta sites are identified as possessing peaks, which is significantly more than what has been observed across California as a whole (e.g. 27%–52%; Wang et al., 2023). One factor which produces peak features is the presence of a strong impedance contrast in the soil column, such as that observed between near-surface soft soils (including peaty layers and soft clays) and stiffer underlying non-organic soils. This strong impedance contrast is expected to give rise to significant site response effects.

Non-ergodic analyses were performed to characterize local site response in the Delta region. This section presents the approach that was applied to support model development.

Total residuals (R_ij) represent differences between data and the median predictions of a GMM:

where $ln(Yij)$ is the natural logarithm of the observed IM at site j from event i, and $μij$ is the natural log median from a GMM conditioned on the indicated parameters. The BSSA14 GMM with subregional anelastic path adjustments as described in Buckreis et al. (2023b) is used here. Total residual R_ij can be partitioned using mixed-effects analyses (Abrahamson and Youngs, 1992; Gelman and Hill, 2006; Sahakian et al., 2018):

where $ck$ represents the model bias for GMM k; $ηE,ik$ and $ηS,jk$ represent random effect terms that quantify the event- and site-specific biases, respectively; and $εijk$ is the remaining residual. For brevity, the subscript k is omitted from subsequent equations and variables. This partitioning is performed on the complete California data set (Buckreis et al., 2023a) to optimize reliability of the random effects and their associated standard error terms.

Site response models (F_S), whether site-specific or ergodic, usually can be expressed in the following general form (e.g. Stewart et al., 2017):

where the first and second terms represent the linear (⁠$f1$⁠) and nonlinear amplification, respectively; $xIMr$ is the amplitude of shaking for a reference site condition (generally rock) for a particular event at a particular site (expressed as an IM, often PGA); $f2$ represents the slope in ln(amplification)-ln(⁠$xIMr$⁠) space for $xIMr>>f3$⁠; and $f3$ represents the transitional value of the reference site IM below which the site response is nearly linear, and above which the trend of amplification with $xIMr$ becomes nearly linear in log-log space (e.g. 0.1 g in SS14).

Non-ergodic site response methods take the site response as the sum of the site response model (used in the GMM applied for residual analyses) and the site term. For this work, the linear site response (⁠$f1$⁠) was estimated as the sum of the site term (established from weak motion data) and the linear portion of the site response GMM:

where $(f1)jo$ represents the data-derived estimate of the linear site response (the “o” superscript is for “observed” and should not be confused with the model estimate of $f1$ from Equation 4) relative to the GMM reference condition (V_S30 = 760 m/s). The subtraction of the nonlinear term in Equation 5 is to remove first-order nonlinear effects in the site response estimate, although this term is nearly zero in the present case.

The resulting site response quantities are plotted in Figure 4 for 24 of the 36 Delta sites, each of which has reliable $ηS,j$ estimates (at least four usable records per period, T) for all periods with T ≤ 9 s, and are arranged by increasing V_S30. All sites exhibit the same general trends: saturation (i.e., flat trend) at short periods, slight decreases at around 0.1 s, and then increasing amplification that sometimes saturates at longer periods. In addition to these general trends, many sites exhibit peak features at a specific oscillator period (e.g. WR_SHER, NP_LVA3, NP_LVA4, YU_SRT, CE_67587, and many more). These features are likely associated with resonance of the upper soft soils (with very low V_S) relative to deeper, stiffer sediments (with higher V_S). The locations of the peaks span over a range of periods (0.2–2 s), which may relate to the natural periods of each site.

The modeling approach outlined in this section aims to develop a linear site response model for the Delta region using weak ground motions (relatively low amplitudes), where significant nonlinear effects are not expected. The computed $(f1)jo$ values from the prior section represent the average site amplification observed over multiple events for station j and are taken as the observations from which site amplification models are developed.

The advantage of using $(f1)jo$ over within-event residuals (⁠$ηS,j+εij$⁠) during model development is that $(f1)jo$ reduces the influence of the number of observations at a given site in regression methods. This is important because sites with the most recordings would otherwise control the regional site response and sites with few recordings would have limited influence. Moreover, $(f1)jo$ provides a direct quantification of site response, whereas within-event residuals quantify the error that remains after accounting for source effects, which may include errors in the site response or path effects. However, the disadvantage is that it can be difficult to accurately estimate $(f1)jo$ at sites with relatively few recordings—even through mixed-effects methods. To address this issue, we use only data which have recorded at least four events (per IM) and consider their associated standard errors (SE_j) during model development through weighted least squares (WLS) regression (Strutz, 2016). Each observation is assigned a weight (w_j) inversely proportional to the square of its SE:

where n represents the number of sites used during regression for each IM. It follows that stations with smaller SE_j values will have greater influence during the regression.

The proposed model was developed in three stages, as described in the following subsections. The first develops the V_S30-scaling model, and the second and third develop additive terms conditioned on mHVSR peak parameters and average mHVSR amplitudes, respectively.

Ergodic site response models are typically conditioned on V_S30 to capture first-order site effects (e.g. SS14; Parker and Stewart, 2022). Data from 36 Delta and 45 surrounding non-Delta stations (Figure 1) are used to develop the V_S30-scaling model. The motivation for including the 45 non-Delta stations at this stage of model development is as follows:

1. The Delta is located within the Central Valley—a large sedimentary basin structure, where we anticipate basin effects that are more readily investigated using a larger population of sites, and

2. A V_S30-scaling model developed using solely data from Delta stations would only be applicable to the parametric range of the data (i.e. 100 < V_S30 < 300 m/s). The resulting model would likely be biased for stiffer soil sites near the outer-edges of the Delta, so data in the general vicinity with similar geologic conditions (i.e. within the Central Valley) are used to constrain the model for moderately-soft-to-stiff site conditions.

In Figure 5, the $(f1)jo$ terms for these stations are plotted against V_S30, and the ergodic SS14 linear V_S30-scaling model is shown with a black line for reference for six IMs. Binned means are included along with their 95% confidence intervals to highlight data trends. In general, the ergodic model reasonably captures the observed amplification and V_S30-scaling (slope) of non-Delta stations for peak ground velocity (PGV), PGA, and pseudospectral acceleration (PSA) at periods shorter than approximately 3.0 s. However, SS14 under-predicts the amplification and V_S30-scaling for non-Delta stations at longer periods. This observation provides evidence that the ergodic model is generally appropriate for typical non-peaty sites in the greater-Delta area, while the misfits at long periods may result from sedimentary basin effects in the Central Valley.

For Delta stations, the ergodic model tends to over-predict amplification for PGA and periods shorter than about 1.0 s. Across all IMs, the V_S30-scaling appears to saturate (i.e. slope goes to zero) for exceptionally soft soil sites (V_S30 < 150 m/s), which is consistent with the findings of recent studies for other regions (Parker et al., 2019 for central and eastern North America; Parker and Stewart, 2022 for global subduction zones; Nweke et al., 2022 for southern California basins). As with non-Delta stations, the ergodic model appreciably under-predicts the site response at periods longer than 3.0 s.

Based on these observations, we modify SS14 as follows:

where c₁, c₂, and c are V_S30-scaling coefficients; V₁, V₂, and V_c are limiting velocities; and V_ref is the reference site condition where F_lin = 0 (V_S30 = 760 m/s). The first two model components (first and second rows) effectively represent the V_S30-scaling in the Delta for exceptionally soft and soft-to-moderately-stiff site conditions, respectively. The third and fourth components match the SS14 model; the only modification is an added lower limiting velocity on the third component (i.e. V₂). The coefficients c and V_c are adopted from SS14, and c₁, c₂, V₁, and V₂ are estimated for the Delta region.

The model is first fit to the data with all four free parameters using WLS regression. V_S30-scaling coefficients c₁ and c₂ are constrained to be non-positive, since site response has been observed to increase as sediment stiffness decreases across many investigations spanning over 50 years (e.g. Borcherdt and Gibbs, 1976; Borcherdt and Glassmoyer, 1994; Idriss, 1990; Seed et al., 1976). V₁ is constrained to be less than V₂, which is constrained to be less than or equal to V_ref to ensure a continuous function and to enforce a reference velocity of V_ref = 760 m/s. From these results, c₁ was identified to be the most stable parameter in terms of lacking sudden between-period fluctuations and was set to a constant value of zero for all IMs except for PGV, effectively removing the first term in the top model component in Equation 7.

A second round of WLS regression was performed with c₂, V₁, and V₂ estimated as free parameters. These results indicated that V₁ was the most stable parameter. Regressed values of V₁ were smoothed using LOESS regression, which is a type of locally weighted scatterplot smoothing (Cleveland, 1979). A third and fourth round of WLS regressions and smoothing were conducted to set V₂ and c₂, respectively. Figure 6 presents the as-regressed and final smoothed-coefficients of the proposed V_S30-scaling model, the results of which for individual IMs are represented by the lines in Figure 5. Coefficients are provided in Table S1 in the supplement.

The proposed Delta-calibrated V_S30-scaling model is implemented as follows:

where the site response is solely a function of V_S30 and PGA at the reference rock condition (V_S30 = 760 m/s; PGA_r) (i.e. F_nl). The “v” superscript indicates that the median prediction is computed using a model that only includes V_S30-scaling.

The variability of $(f1)jo$ data with respect to V_S30 appears from visual inspection of Figure 5 to be quite large (especially for short-to-intermediate periods). This suggests that the site response in the Delta may be influenced by phenomena that are poorly represented by V_S30; subsequently, we consider additional site parameters for their ability to capture additional site response features.

Site resonance effects cannot be predicted by simple V_S30-scaling models, such as that described by Equation 7, because of the necessary smoothing across site-specific features in model development. In this section, we investigate whether resonances observed in site response are related to peaks in mHVSR. There are two general types of HVSR measurements—from microtremors (denoted mHVSR) and earthquake recordings (eHVSR). Most past research relating HVSR to site response has used eHVSR, including Cadet et al. (2012), Zhao and Xu (2013), Ghofrani et al. (2013), and Kwak et al. (2017) for Japan; Hassani and Atkinson (2016, 2018b) for central and eastern North America; Hassani and Atkinson (2018a) for California; and Panzera et al. (2021) for Switzerland. The use of eHVSR is convenient because the data required to measure HVSR are readily available from earthquake ground motion databases at ground motion recording sites; however, the problem is that the recordings used to develop the HVSR site parameters are not independent of those used to compute the site response. As a result, the correlation between eHVSR site parameters and site response is likely to be over-estimated. In addition, eHVSR is not likely to be available for forward applications at most sites. For these reasons, we only use mHVSR as the basis for site parameters used to develop the model. A similar approach was used previously for another local region having peaty-organic soils—Obihiro in Hokkaido, Japan (Wang et al., 2022a).

The Delta is rather unique in California because the geologic environment and depositional history results in a high percentage of sites (∼71%) that have significant mHVSR peaks, which is suggestive of strong impedance contrasts. It follows that site parameters, such as those summarized in Table 1, may have potential to improve site response predictions in the Delta. This subsection addresses the following questions: (1) can mHVSR be used to predict the presence of peaks in the site response, and (2) when mHVSR contain peaks, can the associated site parameters be used to improve predictions of site response? Of the 36 Delta sites with usable data, 34 have measurements of mHVSR and can be used to answer these questions.

To objectively assess whether resonance manifests at a given site, we developed an algorithm similar to that proposed by Wang et al. (2023) for HVSR to identify peak features in site response. In general, “peaks” are defined as features possessing the following key attributes: (1) relatively localized (i.e. the width should not span too large of a period range), (2) have sufficiently large mean amplitude relative to adjacent periods, and (3) have sufficient confidence that the feature is meaningful (i.e. uncertainty in amplitudes or periods should not be too large).

Total residuals (⁠$Rijv$⁠) are computed using the predictions given by Equation 8 and partitioned using mixed-effects methods to estimate site terms (⁠$ηS,jv$⁠), which represent the difference between site-specific observed amplification and the GMM described by Equation 8. These site terms (⁠$ηS,jv$⁠) can be considered to represent the unmodeled site response after V_S30-scaling effects have been accounted for, which is referred to below as relative site response.

As in Wang et al. (2023), a regression tree (Breiman et al., 1984) is used to effectively smooth and simplify the empirical relative site response as a piecewise function of non-overlapping linear segments (i.e. steps), to facilitate peak identification. The peak detection algorithm operates on the stepped results of the regression tree to identify peaks based on the three criteria previously defined. These conditions are assessed using the algorithm presented in Buckreis (2022). A Python code implementation, which includes documented examples, is available at https://github.com/tristanbuckreis/srPeak.

Figure 7 presents plots of $ηS,jv$ versus period and includes the results of the automated peak detection algorithm (“P” is indicated when a peak is identified, “NP” is indicated when no-peak is identified). Of the 36 sites with reliable $ηS,jv$ estimates, 17 (47%) are identified as possessing a resonant peak feature. Included in Figure 7 are plots of mHVSR, in which the typical frequency axis is transformed to period (as the inverse). The linkage between mHVSR peaks and relative site response peaks is addressed in the next subsection.

Figure 7 shows plots of mHVSR curves and fitted Gaussian functions from Equation 1, applied when peaks were identified, for all sites except NP_LVB4 where mHVSR data were not available. A comparison to peaks from site response (previous section) indicates coincident peak features for about 71% of cases. This suggests that the existence of an mHVSR peak is a potentially useful indicator of when site resonances can be anticipated, although the imperfect correlation suggests other (un-parameterized) factors also affect site resonances.

Investigations of correlation between relative site response and mHVSR peaks suggests that sites with clear relative site response peaks generally coincide with mHVSR peaks having tall relative peak amplitudes (a_p ≥ 3) and small tail amplitudes (c₀). This is shown in Figure 8a where sites with and without relative site response peaks are plotted in c₀-a_p space. Using logistic regression (Lottes et al., 1996), we propose a probabilistic model conditioned on c₀ and a_p to predict the presence of relative site response peaks:

where P_p is the probability of a resonant peak (0 ≤ P_p ≤ 1), $β0$ = 19.2471 ± 11.2033, $β1$ = 3.8467 ± 3.0136, and $β2$ = 4.3943 ± 2.3088. Figure 8a illustrates the variation of P_p in c₀ versus a_p space. A relatively narrow band separates sites predicted to have or not have site resonance peaks. The only outlier is site YU_HOL1 (c₀ = 1.150, a_p = 3.766); as shown in Figure 8b, this site was not identified as having a $ηS,jv$ peak but has predicted P_p = 0.85. Different analysts may have different opinions on whether or not YU_HOL1 has a $ηS,jv$ peak; however, the algorithm did not identify a peak because the relative uncertainty is quite large (the site recorded only four events).

The logistic model described by Equation 9 requires mHVSR peak attributes. If the mHVSR is not found to possess a peak, P_p is taken as zero. If we impose a threshold probability of 50% (i.e. if P_p < 0.50 there is no $ηS,jv$ peak, and if P_p > 0.50 there is a $ηS,jv$ peak), the predictive model correctly predicts $ηS,jv$ peaks from mHVSR approximately 91% of the time, which is considered to be a sufficient indicator for application purposes.

For modeling purposes, we parameterize relative site response (⁠$ηS,jv$⁠) peaks for a given site j with a pulse function described as

where $fp1,j$ represents the peak frequency of the pulse, $α1,j$ is the amplitude of the pulse, $ω1,j$ is the width of the pulse, $βj$ represents the amplification level at short periods, and $δ1,j$ represents the amplification level at long periods. The first term represents a Ricker wavelet (Ryan, 1994) that captures the observed response at periods shorter than the peak resonance including a dip centered at about 40% of $fp1,j−1$⁠, and the second term is provided to capture relatively flat responses at longer periods. Term $η¯S,jv(fp1,j−1)$ represents the model prediction at the peak period (⁠$fp1,j−1)$⁠. Example fits of this pulse function are given in Figure 9. Note at this stage $βj$ is taken as a constant (period-independent), despite variations of $ηS,jv$ with period (outside of the width of the pulse) that may be present in some cases.

Equation 10 estimates the properties of site amplification peaks where relative site responses have been derived from observations (i.e. interpreted ground motions). To generalize the model for sites without recorded ground motions (but with mHVSR), we relate fit coefficients $fp1,j$⁠, $α1,j$⁠, $βj$⁠, $ω1,j$⁠, and $δ1,j$ to mHVSR peak site parameters (f_p, c₀, c₁, a_p, and w_p), and examine possible correlations to develop a predictive model. As illustrated in Figure 10, mHVSR peak frequency (f_p) exhibits the strongest correlation with the fit parameters. The proposed functional form of the predictive model for relative site response (denoted f_1,p) is conditioned on f_p from mHVSR, and is consistent with Equation 10 with the exception of two changes: (1) the $βj$ and $δ1,j$ coefficients are dropped because they show no apparent correlation to any predictor variable, and (2) the amplitude is scaled by P_p:

where the f_1,p term in the second row of Eq. 11 is its value at T=f_p^-1 as evaluated using the first row, and

and

where f_p is in Hz; coefficients $a0$⁠, $a2$⁠, and $a5$ control the scaling of their respective terms with $fp$⁠; $a1$⁠, $a3$⁠, $a4$⁠, $a6$⁠, and $a7$ are constants; and $fL1$–$fL2$ are limiting frequencies that represent the transition from scaled-behavior (i.e. the fitted parameter scales linearly with $fp$⁠) to constant-behavior. The values for coefficients $a0$–$a7$ and $fL1$–$fL2$ are given in Table 2.

The updated GMM, which uses the V_S30-scaling model described by Equation 7 and the mHVSR peak–based model in Equations 11–14, is described by

where the linear site response is a function of V_S30, f_p, c₀, and a_p. The “v, p” superscript indicates that the median prediction is computed using a model that includes V_S30-scaling and peak-resonance effects.

Some previous research efforts have proposed to use broadband mHVSR features (not just peak attributes) for the prediction of site response. Senna et al. (2008) proposed a site response model for Japan conditioned on mHVSR and geologic or topographic units, where the shape of a reference site spectrum is modified by a period-dependent factor based on the mHVSR shape. Pinilla-Ramos et al. (2022) propose an mHVSR-based site response model for California, which similarly uses the whole period-dependent spectrum to predict amplification. These approaches differ from the work presented above because they use mHVSR ordinates directly to predict site response instead of using mHVSR primarily as a means by which to identify site resonances and their effects on amplification. Pinilla-Ramos et al. (2022) present a model conditioned on mHVSR amplitudes and V_S30, where mHVSR amplitudes were found to have the greatest influence for periods between 0.5 and 4.0 s. Their findings suggest that mHVSR amplitudes may have value for site response prediction, which is examined in this subsection.

Within-event residuals (⁠$δWijv,p$⁠) are computed using the data and predictions given by Equation 15 and are partitioned using mixed-effects methods to estimate site terms (⁠$ηS,jv,p$⁠). Those site terms represent the difference between site-specific observed amplification and the GMM described by Equation 15. The site response represented by $ηS,jv,p$ is referred to as residual site response, which captures remaining features after V_S30-scaling and peak-resonance effects have been accounted for.

General trends of $ηS,jv,p$ versus period are relatively flat for short and long periods, but those two period ranges have different amplitudes. An amplitude transition occurs across a range of intermediate periods, as illustrated in Figure 11. That behavior is fit using a logistic function as follows:

where $α0$ is a constant which quantifies the average long period amplification in natural log units, $Δα$ represents the difference between average amplification at short and long periods, and T_c is the period on which the transition in amplification is centered. The width of the smooth transition is fixed as 1.923 and was selected based on trial and error by visual inspections of all sites.

For modeling purposes, we relate the arithmetic average of $ηS,jv,p$ (⁠$μη$⁠) across short- and long-period bands to average mHVSR amplitudes (⁠$mHVSR¯$⁠) across different frequency bands. For example, we compare $μη$ computed over the period range between 0.01 and 3 s and $mHVSR¯$ computed over frequency range between 0.33 and 50 Hz, as shown in Figure 12a, which suggests a positive correlation. We repeated this exercise for many different combinations of $μη$ and $mHVSR¯$⁠, and consistently observed meaningful correlation for short period $μη$ with $mHVSR¯$⁠, but did not observe meaningful correlations for long period $μη$ (e.g., between 1 and 10 s) with $mHVSR¯$⁠.

Based on these observations, we propose a model (denoted f_1,m) conditioned on $mHVSR¯$ to capture this behavior:

where $α2$ represents the average amplification at short periods. The model adjusts residual site response over a short-to-intermediate period range which smoothly transitions to zero at long periods, as shown in Figure 12b. Coefficient $α2$ depends on $mHVSR¯$ as

where $a8$ = –0.2696 ± 0.0770, $a9$ = 0.5760 ± 0.2122, $a10$ = 0.7435 ± 0.2516, $a11$ =0.1334 ± 0.0715, $AL1$ = 0.8228, and $AL2$ = 1.5224. The coefficients in Equation 18 correspond to the optimal bandwidths that produced the highest correlation when relating $η¯S,jv,p$ and $mHVSR¯$⁠, which were found to be 0.01–3 s and 0.33–50 Hz, respectively, as shown in Figure 12a.

The final proposed GMM, which implements the local V_S30-scaling model and is informed by mHVSR, can be summarized as

where the linear site response is a function of V_S30, f_p, c₀, a_p, and $mHVSR¯$⁠. The “v,H/V” superscript indicates that the median prediction is computed using a model that includes V_S30-scaling and is informed by HVSR data. The nonlinear model will be additive to Equation 19 and remains under development.

The proposed Delta-specific local model begins with the BSSA14 GMM, updates the anelastic path model as presented in Buckreis et al. (2023b), and updates the linear site response using one of two potential locally calibrated models. The first model requires V_S30 only and is denoted F_lin(V_S30) (Equation 7). The second model uses V_S30 in combination with mHVSR parameters and is denoted F_lin (V_S30, f_p, a_p, c₀, $mHVSR¯$⁠) (Equations 7, 11, and 17). This section presents a comparison of the predictions of these two models and the global ergodic model presented by SS14. Figure 13 presents individual plots for Delta sites showing linear amplification versus oscillator period as derived from non-ergodic analyses and as predicted using each site response model.

The results shown in Figure 13 illustrate the significant misfit when extrapolating SS14 to the soft-soil conditions encountered in the Delta, and the improved fits that are achieved when using the locally calibrated models. The inclusion of the mHVSR terms is observed to better capture site-specific effects (e.g. resonant peaks and modest short-period amplitude adjustments) as compared to the V_S30-scaling-only model. The following subsection quantitatively compares model bias and variability, and present aleatory variability models.

Model bias is evaluated as the mean misfit of site terms (⁠$ηS,j$⁠) computed when implementing each of the site response models. Figure 14 presents plots of average $ηS,j$ versus oscillator period for all Delta sites, sites likely without resonant peaks (P_p < 0.5), and sites likely with resonant peaks (P_p > 0.5). As expected, the locally calibrated models perform better than SS14 (i.e. are less biased). The substantial bias of SS14 at short-to-intermediate and long periods demonstrates the need for local site factors for these exceptionally soft soils.

When examining biases of the two local models for all sites (Figure 14a), there are negligible differences (biases for each are effectively zero). However, when examining P_p < 0.5 and P_p > 0.5 sites separately (Figure 14b and 14c), the model that considers mHVSR (F_lin(V_S30, f_p, a_p, c₀, $mHVSR¯$⁠)) generally performs better than F_lin(V_S30), except when T > 3.0 s, where the results are identical. The F_lin(V_S30) model under-predicts on average ground motions at short-to-intermediate periods (T < 3.0 s) for non-peak sites, and over-predicts ground motions for this same period range for peak sites. When averaged across all sites, the model is unbiased. The F_lin(V_S30, f_p, a_p, c₀, $mHVSR¯$⁠) model reduces these biases for both site types.

Aleatory variability (or dispersion) terms are computed from residuals calculated for each of the three GMMs (original BSSA14 with subregional path and SS14 for site; previous GMM but with Delta-specific V_S30-conditioned site response; previous GMM with both Delta-specific V_S30-scaling and mHVSR-conditioned amplitude adjustment). Dispersion terms computed include total (σ), between-event (τ), and within-event (ϕ) standard deviations, which are related as

The within-event term can be partitioned into site-to-site (⁠$ϕS2S$⁠) and single-station (⁠$ϕSS$⁠) standard deviations as follows:

Figure 15 shows each of the dispersion terms for the California data set and the Delta subset. The local site response models for the Delta region do not affect the results of the California data set, which is expected since the Delta subset represents only a small fraction. However, there are significant changes of within-event variability and its components (⁠$ϕS2S$ and $ϕSS$⁠) for the Delta subset for most oscillator periods. Delta $ϕSS$ values are sometimes larger than the values computed for the California data set; however, their 95% confidence intervals generally capture the ergodic global $ϕSS$ model values proposed by Goulet et al. (2021). Accordingly, we do not consider a Delta-specific $ϕSS$ model to be justified, which is consistent with the conclusions reached by Buckreis et al. (2023b) for the California data set.

Site-to-site variability (⁠$ϕS2S$⁠), which represents the uncertainty in site response modeling, is null when site response is non-ergodic, but is non-zero otherwise (such as for the models presented in this article). A reference model for $ϕS2S$ based on residual analyses of the NGA-West2 data set was developed by Al Atik (2015) and refined by Goulet et al. (2018; Gea18). Such models are available for all NGA-West2 GMMs, including BSSA14. In those models $ϕS2S$ was found to be M-dependent, as shown in Figure 16. We compare computed mean values of $ϕS2S$ and their 95% confidence intervals to the NGA-West2 results reported by Gea18. The local $ϕS2S$ for the Delta are significantly lower than those in the global model (even when using SS14), because the variations of site response within the local Delta region are much smaller than what is observed at a global scale.

The application of F_lin(V_S30) does not significantly change variability when compared to SS14 for T < 1.5 s; however, significant reductions are observed at longer periods. V_S30 is the only dependent variable used in both models; however, the multilinear form of F_lin(V_S30) is better able to capture the trends observed in site response at long periods. The mHVSR-informed model reduces variability with respect to SS14 and F_lin(V_S30) for T < 3.0 s. For longer periods, the reduction of $ϕS2S$ relative to SS14 is identical for the F_lin(V_S30) and mHVSR-informed models, which likely occurs because there are no site resonant periods in this range. These observations are consistent for small and large M (at different levels) and provide evidence warranting a local $ϕS2S$ model for the Delta region.

We present a model for $ϕS2S$ that is M-dependent for each of the two locally calibrated site response models, which is consistent with the functional form proposed by Nweke et al. (2022):

where $ϕS2S,1$ is the site-to-site standard deviation for small M and is modeled for each site response model as follows:

where k is an indicator for which site response model is used (“v” or “$v,H/V$”), and ΔVar represents the change in variation from small-to-large magnitudes, and is modeled as follows:

where $b0$ and $b3$ represent the slope of $ϕS2S,1k$ and $ΔVar$ against T and ln(T), respectively; $b1$⁠, $b2$⁠, $b4$⁠, and $b5$ are constants; and $t1$–$t3$ represent transition periods between the four segments of each model. Values for $b0$–$b5$ and $t1$–$t3$ are given in Table 3.

The dependence of $ϕS2S$ on the underlying site response model is contained within the $ϕS2S,1$ terms, which are plotted in Figure 17a. The models described by Equations 23 and 24 are shown by the black curves. Computed ΔVar values from both models suggest similar M-dependence effects, as shown in Figure 17b. Accordingly, a single smoothed ΔVar model is presented, which is described by Equation 24. The patterns contained in both models suggest that variability is largest at intermediate periods, which coincides with the range of fundamental site periods at Delta sites. ΔVar in the Delta is less than what Nweke et al. (2022) found in southern California.

Site conditions in much of the Delta region of California are appreciably softer than the recommended range for global ergodic site response models, such as SS14. Not surprisingly, extrapolation of global models to these softer site conditions leads to biased predictions and is subject to large epistemic uncertainty. In this article, linear site response models specific to Delta sites and the immediate surrounding region are developed to facilitate more reliable ground motion predictions for hazard and risk studies.

We present local site response models conditioned on V_S30 and site parameters derived from mHVSR (c₀, a_p, f_p, and $mHVSR¯$⁠). The proposed V_S30-scaling model (Equation 7) can be used for any site condition encountered in the Delta, and is a standalone model. Values for period-dependent coefficients are provided in Table S1 of the supplemental content. Higher-order site effects correlated with mHVSR attributes are modeled as separate terms (Equations 11 and 17), and are used in addition to the V_S30-scaling model. Aleatory variability models for $ϕS2S$ were proposed, which should be paired with the selected mean model in applications (i.e. V_S30 only or V_S30 and mHVSR-informed). These models can be used with existing models from literature for $ϕSS$ (Goulet et al., 2021) and τ (BSSA14).

These models were developed using data from 36 seismic stations with V_S30 between 100 and 390 m/s, peat thicknesses between 0.4 and 10.1 m, and mHVSR peak frequencies (f_p) between 0.6 and 4.0 Hz. Bias could be expected for sites in the study region if they possess site characteristics outside the range of those used during model development. The Delta-calibrated local models should not be considered as applicable to other soft-soil regions without validation. However, the modeling approach outlined herein can be applied to other regions with unusual geologic conditions that may substantially impact site response, or can be generalized for ergodic model development.

This article highlights the utility of mHVSR as a predictor of site response. For Delta sites, the presence or lack thereof of peak features in mHVSR provided a useful indicator of when resonance effects, which act over a relatively narrow range of periods, were manifest. Moreover, we found that site response resonances could be modeled by a parameterization of mHVSR peak features, and general levels of residual site response could be modeled with simple broadband adjustments based on average mHVSR amplitudes. These types of models are advantageous because their functional forms can be easily understood and applied in forward application. It should not be misconstrued that mHVSR on its own can replace V_S30, which captures the average response over a wide period range and is generally readily available. Rather, our recommendation is to pair mHVSR with V_S30 and other site parameters (e.g. z_1.0) to improve site response estimates.

The authors thank Fabian Bonilla, Jim Kaklamanos, and one anonymous reviewer for their input on this article.