All average velocities were calculated as the time-weighted average, i.e., depth to the horizon of interest divided by the total travel time for a vertically ascending S-wave, where the total travel time is the sum of the layer travel times. To calculate layer impedances, we first estimated layer densities (ρ) using the ρ(Vs) relationship in Boore (2016). Then impedance ratios were calculated at each Vs-profile interface as the ratio of the impedance of the underlying layer to average impedance of the overlying layers (Wang and Carpenter, 2023): I R j = V s j − ∙ ρ j − V s j + ¯ ∙ ρ j + ¯ (3) where V s j + ¯ and ρ j + ¯ are the time-averaged Vs and thickness-average ρ of the layers above the j ^th interface and V s j − and ρ j − are the Vs and ρ of the layer immediately below the j ^th interface. IRmax is defined as the largest impedance ratio at a site, Z _IRmax is the depth to that interface, and Vs _IRmax is time-averaged Vs above that interface.

Using the collected Vs-profiles as input, we calculated 1D site response using the matrix propagation method of Haskell (1953); Haskell, (1960). This methodology has been used widely and in the CEUS (Carpenter et al., 2018; Carpenter et al., 2020) to calculate one-dimensional linear elastic SH-wave amplification functions for input velocity profiles and other dynamic parameters. Thus, we used this method to calculate theoretical site responses at each site in this investigation. We modeled vertically incident SH-waves, incorporating viscoelastic effects through complex shear moduli in each layer, similar to Joyner et al. (1976), and we refer to the resultant theoretical amplification functions as the 1D responses.

The density profiles developed to calculate the impedance ratio profiles were also used in the 1D site response calculations. Viscoelastic effects were accounted for using Q _ef , the effective S-wave quality factor. We follow Campbell (2009), who interpreted Q _ef in eastern North America as frequency-independent and calculated as the inverse of the summed intrinsic and scattering S-wave attenuations in near-surface layers. We used relationships between Q _ef and Vs developed for the central and eastern U.S. by Wang et al. (1994), equations 20 and 21, therein, and Campbell (2009), equations 14 and 15 therein, to estimate Q _ef for each layer in each Vs profile. We then calculated γ s , the S-wave damping ratio, for each layer from Q e f using (Campbell, 2009) γ s = 1 2   Q e f . (4)

Carpenter et al. (2020) demonstrated that the range of γ _s estimated from the four Q _ef (Vs) models and Equation 4 is small at the 11 CEUS sites they evaluated, which were on a variety of underlying geologies. Thus, we calculated the final γ s for each layer at each site, except the two deep vertical seismic arrays in the upper Mississippi Embayment, VSAP and CUSSO, by applying Equation 4 on the arithmetic mean of the four Q _ef (Vs) estimates determined for a particular layer. For VSAP and CUSSO (further discussed below), we estimated Q _ef for each layer using equations 20 and 21 in Wang et al. (1994), which respectively were developed for those sites, and then applying Equation 4.

We extracted the first and peak frequencies, f_0,1D and f_p,1D , and corresponding amplifications, A_0,1D and A_p,1D , from the 1D responses, as shown in Figure 1, for evaluation.

Earthquake HVSR, which has been shown to approximate site response at some sites (Lermo and Chávez-García, 1993; Sedaghati et al., 2018; Carpenter et al., 2020), is the ratio of the horizontal-component amplitude spectrum divided by the vertical-component’s amplitude spectrum calculated from seismograms of local and regional earthquakes. For each of the stations with identified Vs-profiles except CUSSO and VSAP we acquired and processed available recordings of magnitude 2.5 and greater earthquakes within 3.0 degrees of a given station from EarthScope (www.ds.iris.edu/ds/nodes/dmc/). To avoid any effects of nonlinearity on the empirical spectral ratios, we used only three-component recordings with peak ground accelerations of less than 25 cm/s². For CUSSO and VSAP, we used the same weak-motion dataset that Carpenter et al. (2018) used, consisting of triggered earthquake recordings. Our S-wave HVSR processing involved estimating site fundamental frequencies from preliminary HVSR curves, which we used to establish site-specific processing parameters for the final S-wave HVSRs. We evaluated resonance frequencies and spectral ratio values measured on the final curves.

