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

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 (VS30). 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 (fp), peak amplitude (ap), and average mHVSR amplitude over some frequency bandwidth ( 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 fp is used to model resonance effects, while a logistic function whose amplitude is correlated with 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 VS30-scaling model. Relative to a global ergodic model, the VS30-scaling model does not appreciably change the site-to-site aleatory variability ( ϕ S 2 S ) for periods shorter than about 1.5 s but reduces ϕ S 2 S by ~0.1 (natural log units) at long periods. When the mHVSR-informed model components are used, ϕ S 2 S 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 andAtkinson, 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.

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

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.
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 1 is peak amplitude, w p is peak width, c 0 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 0 and c 1 .The fit operation estimates f p , c 0 , c 1 , 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 0 ) for various reasons, including (1) f 0 could be out of the mHVSR usable frequency range for sites in deep basins, and (2) horizontal amplitudes near f 0 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: where ln Y ij À Á is the natural logarithm of the observed IM at site j from event i, and m 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 c k represents the model bias for GMM k; h E, ik and h S, jk represent random effect terms that quantify the event-and site-specific biases, respectively; and e 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): where the first and second terms represent the linear (f 1 ) and nonlinear amplification, respectively; x 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 IMr ) space for x 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 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: where f 1 ð Þ o j 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 h 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.

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 f 1 ð Þ o j 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 f 1 ð Þ o j over within-event residuals (h S, j + e ij ) during model development is that f 1 ð Þ o j 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, f 1 ð Þ o j 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 f 1 ð Þ o j 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.
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 f 1 ð Þ o j 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 underpredicts 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  Based on these observations, we modify SS14 as follows: where c 1 , c 2 , and c are V S30 -scaling coefficients; V 1 , V 2 , 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 2 ).The coefficients c and V c are adopted from SS14, and c 1 , c 2 , V 1 , and V 2 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 1 and c 2 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 1 is constrained to be less than V 2 , 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 1 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 2 , V 1 , and V 2 estimated as free parameters.These results indicated that V 1 was the most stable parameter.Regressed values of V 1 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 2 and c 2 , 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 f 1 ð Þ o j 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 andAtkinson (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 v ij ) are computed using the predictions given by Equation 8 and partitioned using mixed-effects methods to estimate site terms (h v S, j ), which represent the difference between site-specific observed amplification and the GMM described by Equation 8.These site terms (h v S, j ) 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 nonoverlapping 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 h v S, j 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 h v S, j 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.
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 0 ).This is shown in Figure 8a where sites with and without relative site response peaks are plotted in c 0 -a p space.Using logistic regression (Lottes et al., 1996), we propose a probabilistic model conditioned on c 0 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), b 0 = 19.24716 11.2033, b 1 = 3.8467 6 3.0136, and b 2 = 4.3943 6 2.3088.Figure 8a illustrates the variation of P p in c 0 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 0 = 1.150, a p = 3.766); as shown in Figure 8b, this site was not identified as having a h v S, j peak but has predicted P p = 0.85.Different analysts may have different opinions on whether or not YU_HOL1 has a h v S, j 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 h v S, j peak, and if P p .0.50 there is a h v S, j peak), the predictive model correctly predicts h v S, j 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 (h v S, j ) peaks for a given site j with a pulse function described as where f p1, j represents the peak frequency of the pulse, a 1, j is the amplitude of the pulse, v 1, j is the width of the pulse, b j represents the amplification level at short periods, and d 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 À1 p1, j , and the second term is provided to capture relatively flat responses at longer periods.Term h v S, j (f À1 p1, j ) represents the model prediction at the peak period (f À1 p1, j ).Example fits of this pulse function are given in Figure 9.Note at this stage b j is taken as a constant (period-independent), despite variations of h v S, j 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 f p1, j , a 1, j , b j , v 1, j , and d 1, j to mHVSR peak site parameters (f p , c 0 , c 1 , 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 b j and d 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 Figure 8.(a) Occurrences and non-occurrences of relative site response peaks as a function of mHVSR tail amplitude (c 0 ) and peak amplitude (a p =c 0 + c 1 ) (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, h v S, j (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).
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.
The updated GMM, which uses the V S30 -scaling model described by Equation 7and 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 0 , 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 (dW v, p ij ) are computed using the data and predictions given by Equation 15 and are partitioned using mixed-effects methods to estimate site terms (h v, p S, j ).Those site terms represent the difference between site-specific observed amplification and the GMM described by Equation 15.The site response represented by h v, p S, j is referred to as residual site response, which captures remaining features after V S30 -scaling and peakresonance effects have been accounted for.
General trends of h v, p S, j 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: 1.5527 Hz 0.4221 h v, p S, j = Da where a 0 is a constant which quantifies the average long period amplification in natural log units, Da 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 h v, p S, j (m h ) across short-and long-period bands to average mHVSR amplitudes (mHVSR) across different frequency bands.For example, we compare m h 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 m h and mHVSR, and consistently observed meaningful correlation for short period m h with mHVSR, but did not observe meaningful correlations for long period m h (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 a 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 a 2 depends on mHVSR as where a 8 = -0.26966 0.0770, a 9 = 0.5760 6 0.2122, a 10 = 0.7435 6 0.2516, a 11 = 0.1334 6 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 h v, p S, j 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 0 , 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.

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 0 , 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
Model bias is evaluated as the mean misfit of site terms (h S, j ) computed when implementing each of the site response models.Figure 14 presents plots of average h 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 0 , 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 overpredicts 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 0 , 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 (s), between-event (t), and within-event (f) standard deviations, which are related as The within-event term can be partitioned into site-to-site (f S2S ) and single-station (f 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 (f S2S and f SS ) for the Delta subset for most oscillator periods.Delta f SS values are sometimes larger than the values computed for the California data set; however, their 95% confidence intervals generally capture the ergodic global f SS model values proposed by Goulet et al. (2021).Accordingly, we do not consider a Delta-specific f 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 (f 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 f 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 f S2S was found to be M-dependent, as shown in Figure 16.We compare computed mean values of f S2S and their 95% confidence intervals to the NGA-West2 results reported by Gea18.The local f 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 f 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 f S2S model for the Delta region.
We present a model for f 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 f 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 DVar represents the change in variation from small-to-large magnitudes, and is modeled as follows: where b 0 and b 3 represent the slope of f k S2S, 1 and DVar 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.
The dependence of f S2S on the underlying site response model is contained within the f S2S, 1 terms, which are plotted in Figure 17a.The models described by Equations 23 and 24 are shown by the black curves.Computed DVar values from both models suggest similar M-dependence effects, as shown in Figure 17b.Accordingly, a single smoothed DVar 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.DVar in the Delta is less than what Nweke et al. (2022) found in southern California.

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 0 , 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 f 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 f SS (Goulet et al., 2021) and t (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 Deltacalibrated 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.

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

Figure 2 .
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.25fc,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.

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

Figure 4 .
Figure 4. Plots of observed linear mean amplification, f 1 ð Þ o j , 6one standard error (SE) versus oscillator period for stations with reliable h S, j estimates for all periods < 9 s, ordered by increasing V S30 from topleft to bottom-right.The shading around the mean indicates 6one standard error.
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.

Figure 6 .
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.

Figure 7 .
Figure 7. Plots of relative site response, h v S, j (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.

Figure 9 .
Figure 9. Relative site response (h v S, j ) 6 one standard error and example fits of Equation 10 to two Delta sites with resonance peaks: (a) CE 67587 and (b) NP LVB4.

Figure 10 .
Figure 10.Relationships of mHVSR peak frequency (f p ) and relative site response resonance model coefficients (a) fp , (b) â1 , and (c) v.

Figure 11 .
Figure 11.Plots of residual site response, h v, p S, j (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).

Figure 12 .
Figure 12.(a) Relationship between average residual site response (m h ) computed for the 0.01-3 s period range versus average mHVSR amplitude (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 ).

Figure 13 .
Figure 13.Plots of observed mean linear amplification, f 1 ð Þ 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).

Figure 14 .
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%).

Figure 15 .
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 (s), (b) between-event variability (t), (c) within-event variability (f), (d) site-to-site variability (f S2S ), and (e) single-station within-event variability (f SS ); results are shown for all California data (solid) and Delta subset used in this study (dashed).

Figure 17 .
Figure 17.(a) Site-to-site standard deviations (f S2S ) for Delta data with M < 5 calculated using the proposed Delta-specific linear site response models, and their recommended f S2S, 1 models; (b) DVar results and recommended model.

Table 1 .
Metadata for the 36 Delta stations used in this study b NA = No mHVSR available.

Table 2 .
Estimated values and standard errors for coefficients in Equations 12-14

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