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

Abstract

Sites located in the Sacramento-San Joaquin Delta region of California typically have peaty-organic soils near the ground surface, which are characteristically soft, with shear wave velocities as low as 30 m/s. These unusually soft geotechnical conditions, which are outside the range of applicability of existing ergodic site amplification models, can be anticipated to produce significant site effects during earthquake shaking. We evaluate site response for 36 seismic stations in the Delta region using non-ergodic methods with low-amplitude ground motion data. We model first-order site effects using a period-dependent relation conditioned on the 30 m time-averaged shear wave velocity (V_S30). Relative to extrapolations of global ergodic models, this Delta-specific model provides lower and higher levels of amplification for short and long periods, respectively. While the local model provides unbiased predictions for Delta sites as a whole, it smooths over site-specific effects such as resonant peaks. Microtremor-based horizontal-to-vertical spectral ratios (mHVSRs) were measured for 34 sites, from which additional site parameters such as peak frequency (f_p), peak amplitude (a_p), and average mHVSR amplitude over some frequency bandwidth ( $\overline{\mathrm{mHVSR}}$ ) are derived. Sites with prominent mHVSR peaks are found to often exhibit site resonance effects, while sites without prominent peak features generally do not. A modified Ricker wavelet pulse function conditioned on f_p is used to model resonance effects, while a logistic function whose amplitude is correlated with $\overline{\mathrm{mHVSR}}$ is used to model general broadband amplification that transitions to zero at long periods. These mHVSR-informed models are implemented as additive components to the V_S30-scaling model. Relative to a global ergodic model, the V_S30-scaling model does not appreciably change the site-to-site aleatory variability ( ${\varphi }_{S2S}$ ) for periods shorter than about 1.5 s but reduces ${\varphi }_{S2S}$ by ~0.1 (natural log units) at long periods. When the mHVSR-informed model components are used, ${\varphi }_{S2S}$ is reduced by ~0.05 to 0.1 for short-to-intermediate periods (up to ~2 s). A model for nonlinear effects, derived from nonlinear ground response analyses, will be presented in a separate paper.

Introduction

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).

Figure 1.

(a) Locations of earthquakes in California and western Nevada recorded by Delta stations (boundary of Central Valley is from Buckreis et al., 2023b); (b) inset map showing locations of Delta stations within the Delta jurisdictional boundaries given by Department of Water Resources (2018); and (c) distribution of V_S30 from “Delta” (yellow) and “non-Delta” (blue) stations used in this study, from V_S profiles measured in the Delta irrespective of whether they are located near a ground motion station (“Delta V_S Profile”; red-outline), and from all California sites (gray) in the national database (Kwak et al., 2021).

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.

Data sources

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.

Ground motions and related metadata

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.

Figure 2.

(a) Distribution of Delta ground motion records in magnitude versus distance space; (b) distribution of median-component (RotD50) peak ground accelerations (PGA); (c) number of usable records (T < 1/(1.25f_c,hp), where f_c,hp is the greater of the two high-pass corner frequencies selected during signal processing) and usable stations (each with at least four records) versus period.

Site conditions

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.

Time-averaged shear wave velocity in the upper 30 m (V_S30)

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.

Table 1.

Metadata for the 36 Delta stations used in this study