Preliminary and final S-wave HVSRs were processed following nearly the same steps and as follows: 1. Calculate S-wave and noise window lengths, which, as discussed below, differed between the preliminary and the final processing. 2. Calculate the individual-earthquake amplitude spectra from the windowed, de-trended (i.e., linear trends removed), and tapered three-component S-wave time series after removal of the instrument responses. 3. Calculate pre-P-wave noise amplitude spectra using identical processing parameters and procedure. 4. Smooth the S-wave and noise amplitude spectra; the smoothing approaches differed between the preliminary and final processing as expounded below. 5. Determine reliable S-wave spectral amplitudes at each frequency via signal-to-noise ratios (SNR), where SNR is a function of frequency and is calculated by division of the signal spectrum by the noise spectrum; spectral components with SNR < 2.0 were rejected (following Carpenter et al., 2020). 6. Calculate individual-event HVSRs as the ratios of the horizontal-to-vertical S-wave amplitude spectra but only at frequencies where each component’s SNR is ≥ 2.0. 7. Form the site’s S-wave HVSR as the arithmetic mean from a minimum of three but up to 50 individual-event spectral ratios at each frequency, following the procedure described in Carpenter et al. (2020). 8. Measure resonance frequencies and spectral ratios from the mean HVSR peaks that have prominences of at least 1.0. In other words, HVSR peaks were ignored if their heights above either of their bounding troughs were less than 1.0.

To estimate sites’ fundamental resonance frequencies for final HVSR processing, we attempted to determine preliminary S-wave HVSRs at each site. As expounded in the Supplementary Material to this article, we used fixed-length 60 s S-wave windows to process the preliminary HVSRs, where individual-component S-wave amplitude spectra were smoothed using a Konno-Ohmachi smoother (Konno and Ohmachi, 1998). At some sites, no clear peaks were observed in the preliminary HVSRs. At some other sites in or near sedimentary basins, peaks were observed at unreasonably low frequencies when compared with the predicted 1D responses. Such low-frequency peaks have been attributed to S-wave resonances between the surface and rock layers well below the bases of the Vs profiles determined at these or nearby stations (Mendoza and Hartzell, 2019; Carpenter et al., 2020). The presence of such peaks rendered identifying the first peaks due to the characterized near-surface layers as ambiguous. Figure 4 shows examples of these types of preliminary HVSR curves from seismic stations in Oklahoma, United States. The HVSR curve at OK009 has an unambiguous peak that can be associated with the characterized near-surface Vs structure whereas the spectral ratio at OK032 has two clear peaks. The first peak occurs at a frequency (∼1 Hz) and is much lower than the f _0,1D predicted by 1D modeling (∼5.6 Hz), suggesting site response is affected by deeper geological layers at this site. The HVSR curve for OK001 has no peaks with significant prominence.

To process the final HVSRs, we used site-specific window lengths of at least 10/f _0,HV seconds (similar to the SESAME standards (SESAME project, 2004)), where f _0,HV is the frequency of the first peak on the preliminary site HVSRs, to capture any S-wave reverberations in the near-surface layers. For sites with no clear peaks or with ambiguous first peaks, we assigned a default S-wave window of 10.0 s under the assumption that near-surface, characterized layers would not significantly amplify SH-waves at frequencies lower than 1.0 Hz. The selection of 1.0 Hz is consistent with the peaks seen in HVSRs determined for stations in the same region whose HVSRs appear to be unaffected by deeper structure. S-wave spectra for the final HVSRs were smoothed with running Hanning windows of site-specific lengths f _0,HV /2 Hz. For sites processed with the default window length of 10.0 s, we used Hanning windows of length 0.5 Hz.

The fundamental-mode and peak frequencies, f _p,HV , and corresponding spectral ratios, A _0,HV and A _p,HV , respectively, were measured from the site S-wave HVSR curves by selecting the first and highest peaks with prominences of at least 1.0. As discussed below, although we conducted HVSR analyses at all stations possible, including those whose preliminary HVSR curves lacked clear peaks, we decided ultimately to exclude the HVSR curves with ambiguous first peaks from the final analyses.

The theoretical primary site response parameters, f_0,1D , f_p,1D , A_0,1D , and A_p,1D , were gathered from the 1D responses (Supplementary Table S2), and the proxies were derived from the site Vs profiles (Supplementary Table S2; Figure 3). To explore the parameter space, we first generated scatter plots of each proxy-parameter pair, shown in Figures 6, 7. We evaluated fits to the parameter-proxy pairs, also shown in Figures 6, 7, using linear and two- and three-parameter power laws of the form y = c ∙ x d and y = c ∙ x d + e , respectively, where c, d, and e are free model parameters. The frequency-based parameters and depth-based proxies are in logarithmic (base-10) units for the linear fits.

Site resonances can occur at lower (<1 Hz) and at higher frequencies (≥1 Hz) at deep- and at shallow-soil sites. For example, as listed in Supplementary Table S2, EVIN and VSAP have equal f _0,1D (1.2 Hz) and similar Vs30 of 332 m/s and 289 m/s, respectively, yet these sites have very different Z _1.0 of 36 m and 100 m, respectively and different f _p,1D of 1.2 Hz and 5.1 Hz, respectively. As another example, Z _1.0 at LPAR and MCIL differ by approximately an order of magnitude (657 m versus 20.5 m, respectively) as do their f _0,1D (0.3 Hz and 1.9 Hz, respectively), yet these sites also have nearly equal Vs30 (205 m/s and 218 m/s, respectively) and f _p,1D (1.8 Hz and 1.9 Hz, respectively). Additionally, A_0,1D and A_p,1D , respectively are factors of approximately 2 and 2.5 greater at MCIL than at LPAR. Thus, sites with similar near-surface characterizations may have very different deeper subsurface characteristics that affect site resonance. Therefore, separately regressing subsets of proxy-parameter pairs based on, e.g., depth, frequency, or other, may inadequately characterize the relationships and associated uncertainties. In addition, some recent studies have argued that the near-surface attribute Vs30 is suitable as a site proxy even in settings such as the CEUS where deep geologic layers contribute to site response (e.g., McNamara et al., 2015). To provide a more reliable characterization of the site-response parameter uncertainties that could be encountered across the CEUS, we felt it important to evaluate all proxies using the full range of site-response parameters.

As with the linear scatter plots, we took the base-10 logarithm of the frequency- and depth-based variables prior to fitting the linear models. Table 1 lists goodness-of-fit metrics for the linear and power-law models. The best-fitting linear and power-law models’ parameters are listed in Supplementary Tables S3, S4, respectively. Compared to the three-parameter power-law model, the two-parameter power-law model generally yielded a similar or poorer fit to the parameter-proxy pairs and was not considered for subsequent comparisons and discussion.

Based on the greatest R² values (we used adjusted R², as opposed to non-adjusted R², which accounts for the number of degrees of freedom and allows more meaningful comparisons of the goodness-of-fit between models with differing numbers of free parameters (James et al., 2013)), the linear model is preferred to the power-law models for frequency-depth parameter-proxy pairs, whereas the power-law model is preferred for almost all other pairs. The strongest correlations (i.e., the Pearson’s correlation coefficient) are between depth-based proxies and f _0,1D , with Z _2.5 being the strongest. The same is true for peak frequency, although the correlations are slightly weaker than for fundamental frequency. Apparently, Vs30 also relates to f _0,1D and f _p,1D , with the strength of the former power-law relationship being much stronger than that of the latter.

Amplifications generally have weaker relationships with the various proxies than frequencies, and only IRmax can be considered to have robust and strong relationships (R² ≥ 0.5) with amplification. Although Z _2.5 appears to linearly correlate with A _0,1D , the 95% prediction bounds are quite large (∼±5) at all depths (in this study, the prediction intervals quantify the confidence associated with the fitted curve and depend on the inverse of the Student’s t cumulative distribution function and the variability of a new observation as quantified by the mean squared error). More data may assist with improving the robustness of the Z _2.5 versus amplification relationships. Our analyses indicate that Vs30 does not correlate with amplification in our study area.