Station code	Latitude (o)	Longitude (o)	VS30 (m/s) a	No. of mHVSR peaks b	c0	c1	fp (Hz)	wp	No. of records
BK_JEPS	38.2579	−121.8252	(389.8)	1	0.833	1.457	0.661	0.495	20
BK_TWIT	38.0971	−121.6832	(157.6)	1	0.944	2.940	0.691	0.313	21
CE_57195	37.9753	−121.3140	306.2	0	—	—	—	—	19
CE_57531	37.9332	−121.6956	(250.6)	2	0.814	1.507	0.750	0.500	18
CE_57534	37.9119	−121.6219	(238.3)	0	—	—	—	—	28
CE_67215	37.9925	−121.6384	(227.5)	1	1.100	0.914	0.695	0.383	9
CE_67533	38.0155	−121.6396	(188.4)	0	—	—	—	—	16
CE_67587	38.0390	−121.8986	(226.8)	2	1.100	4.150	1.105	0.289	20
CE_67615	38.0005	−121.7830	(287.1)	0	—	—	—	—	13
NC_C057	37.9555	−121.6979	(344.9)	0	—	—	—	—	13
NP_DIX	38.3771	−121.8424	284.8	1	1.200	1.026	0.689	0.313	39
NP_EMR	38.0605	−121.4993	(134.1)	1	1.000	3.830	1.030	0.285	17
NP_KIR	38.0551	−121.4582	(158.5)	1	1.456	7.900	2.082	0.347	29
NP_LVA3	38.0328	−121.7631	(154.6)	1	1.750	3.430	1.000	0.251	6
NP_LVA4	38.0335	−121.7634	(139.8)	2	1.106	5.088	0.855	0.343	11
NP_LVB3	37.9102	−121.5649	(195.0)	NA	NA	NA	NA	NA	4
NP_LVB4	37.9097	−121.5656	(190.1)	NA	NA	NA	NA	NA	6
NP_MCD	37.9802	−121.4735	(105.8)	0	—	—	—	—	11
NP_PLA	37.7987	−121.4632	(184.7)	0	—	—	—	—	29
NP_SIA	38.0503	−121.7367	135.8	0	—	—	—	—	38
WR_CKR	38.3145	−121.4920	(319.9)	0	—	—	—	—	26
WR_CLFN	37.8570	−121.5730	160.6	0	—	—	—	—	9
WR_CLFS	37.8273	−121.5604	185.2	0	—	—	—	—	14
WR_HOLT	37.9551	−121.4227	(257.0)	1	1.230	1.438	0.647	0.320	7
WR_MONN	38.0929	−121.8855	152.0	2	1.230	4.248	1.314	0.176	7
WR_SHER	38.0306	−121.7450	122.4	1	1.152	4.455	0.732	0.402	7
WR_SIFF	38.0311	−121.7450	(103.4)	1	1.019	2.645	0.735	0.245	4
WR_STNI	38.1197	−121.5401	165.8	1	1.790	1.320	1.822	0.100	4
YU_CEC	37.9381	−121.5968	228.3	0	—	—	—	—	9
YU_HMT	37.9749	−121.5813	(191.5)	1	1.290	2.856	1.553	0.239	4
YU_HOL1	37.9748	−121.5818	(134.8)	1	1.150	2.616	2.179	0.385	4
YU_HOL2	37.9733	−121.5823	220.7	1	1.190	4.151	2.124	0.266	9
YU_HOL3	37.9733	−121.5844	235.5	1	1.050	4.195	4.653	0.272	4
YU_SMT	37.9982	−121.6252	(136.8)	1	1.520	1.150	3.158	0.200	5
YU_SRT	38.0107	−121.6246	(166.5)	1	1.592	2.804	2.879	0.223	13
YU_WHR	37.8984	−121.4508	197.0	1	0.869	1.076	3.071	0.375	13

mHVSR: microtremor-based horizontal-to-vertical spectral ratio.

a

Values in parenthesis are estimated. Values not in parenthesis are measured.

b

NA = No mHVSR available.

Microtremor-based horizontal-to-vertical spectral ratios (mHVSR)

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):

$$F_{\frac{H}{V},i}=c_{0,i}+c_{1,i}\exp [-\frac{1}{2}{\left(\frac{\ln \left(\frac{f}{f_{p,i}}\right)}{2w_{p,i}}\right)}^2]$$

(1)

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.

Figure 3.

Examples of peak identification and fitting algorithm proposed by Wang et al. (2023) for Delta sites with (a) no peak (WR CLFS) and (b) clear peak (CE 67265).

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.

Ground motion data analysis

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.

Residual calculations

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

$$R_{\mathrm{ij}}=\ln \left(Y_{\mathrm{ij}}\right)-{\mu }_{\mathrm{ij}}({\mathbf{M}}_i,R_{\mathrm{JB},\mathrm{ij}},V_{S30,j},z_{1,j})$$

(2)