We refer to the first and largest peaks in the S-wave HVSR curves as A _0,HV and A _p,HV , respectively, which occur at frequencies of f _0,HV , and f _p,HV , respectively. To avoid the possible influence of deeper structure, such as in or near basins, on the empirical spectral ratios we opted not to include the peaks selected at sites where the first peak was ambiguous. In addition to the examples given in Figure 4, Supplementary Figure S4 presents additional HVSR curves with significant peaks that occur at frequencies well below those predicted by the 1D responses. Although such peaks may represent a lack of Vs-profile resolution in the near surface, they have been observed in Oklahoma by Mendoza and Hartzell (2019), and elsewhere in the CEUS (Yassminh et al., 2019; Carpenter et al., 2020) where they were attributed to deeper sedimentary strata. Additionally, no peaks were observed on some S-wave HVSRs most likely because no significant impedance contrasts underlie the corresponding sites (e.g., OK001 in Figure 4). In total, we measured clear peaks in the S-wave HVSR curves developed at 36 of the 72 stations with sufficient earthquake recordings, which are listed in Supplementary Table S5.

Figure 8 compares the first- and peak-frequencies measured from the S-wave HVSR curves with those determined through the 1D analyses at the 36 seismic stations. In general, the empirical measurements agree well with the first- and peak-frequencies calculated through the 1D analyses, particularly after excluding the frequencies of the first peaks measured at sites suspected to be influenced by deeper structure (Figure 8A). The comparison between f _0,1D and f _0,HV is particularly strong (r=0.90; p=0.00), showing that S-wave HVSR has skill to estimate the fundamental resonance frequencies at these sites, assuming the Vs profiles are reasonably accurate. The comparison of the frequencies of maximum amplifications, f _p , is more scattered (r=0.47; p=0.00) and shows that the calculated peak frequencies tend to be greater than those measured with HVSR. A likely explanation of this is that at frequencies greater than f ₀ , HVSR tends to decrease relative to the site transfer function (Rong et al., 2017; Carpenter et al., 2018; Zhu et al., 2020b). Thus, at sites where f _p,1D > f _0,1D , the first HVSR peaks are more likely to be the largest peaks. In fact, at all except for seven sites, f _p,HV = f _0,HV (Figure 8B), whereas f _p,1D = f _0,1D at only 18 of the same 36 sites. Therefore, the utility of HVSR to estimate the site-response parameters may be limited to the fundamental modes.

Figure 8 also compares the spectral ratios at f ₀ . There is much more scatter in the A ₀ comparison than in the frequency comparisons, much of which could arise from large variabilities observed in HVSR at site resonance frequencies (e.g., Carpenter et al., 2018). Nevertheless, there is a positive, albeit weak correlation (r=0.35; p=0.04) between the empirical and theoretical spectral ratios at f _o . We do not compare A _p,HV and A _p,1D since we observed that f _p,HV and f _p,1D do not correspond at numerous sites (Figure 8B, Supplementary Figure S2).

The results reveal site fundamental- and peak-resonance frequencies can be well predicted from the depths to the 1.0 km/s velocity horizon and from the depth to the maximum impedance ratios in the Vs profiles. The depth of the 2.5 km/s horizon has the strongest relationship with these frequency parameters. These favorable frequency-depth results are consistent with findings in other studies that demonstrated strong correlations between depth-f pairs of variables globally (e.g., Yamanaka et al., 1994) and in the CEUS (e.g., Schleicher and Pratt, 2021; Zhu et al., 2021) and are expected based on Equation 1. Although, in reference to Equation 1, it was also expected that Vs _IRmax would be related to frequency, the regression results indicate that there is no relationship between fundamental or peak frequency and Vs _IRmax . We considered the possibility that the average Vs to other horizons might result in stronger correlations, and thus also evaluated the average Vs to depths of Z _1.0 and to Z _2.5 against f _0,1D and f _p,1D . However, the coefficients of determination were lower than those from the Vs _IRmax regressions indicating that these average-Vs attributes are unreliable as resonance proxies.

Interestingly, although each depth-based proxy has a large value range of more than two orders of magnitude (Figure 3), Vs30 — which characterizes just the upper 30 m—correlates with f _0,1D . The strength of the relationship, which has been observed in other studies in the region (McNamara et al., 2015; Hassani and Atkinson, 2016) and has a relatively large R² of 0.73 in our dataset, which may be due at least in part to most sites having Vs characterized to 30 m depth or less (n = 55). However, the strength of the relationship belies potential problems with reliably predicting f ₀ at low Vs30 (166 m/s to ∼400 m/s) values, where amplification effects may be of greatest concern.

As Figure 9 shows, the f _0,1D -Vs30 residuals are relatively small at low-Vs30 (Figure 9B), but as shown in Figure 9C, the percentage of f _0,1D that the residuals constitute is the very large at Vs30 less than ∼400 m/s. Thus, predicting resonance frequencies using Vs30 at low-Vs30 sites is highly uncertain. Also, as Figure 9 shows Vs30 alone does not distinguish shallow soil from deep soil sites, where Vs30 is most often low and wave propagation effects, particularly at f ₀ , differ greatly from shallow sites (e.g., Hashash et al., 2008). As Figure 9 shows, all but one of the 15 sites with thicker-sediments (depth to the base of the Vs-profile layers ≥ 50 m) have low Vs30 and these same sites have the largest percent residuals (up to 1,100% of f _0,1D ).

Plots analogous to those in Figure 9 for linear f _0,1D –Z _1.0 and f _0,1D –Z _IRmax regressions are given in Supplementary Figure S5. Notably, not only are the residuals for these regressions lower overall than those derived from Vs30, but the residuals are much smaller in terms of percent f _0,1D , indicating that both Z _1.0 and Z _IRmax are much more appropriate for predicting fundamental resonance frequencies than Vs30.

As listed in Table 1, strong (R² ≥ 0.5) relationships between the same proxies—Vs30, Z _1.0 , and Z _IRmax —were found for f _p , also with power (Vs30 only) and linear relationships. Also, we observed the same residual trends for f _p versus Vs30 compared with f _p versus the depth-based proxies: chiefly that peak-frequencies at deep-soil sites are uncertain for the Vs30 relationship and the depth-proxies yield improved residuals overall relative to Vs30. However, the goodness-of-fit metrics for all relationships with f _p were less than those for f ₀ (Table 1).

Although the relationships between the frequency parameters and Z _1.0 (f ₀ and f _p ) and Z _2.5 (f ₀ only) are strong, the depth to the maximum impedance ratio is as well strong but also has an important advantage over the other two: Z _IRmax can be determined at any site with a measured Vs-profile, whereas, Z _1.0 and Z _2.5 require estimations of the depths to those Vs horizons that may be unknown or can be estimated only with large uncertainties.

We found that only one parameter, the maximum impedance ratio, has a strong relationship (R² ≥ 0.5) with fundamental and peak amplifications. This is expected when considering the expression for the 1D amplification at f ₀ for a single layer over a half-space given in Eq. 2: Equation 2 predicts that A _0,1D would likely relate strongly to IRmax, particularly for sites underlain by a single strong impedance ratio. Our 1D calculations include viscoelastic effects and thus, as expected (because most sites are underlain by a single strong impedance contrast), the power-law relationship fits the A _0,1D —IRmax ordered pairs better than the linear model (Table 1). The same is true for the A _p,1D —IRmax ordered pairs, however the power-law fit for this pair is of slightly lower quality.