where $\ln \left(Y_{\mathrm{ij}}\right)$ is the natural logarithm of the observed IM at site j from event i, and ${\mu }_{\mathrm{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):

$$R_{\mathrm{ijk}}=c_k+{\eta }_{E,\mathrm{ik}}+{\eta }_{S,\mathrm{jk}}+{\epsilon }_{\mathrm{ijk}}$$

(3)

where $c_k$ represents the model bias for GMM k; ${\eta }_{E,\mathrm{ik}}$ and ${\eta }_{S,\mathrm{jk}}$ represent random effect terms that quantify the event- and site-specific biases, respectively; and ${\epsilon }_{\mathrm{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.

Non-ergodic site response for Delta region stations

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):

$$F_S=f_1+f_2\ln \left(\frac{x_{\mathrm{IMr}}+f_3}{f_3}\right)$$

(4)

where the first and second terms represent the linear ( $f_1$ ) and nonlinear amplification, respectively; $x_{\mathrm{IMr}}$ 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); $f_2$ represents the slope in ln(amplification)-ln( $x_{\mathrm{IMr}}$ ) space for $x_{\mathrm{IMr}}>>f_3$ ; and $f_3$ 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 $x_{\mathrm{IMr}}$ 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 ( $f_1$ ) was estimated as the sum of the site term (established from weak motion data) and the linear portion of the site response GMM:

$${\left(f_1\right)}_j^o={\eta }_{S,j}+[F_{S,j}-f_2\ln \left(\frac{x_{\mathrm{IMr}}+f_3}{f_3}\right)]$$

(5)

where ${\left(f_1\right)}_j^o$ 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 $f_1$ 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 ${\eta }_{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.

Figure 4.

Plots of observed linear mean amplification, ${\left(f_1\right)}_j^o$ , ±one standard error (SE) versus oscillator period for stations with reliable ${\eta }_{S,j}$ estimates for all periods ≤ 9 s, ordered by increasing V_S30 from top-left to bottom-right. The shading around the mean indicates ±one standard error.

Site amplification model development

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 ${\left(f_1\right)}_j^o$ 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 ${\left(f_1\right)}_j^o$ over within-event residuals ( ${\eta }_{S,j}+{\epsilon }_{\mathrm{ij}}$ ) during model development is that ${\left(f_1\right)}_j^o$ 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, ${\left(f_1\right)}_j^o$ 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 ${\left(f_1\right)}_j^o$ 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:

$$w_j=\frac{n{\left(S{E}_j\right)}^{-2}}{\sum _{j=1}^n{\left(S{E}_j\right)}^{-2}}$$

(6)

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.

Stage 1: Analysis of V_S30-scaling

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 ${\left(f_1\right)}_j^o$ 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.

Figure 5.

V_S30-scaling for observed linear amplification, ${\left(f_1\right)}_j^o$ , and predictions from SS14 (black) and the proposed local model (red) for PGA, PGV, and PSA at 0.05, 0.3, 1.0, and 5.0 s.

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:

$$F_{lin}\left(V_{S30}\right)=\{\begin{array}{cc}{c}_1\ln \left(\frac{V_{S30}}{V_1}\right)+c_2\ln \left(\frac{V_1}{V_2}\right)+c\ln \left(\frac{V_2}{V_{ref}}\right)& V_{S30}<V_1\\ c_2\ln \left(\frac{V_{S30}}{V_2}\right)+c\ln \left(\frac{V_2}{V_{ref}}\right)& V_1\le V_{S30}<V_2\\ c\ln \left(\frac{V_{S30}}{V_{ref}}\right)& V_2\le V_{S30}<V_c\\ c\ln \left(\frac{V_c}{V_{ref}}\right)& V_c\le V_{S30}\end{array}$$

(7)

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.

Figure 6.

As-regressed (discrete symbols) and recommended-smoothed coefficients (lines) for (a) V_S30-scaling coefficients and (b) limiting velocities. The impacts of the smoothing operation are indicated by the differences between symbols and lines.

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

$${\mu }_{\mathrm{ij}}^v({\mathit{M}}_i,R_{\mathrm{JB},\mathrm{ij}},V_{S30,j})={\mu }_{\mathrm{ij}}({\mathit{M}}_i,R_{\mathrm{JB},\mathrm{ij}},V_{S30}=760m/s,\delta z_1=0)+F_{\mathrm{lin},j}\left(V_{S30,j}\right)+F_{\mathrm{nl},\mathrm{ij}}$$

(8)

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 ${\left(f_1\right)}_j^o$ 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.

Stage 2: Analysis of site resonance effects

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.

Identification of site resonances

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 ( $R_{\mathrm{ij}}^v$ ) are computed using the predictions given by Equation 8 and partitioned using mixed-effects methods to estimate site terms ( ${\eta }_{S,j}^v$ ), which represent the difference between site-specific observed amplification and the GMM described by Equation 8. These site terms ( ${\eta }_{S,j}^v$ ) 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 ${\eta }_{S,j}^v$ 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 ${\eta }_{S,j}^v$ 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.

Plots of relative site response, ${\eta }_{S,j}^v$ (black), versus (oscillator) period for stations shown in Figure 1. Overlain on the figure are plots of mHVSR (blue) versus period (inverse of frequency) and fitted peaks (red curves; no curve indicates lack of peak). The presence of resonance peaks in relative site response is indicated in the upper right of each plot by “P” for peak, or “NP” for no-peak.

Predicting site resonances using mHVSR

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:

$$P_p(c_0,a_p)={[1+\exp ({\beta }_0-{\beta }_1{c}_0-{\beta }_2{a}_p)]}^{-1}$$

(9)

where P_p is the probability of a resonant peak (0 ≤ P_p ≤ 1), ${\beta }_0$  = 19.2471 ± 11.2033, ${\beta }_1$  = 3.8467 ± 3.0136, and ${\beta }_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 ${\eta }_{S,j}^v$ peak but has predicted P_p = 0.85. Different analysts may have different opinions on whether or not YU_HOL1 has a ${\eta }_{S,j}^v$ peak; however, the algorithm did not identify a peak because the relative uncertainty is quite large (the site recorded only four events).

Figure 8.

(a) Occurrences and non-occurrences of relative site response peaks as a function of mHVSR tail amplitude (c₀) and peak amplitude (a_p =c₀ + c₁) (discrete symbols) and predictions of logistic model (lines) for peak probabilities P_p = 0.15, 0.5, and 0.85. The shaded region corresponds to conditions that are not possible for mHVSR peaks. (b) Plot of relative site response, ${\eta }_{S,j}^v$ (black), versus (oscillator) period for station YU HOL1, which is an outlier. Overlain on the plot is mHVSR (blue) versus period (inverse of frequency) and fitted peak (red).

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 ${\eta }_{S,j}^v$ peak, and if P_p > 0.50 there is a ${\eta }_{S,j}^v$ peak), the predictive model correctly predicts ${\eta }_{S,j}^v$ peaks from mHVSR approximately 91% of the time, which is considered to be a sufficient indicator for application purposes.

Using mHVSR peak parameters to predict relative site response resonances

For modeling purposes, we parameterize relative site response ( ${\eta }_{S,j}^v$ ) peaks for a given site j with a pulse function described as

$${\overline{\eta }}_{S,j}^v\left(T\right)=\{\begin{array}{cc}\frac{2{\alpha }_{1,j}}{\sqrt{3{\omega }_{1,j}}{\pi }^{1/4}}[1-{\left(\frac{\ln \left(T{f}_{p1,j}\right)}{{\omega }_{1,j}}\right)}^2]\exp {[-\frac{1}{2}\left(\frac{\ln \left(f_{p1,j}\right)}{{\omega }_{1,j}}\right)]}^2+{\beta }_j\hfill & T\le f_{p1,j}^{-1}\\ [{\overline{\eta }}_{S,j}^v\left(f_{p1,j}^{-1}\right)-{\delta }_{1,j}]\exp [-{\left(\frac{\ln \left(T{f}_{p1,j}\right)}{2{\omega }_{1,j}/3}\right)}^2+{\delta }_{1,j}]\hfill & T>f_{p1,j}^{-1}\end{array}$$

(10)

where $f_{p1,j}$ represents the peak frequency of the pulse, ${\alpha }_{1,j}$ is the amplitude of the pulse, ${\omega }_{1,j}$ is the width of the pulse, ${\beta }_j$ represents the amplification level at short periods, and ${\delta }_{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 $f_{p1,j}^{-1}$ , and the second term is provided to capture relatively flat responses at longer periods. Term ${\overline{\eta }}_{S,j}^v(f_{p1,j}^{-1})$ represents the model prediction at the peak period ( $f_{p1,j}^{-1})$ . Example fits of this pulse function are given in Figure 9. Note at this stage ${\beta }_j$ is taken as a constant (period-independent), despite variations of ${\eta }_{S,j}^v$ with period (outside of the width of the pulse) that may be present in some cases.

Figure 9.

Relative site response ( ${\eta }_{S,j}^v$ ) ± one standard error and example fits of Equation 10 to two Delta sites with resonance peaks: (a) CE 67587 and (b) NP LVB4.

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 $f_{p1,j}$ , ${\alpha }_{1,j}$ , ${\beta }_j$ , ${\omega }_{1,j}$ , and ${\delta }_{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 ${\beta }_j$ and ${\delta }_{1,j}$ coefficients are dropped because they show no apparent correlation to any predictor variable, and (2) the amplitude is scaled by P_p:

$$f_{1,p}(T,f_p)=\{\begin{array}{cc}{P}_p\times \frac{\widehat{2{\alpha }_1}}{\sqrt{3\hat{\omega }}{\pi }^{1/4}}[1-{\left(\frac{\ln \left(T{\hat{f}}_p\right)}{\hat{\omega }}\right)}^2]\exp {[-\frac{1}{2}\left(\frac{\ln \left({\hat{f}}_p\right)}{\hat{\omega }}\right)]}^2\hfill & T\le f_p^{-1}\\ P_p\times f_{1,p}(T=f_p^{-1},f_p)\exp [-{\left(\frac{\ln \left(T{\hat{f}}_p\right)}{2\hat{\omega }/3}\right)}^2]\hfill & T>f_p^{-1}\end{array}$$

(11)

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

$${\hat{f}}_p\left(f_p\right)=\exp [a_0\ln \left(f_p\right)+a_1]$$

(12)

$${\hat{\alpha }}_1\left(f_p\right)=\{\begin{array}{cc}{a}_2\ln \left(f_p\right)+a_3\hfill & f_p<f_{L1}\\ a_4\hfill & f_p\ge f_{L1}\end{array}$$

(13)

and

$$\hat{\omega }\left(f_p\right)=\{\begin{array}{cc}{a}_5\ln \left(f_p\right)+a_6\hfill & f_p<f_{L2}\\ a_7\hfill & f_p\ge f_{L2}\end{array}$$

(14)

where f _p is in Hz; coefficients $a_0$ , $a_2$ , and $a_5$ control the scaling of their respective terms with $f_p$ ; $a_1$ , $a_3$ , $a_4$ , $a_6$ , and $a_7$ are constants; and $f_{L1}$ – $f_{L2}$ are limiting frequencies that represent the transition from scaled-behavior (i.e. the fitted parameter scales linearly with $f_p$ ) to constant-behavior. The values for coefficients $a_0$ – $a_7$ and $f_{L1}$ – $f_{L2}$ are given in Table 2.

Figure 10.

Relationships of mHVSR peak frequency (f_p) and relative site response resonance model coefficients (a) ${\hat{f}}_p$ , (b) ${\hat{\alpha }}_1$ , and (c) $\hat{\omega }$ .

Table 2.

Estimated values and standard errors for coefficients in Equations 12–14

Coefficient	Estimate	SE
$a_0$	0.9265	0.0212
$a_1$	0.0978	0.0235
$a_2$	−0.1790	0.0420
$a_3$	0.2355	0.0157
$a_4$	0.1044	0.0465
$a_5$	−0.3378	0.1467
$a_6$	0.5213	0.0344
$a_7$	0.3733	0.0987
$f_{L1}$	2.0821 Hz	0.4934
$f_{L2}$	1.5527 Hz	0.4221

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

$$\begin{array}{c}\multicolumn{1}{c}{{\mu }_{\mathrm{ij}}^{v,p}({\mathit{M}}_i,R_{\mathrm{JB},\mathrm{ij}},V_{S30,j},f_{p,j},c_{0,j},a_{p,j})={\mu }_{\mathrm{ij}}({\mathbf{M}}_i,R_{\mathrm{JB},\mathrm{ij}},V_{S30}=760\mathrm{m}/\mathrm{s},\delta z_1=0)+F_{\mathrm{lin},j}\left(V_{S30,j}\right)}\\ +f_{1,p,j}(T,f_{p,j},c_{0,j},a_{p,j})+F_{\mathrm{nl},\mathrm{ij}}\hfill \end{array}$$

(15)

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.

Stage 3: Analysis of remaining amplification 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 ( $\delta W_{\mathrm{ij}}^{v,p}$ ) are computed using the data and predictions given by Equation 15 and are partitioned using mixed-effects methods to estimate site terms ( ${\eta }_{S,j}^{v,p}$ ). Those site terms represent the difference between site-specific observed amplification and the GMM described by Equation 15. The site response represented by ${\eta }_{S,j}^{v,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 ${\eta }_{S,j}^{v,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:

$${\overline{\eta }}_{S,j}^{v,p}=\frac{\mathrm{\Delta }\alpha }{1+\exp [1.923\ln \left(\frac{T}{T_c}\right)]}+{\alpha }_0$$

(16)

where ${\alpha }_0$ is a constant which quantifies the average long period amplification in natural log units, $\mathrm{\Delta }\alpha$ 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.

Figure 11.

Plots of residual site response, ${\eta }_{S,j}^{v,p}$ (black), versus (oscillator) period for stations shown in Figure 1. Fits of Equation 16 are represented by the red curves. Overlain on the figure are plots of mHVSR (blue) versus period (inverse of frequency).

For modeling purposes, we relate the arithmetic average of ${\eta }_{S,j}^{v,p}$ ( ${\mu }_{\eta }$ ) across short- and long-period bands to average mHVSR amplitudes ( $\overline{\mathrm{mHVSR}}$ ) across different frequency bands. For example, we compare ${\mu }_{\eta }$ computed over the period range between 0.01 and 3 s and $\overline{\mathrm{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 ${\mu }_{\eta }$ and $\overline{\mathrm{mHVSR}}$ , and consistently observed meaningful correlation for short period ${\mu }_{\eta }$ with $\overline{\mathrm{mHVSR}}$ , but did not observe meaningful correlations for long period ${\mu }_{\eta }$ (e.g., between 1 and 10 s) with $\overline{\mathrm{mHVSR}}$ .

Figure 12.

(a) Relationship between average residual site response ( ${\mu }_{\eta }$ ) computed for the 0.01–3 s period range versus average mHVSR amplitude ( $\overline{\mathrm{mHVSR}}$ ) computed for the 0.33–50 Hz frequency range; (b) schematic illustrating the period dependence of the mHVSR-based constant amplification adjustment (f_1,m).

Based on these observations, we propose a model (denoted f_1,m) conditioned on $\overline{\mathrm{mHVSR}}$ to capture this behavior:

$$f_{1,m}(T,\overline{\mathrm{mHVSR}})=\frac{{\alpha }_2}{1+\exp [1.923\ln \left(\frac{T}{1.5s}\right)]}$$

(17)

where ${\alpha }_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 ${\alpha }_2$ depends on $\overline{\mathrm{mHVSR}}$ as

$${\alpha }_2\left(\overline{\mathrm{mHVSR}}\right)=\{\begin{array}{cc}{a}_8\hfill & \overline{\mathrm{mHVSR}}<A_{L1}\\ a_9\overline{\mathrm{mHVSR}}-a_{10}\hfill & A_{L1}\le \overline{\mathrm{mHVSR}}\le A_{L2}\\ a_{11}\hfill & \overline{\mathrm{mHVSR}}>A_{L2}\end{array}$$

(18)

where $a_8$  = –0.2696 ± 0.0770, $a_9$  = 0.5760 ± 0.2122, $a_{10}$  = 0.7435 ± 0.2516, $a_{11}$  =0.1334 ± 0.0715, $A_{L1}$  = 0.8228, and $A_{L2}$  = 1.5224. The coefficients in Equation 18 correspond to the optimal bandwidths that produced the highest correlation when relating ${\overline{\eta }}_{S,j}^{v,p}$ and $\overline{\mathrm{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

$$\begin{array}{c}\multicolumn{1}{c}{{\mu }_{\mathrm{ij}}^{v,H/V}({\mathit{M}}_i,R_{\mathrm{JB},\mathrm{ij}},V_{S30,j},f_{p,j},c_{0,j},a_{p,j},{\overline{\mathrm{mHVSR}}}_j)={\mu }_{\mathrm{ij}}({\mathbf{M}}_i,R_{\mathrm{JB},\mathrm{ij}},V_{S30}=760m/s,\delta z_1=0)}\\ \multicolumn{1}{c}{+F_{\mathrm{lin},j}\left(V_{S30,j}\right)+f_{1,p,j}(T,f_{p,j},c_{0,j},a_{p,j})+f_{1,m,j}\left({\overline{\mathrm{mHVSR}}}_j\right)+F_{\mathrm{nl},\mathrm{ij}}}\end{array}$$

(19)

where the linear site response is a function of V_S30, f_p, c₀, a_p, and $\overline{\mathrm{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.

Model performance and comparison to global ergodic model

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₀, $\overline{\mathrm{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.

Figure 13.

Plots of observed mean linear amplification, ${\left(f_1\right)}^o$ , and model predictions provided by SS14 (ergodic model; gray) and two proposed local models for the Delta: simple V_S30-scaling (blue) and mHVSR-informed (red).

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

Model bias is evaluated as the mean misfit of site terms ( ${\eta }_{S,j}$ ) computed when implementing each of the site response models. Figure 14 presents plots of average ${\eta }_{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.

Figure 14.

Comparison of bias in average site amplification predicted by SS14 (ergodic model; gray) and two proposed local models for the Delta: simple V_S30-scaling (blue) and mHVSR-informed (red) for (a) all sites, (b) sites predicted to have minimal site resonance (P_p < 50%), and (c) sites predicted to have impactful site resonances (P_p ≥ 50%).

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₀, $\overline{\mathrm{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₀, $\overline{\mathrm{mHVSR}}$ ) model reduces these biases for both site types.

Aleatory variability

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

$$\sigma =\sqrt{{\tau }^2+{\varphi }^2}$$

(20)

The within-event term can be partitioned into site-to-site ( ${\varphi }_{S2S}$ ) and single-station ( ${\varphi }_{\mathrm{SS}}$ ) standard deviations as follows:

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

(21)

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 ( ${\varphi }_{S2S}$ and ${\varphi }_{\mathrm{SS}}$ ) for the Delta subset for most oscillator periods. Delta ${\varphi }_{\mathrm{SS}}$ values are sometimes larger than the values computed for the California data set; however, their 95% confidence intervals generally capture the ergodic global ${\varphi }_{\mathrm{SS}}$ model values proposed by Goulet et al. (2021). Accordingly, we do not consider a Delta-specific ${\varphi }_{\mathrm{SS}}$ model to be justified, which is consistent with the conclusions reached by Buckreis et al. (2023b) for the California data set.

Figure 15.

Standard deviations calculated using SS14 (ergodic model; black) and two proposed local models for the Delta: simple V_S30-scaling (blue) and mHVSR-informed (red) for (a) total variability ( $\sigma$ ), (b) between-event variability ( $\tau$ ), (c) within-event variability ( $\varphi$ ), (d) site-to-site variability ( ${\varphi }_{S2S}$ ), and (e) single-station within-event variability ( ${\varphi }_{\mathrm{SS}}$ ); results are shown for all California data (solid) and Delta subset used in this study (dashed).

Site-to-site variability ( ${\varphi }_{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 ${\varphi }_{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 ${\varphi }_{S2S}$ was found to be M-dependent, as shown in Figure 16. We compare computed mean values of ${\varphi }_{S2S}$ and their 95% confidence intervals to the NGA-West2 results reported by Gea18. The local ${\varphi }_{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.

Figure 16.

Site-to-site standard deviations ( ${\varphi }_{S2S}$ ) for global data calculated using BSSA14 (Goulet et al., 2018: Gea18; gray) and Delta data calculated using SS14 (ergodic model; black) and two proposed Delta-specific local models: simple V_S30-scaling (blue) and mHVSR-informed (red) for (a) M < 5 and (b) M ≥ 5.

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 ${\varphi }_{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 ${\varphi }_{S2S}$ model for the Delta region.

We present a model for ${\varphi }_{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):

$${\varphi }_{S2S}\left(\mathit{M}\right)=\{\begin{array}{cc}{\varphi }_{S2S,1}\hfill & \mathit{M}\le 5\\ \sqrt{{\varphi }_{S2S,1}^2-\mathrm{\Delta }\mathrm{Var}(\mathit{M}-5)}\hfill & 5<\mathit{M}<6\\ \sqrt{{\varphi }_{S2S,1}^2-\mathrm{\Delta }\mathrm{Var}}\hfill & \mathit{M}\ge 6\end{array}$$

(22)

where ${\varphi }_{S2S,1}$ is the site-to-site standard deviation for small M and is modeled for each site response model as follows:

$${\varphi }_{S2S,1}^k\left(T\right)=\{\begin{array}{cc}{b}_{0}T+b_1\hfill & T<t_1\\ b_2\hfill & t_1\le T<t_2\\ b_3\ln \left(T\right)+b_4\hfill & t_2\le T<t_3\\ b_5\hfill & t_3\le T\end{array}$$

(23)

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:

$$\mathrm{\Delta }\mathrm{Var}(T)=\{\begin{array}{cc}{b}_{0}T+b_1\hfill & T<t_1\\ b_2\hfill & t_1\le T<t_2\\ b_3\ln \left(T\right)+b_4\hfill & t_2\le T<t_3\\ b_5\hfill & t_3\le T\end{array}$$

(24)

where $b_0$ and $b_3$ represent the slope of ${\varphi }_{S2S,1}^k$ and $\mathrm{\Delta }\mathrm{Var}$ against T and ln(T), respectively; $b_1$ , $b_2$ , $b_4$ , and $b_5$ are constants; and $t_1$ – $t_3$ represent transition periods between the four segments of each model. Values for $b_0$ – $b_5$ and $t_1$ – $t_3$ are given in Table 3.

Table 3.

Estimated coefficient values for Equations 23 and 24 for each site response model

Model	$b_0$	$b_1$	$b_2$	$b_3$	$b_4$	$b_5$	$t_1$ (s)	$t_2$ (s)	$t_3$ (s)
${\varphi }_{S2S,1}^v$	0.1793	0.3211	0.3839	−0.1736	0.3809	0.2151	0.35	0.96	2.54
${\varphi }_{S2S,1}^{H/V}$	0.2731	0.2558	0.3144	−0.0964	0.2981	0.2127	0.27	0.68	1.95
$\mathrm{\Delta }\mathrm{Var}$	0.2012	0.0498	0.0817	−0.0750	0.0399	0.0009	0.16	0.57	1.69

The dependence of ${\varphi }_{S2S}$ on the underlying site response model is contained within the ${\varphi }_{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.

Figure 17.

(a) Site-to-site standard deviations ( ${\varphi }_{S2S}$ ) for Delta data with M ≤ 5 calculated using the proposed Delta-specific linear site response models, and their recommended ${\varphi }_{S2S,1}$ models; (b) $\mathrm{\Delta }\mathrm{Var}$ results and recommended model.

Conclusions and recommendations

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 $\overline{\mathrm{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 ${\varphi }_{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 ${\varphi }_{\mathrm{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.

Supplemental Material

sj-xlsx-1-eqs-10.1177_87552930231217165 – Supplemental material for Small-strain site response of soft soils in the Sacramento-San Joaquin Delta region of California conditioned on VS30 and mHVSR

Supplemental material, sj-xlsx-1-eqs-10.1177_87552930231217165 for Small-strain site response of soft soils in the Sacramento-San Joaquin Delta region of California conditioned on VS30 and mHVSR by Tristan E Buckreis, Jonathan P Stewart, Scott J Brandenberg and Pengfei Wang in Earthquake Spectra