Carpenter et al. (2020) identified situations where Eqs. 1, 2 provide poorer approximations of fundamental-mode resonance frequencies and amplifications calculated by 1D analyses. They found that f ∼ 0 and A ∼ 0 are respectively less consistent with f _0,1D and A _0,1D at sites underlain by more than one strong impedance contrast. Wang and Carpenter (2023) showed that the consistency further decreases for sites whose Vs-profiles are best characterized as gradient-like, i.e., lacking strong impedance ratios. To expand the findings in those studies and to explore further the effects that the number and magnitudes of underlying impedance ratios have on the calculated primary site response parameters, we categorized each Vs profile based on the number of strong impedance ratios, i.e., 3.0 or greater, although other ratios could be evaluated. We assigned a Vs-type of 1L, ML, UL, G, or B to each profile: 1L–single layer, a single strong impedance ratio; ML–multilayer, two or more strong ratios; UL–upper layer, a single strong ratio, but which occurs at significantly shallower depth than the base of the Vs-profile, namely, less than half the depth; G–gradient, multiple layers in the Vs profile but no strong impedance ratios; and B–bedrock, no strong ratios and all layers’ Vs of at least 760 m/s. Figure 10 shows examples of each Vs-type and gives the distribution of assigned Vs-types. Because there are so few UL (2) and B (1) sites, we only investigated regressions using site response parameters and proxies at 1L, ML, and G sites.

We repeated the linear and power-law model regressions for the 1L, ML, and G Vs-types and present the models’ parameters in Supplementary Tables S6—S11. Comparisons of the linear versus power-law fits for the Vs-type subsets are given in Supplementary Table S12 and a summary of the best-fitting models’ parameters is given in Table 2. The results reveal that analyzing subsets of the dataset based on Vs-type may improve the regressions as quantified by both increasing R² values and decreasing the standard errors. Notably, all regressions are strengthened for single-layer Vs-types. About half of the regressions are strengthened for multilayer sites. For gradient sites, only the single A _p,1D - IRmax pair has a strong relationship.

In Figure 11, we show ordered pair subsets by Vs-type and corresponding best fit relationships for the IRmax-based proxies and each site response parameter. The trends of the fits are consistent regardless of Vs-type, but the ordered pairs are more scattered for multilayer and gradient sites than for single-layer sites. This figure also illustrates that the largest amplifications predicted for this dataset are produced at sites underlain by multiple strong impedance contrasts and that the longest-period resonances occur at single-layer sites.

Because many of the regressions improved when calculated on data subsets based on Vs-type, we reevaluated the parameters picked on the S-wave HVSR curves for the various types. Because some site characterization surveys were conducted at large distances from the seismic stations, and because Vs-structure can vary rapidly over short distances (e.g., Thompson et al., 2009; Hallal and Cox, 2021), we first attempted to identify sites with accurate average Vs profiles. We applied a strategy similar to that used by Pilz and Cotton (2019), who identified KiK-Net sites with reasonable Vs profiles as those with correlation coefficients between empirical and theoretical borehole transfer functions at the fundamental mode of 0.5 and greater. Carpenter et al. (2018) observed in our study area that S-wave HVSR curves are similar to empirical and theoretical site response functions at the fundamental frequencies at sites with accurate Vs profiles. Under the assumption that this observation from Carpenter et al. (2018) is applicable to each of the 36 sites for which we developed S-wave HVSR curves, we calculated Pearson (r) and Spearman (ρ) correlation coefficients for the HVSR curves against the corresponding 1D responses in site-specific frequency bands from 0.5*f _0,HV to 3*f _0,HV . By rejecting sites where either r < 0.5 or ρ < 0.5, we gained confidence that the sites we analyzed have reasonable Vs profiles on average and thus the extracted frequencies and spectral ratios are suitable for comparison against the proxies at these sites. We found only 19 sites that were suitable for such analyses: 15 single-layer, 2 multilayer, 1 upper-layer, and 1 gradient.

Figure 12 compares the fundamental-mode site-response parameters from HVSR and 1D. It is no surprise that the f ₀ parameters are consistent since we required that the theoretical and empirical spectral ratios correlate well around f _0,HV . Therefore, drawing frequency-based conclusions requires additional, independent analyses to confirm the reliability of the site Vs profiles. For this reason, and because f _0,HV has been shown to reliably estimate f _0,1D in numerous studies, as discussed in the Introduction, our HVSR discussion focuses on A ₀ . The scatter in the A ₀ ordered pairs is reduced compared to all the pairs (Figure 8) and there is an improvement in the correlation (r=0.52; p=0.02). These results suggest that S-wave HVSR curves can provide a first-order estimation of A _0,1D and the best-fitting linear relationship is A 0,1 D = 0.4 ∙ A 0 , H V + 3.1 . (5)

Figure 13 presents example S-wave HVSR curves for each Vs type at sites where Vs is well characterized. The close match between the empirical and theoretical curves suggests that S-wave HVSR reflects characteristics of the 1D site response functions. The only example of a bedrock site (station WMOK) was included even though it was not compared in Figure 13 (since f _0,HV could not be determined for this station).

The examples shown in Figure 13 also suggest that S-wave HVSR curves differ by site type. The HVSR curve at SIUC (1L site) has regularly spaced peaks—i.e., the peaks occur at approximately odd multiples of the fundamental frequency, as predicted by Equation 1—that diminish with increasing frequency, as predicted for a single layer over a half-space. The HVSR curve at U43A (ML site) is more complex, with irregularly spaced peaks—i.e., the HVSR peaks are not separated by an approximately constant frequency interval—and with a relatively high peak (≥5) indicating the presence of at least one larger impedance ratio beneath the station. For OK033, the UL example, only a single clear HVSR peak of value 4.0 is apparent, which occurs at a relatively high frequency of 11 Hz. Peaks at higher frequencies are expected at OK033 but may not have been observed due to data sampling limitations. The HVSR curves, as well as the site transfer functions, at UL sites may resemble those at 1L or ML sites. The largest peak at the G site V35A is much less than the peaks observed at the sites underlain by at least one strong impedance ratio, corroborating that the site lacks a strong underlying impedance contrast. Unfortunately, the Nyquist frequency limits evaluation of the HVSR at frequencies greater than 20 Hz at V35A, but comparing curves at another G site at f >> f _0,HV and with a higher Nyquist frequency, OK009 (Supplementary Figure S4; OK009 was categorized as a basin site and the f _0,HV had to be excluded from the discussed analyses), indicates that the empirical-theoretical spectral-ratio curves are indeed comparable for this site type. For the B site, WMOK, the HVSR curve is flat, as expected for a site at which little-to-no site response is predicted.

Hashash et al. (2014) recommend a relatively fast velocity as reference site conditions for ground motion models and for predicting site amplifications for the CEUS, in the range 2,700 to 3,300 m/s. However, only nine of the Vs profiles we collected include velocities in this range or faster (the median of the sites’ maximum Vs is 1,524 m/s and 1,134 m/s and 2,129 m/s bound the second and third quartiles of the maximum Vs). Supplementary Figure S3 demonstrates that the amplifications at these nine sites lie in the upper ranges of those determined for all sites. Thus, the distribution of amplifications we determined may not be representative of CEUS amplifications that are referenced to the 2,700 to 3,300 m/s range. Nevertheless, we emphasize that site resonances do not depend on absolute Vs maxima, but rather on impedance ratios as indicated by Eq. 2 and suggested by the amplification versus IRmax relationships, and thus note that even the recommendations from Hashash et al. (2014) may be insufficient to account for site resonance in the CEUS.

Another potential limitation of our study results is related to the IRmax-based relationships we developed. The related concern is that the Vs profiles used in this study may not have reached the depths to all impedance contrasts that contribute to site resonance in the frequency range of engineering interest. Thus, it is possible that either or both IRmax and Z _IRmax are too small at some sites, in which case Vs _IRmax may need to be recalculated. As observed in this study, S-wave HVSR can assist with identifying sites that have deeper unmodeled impedance contrasts of relevance (Figure 13, Supplementary Figure S4). And based on the consistencies between the empirical and theoretical spectral ratios in Supplementary Figure S2, it appears that the Vs profiles at sites with underlying thick sediment deposits in the upper Mississippi Embayment and Atlantic Coastal Plain extend to adequate depths to account for all relevant impedance contrasts. Outside of these domains and as discussed in Sections 4.2 and 5.4, the spectral ratio comparisons (Supplementary Figures S2, S4) suggest that some sites lack the detailed structural models needed to assess important resonances from deeper layers.

Because some sites appear to be subject to lower-frequency, unmodeled resonances, our results will be corroborated or strengthened once additional Vs profiles that account for all relevant impedance ratios become available. Nevertheless, we anticipate that the relationships developed using IRmax-based proxies in our dataset will be useful in the CEUS and perhaps elsewhere, as discussed in what follows.

Relevant to the A _0,1D and A _p,1D versus IRmax relationships, Z _IRmax and IRmax are uncorrelated. Also, the distribution of IRmax used in the regressions is similar to the distribution for high-Vs-max subsets of the dataset (Supplementary Figure S6). Thus, we do not expect that very different relationships for A _0,1D and A _p,1D versus IRmax would be developed if deeper depths and higher Vs (e.g., in the range of the Hashash et al. (2014) recommendation) were reported for all sites.

Related to f _0,1D and f _p,1D versus Z _IRmax , at the 73 sites where Z _1.0 was determined, the maximum impedance contrast occurred at that depth at most sites (53 sites; Supplementary Figure S7) or at shallower depths (14 sites). At just six sites, Z _IRmax was found to be deeper than Z _1.0 . Also, at the 14 sites where Z _2.5 was measured, Z _IRmax corresponds with Z _2.5 at all but three sites. Therefore, it is possible that stronger impedance contrasts exist at depths corresponding to Z _1.0 , Z _2.5 , or deeper, at sites where these depths were not determined.

To assess the extent to which our results involving IRmax and Z _IRmax may differ from those developed at adequately characterized sites, we recomputed the regressions using the best-fitting model types listed in Table 2 for “All” Vs profiles (linear versus power-law) on a subset of 39 sites where S-wave HVSRs were determined at frequencies down to 0.1 Hz (Supplementary Figure S2) and which lacked ambiguous low-frequency peaks. The resultant regressions for this subset (Supplementary Figure S8)—all of which have R² ≥ 0.5—nearly coincide with the relationships derived from the entire dataset. The consistencies of these relationships suggest that the results involving the IRmax and Z _IRmax proxies at all sites are reliable. Nevertheless, it would be beneficial to reassess the strengths of these IRmax-based relationships as additional Vs profiles that account for all relevant impedance ratios become available.

Another potential limitation of our study is that we only considered single proxies to predict single site-resonance parameters. It is possible that multiple regression analyses involving two or more of the proxies we evaluated, and perhaps others we did not consider, would yield improved relationships with the primary site response parameters.

The strong correlations revealed in this study support the findings of previous studies that the Z‐based proxies can predict resonance frequencies and thus also support the use of Z_1.0 and Z_2.5 as site terms (e.g., Day et al., 2008) in CEUS ground-motion models. The results also demonstrate that ZIRmax may serve the same purpose for locations where Z_1.0 and Z_2.5 are unavailable, although S‐wave HVSR observations should be used to identify unknown impedance contrasts that produce unmodeled low‐frequency resonances of engineering concern. IRmax correlates strongly with fundamental‐mode and peak amplifications, and thus serves as another site term for CEUS ground‐motion models. As with Z_IRmax , S‐wave HVSRs can be used to ensure that additional large impedance ratios are not present beneath the base of site Vs profiles. The Z‐based proxies may be useful in developing design spectra for linear ground motions across appropriate site-specific periods (i.e., the inverse of frequency) since fundamental frequencies determined from weak-motion response spectral ratios and those determined from 1D linear response analyses have been shown to be consistent (Wang and Carpenter, 2023). Also, because the Z‐based proxies correlate strongly with resonance frequencies and IRmax correlates strongly with amplification, pairs of Z‐based proxies and IRmax can be used to modify predicted bedrock spectral amplitudes at a given site, using the relationships we developed, to include linear site effects. Table 2 indicates that, when possible, relationships developed for “1L” sites should be used for these purposes at sites with simple Vs structures and “ML” should be used at sites with multiple underlying strong impedance ratios to predict f_0,1D using Z_1.0 or Z_2.5 , and for f_p,1D using Vs30, but only at higher Vs sites.