Common procedures used to account for spatial variability of shear wave velocity (V s) in one-dimensional (1D) ground response analyses (GRAs) , such as stochastic randomization of V s or increasing small -strain damping, have been shown to improve seismic site response predictions relative to 1D GRAs where no attempts are made to account for spatial variability. However , even after attempting to account for spatial variability using common procedures, 1D GRAs often still yield results that are different than ground motions recorded at many downhole array sites. When 1D predictions differ from observations, the site is typically considered to be too spatially variable to effectively use 1D GRAs. While there is no doubt that some sites are indeed too variable for 1D GRAs , it is also possible that simple 1D analyses could still be effectively used at many sites if spatial variability is accounted for via a more rational, site -specific approach . In this study, an H/V geostatistical approach for building pseudo- 3D Vs models is implemented to account for spatial variability in 1D GRAs . The geostatistical approach is used to generate a uniform grid of Vs profiles that have been scaled to match fundamental site frequency estimates from horizontal -to-vertical spectral ratio (H/V) noise measurements . In this paper, 1D GRAs are performed for each grid -point and the results are statistically combined to reflect the average site response and its variability. This 1D application is demonstrated at the Treasure Island and D elaney Park Downhole Array sites, where it is shown to produce superior fit s to the small- strain recorded site response relative to existing approaches used to account for spatial variability in 1D GRAs . Using the proposed approach, we also investigate the lateral area that is likely influencing site response at each site and show that it could extend to significan t distances (as much as 1 km) from the boreholes . 2 Introduct ion The site -specific influence of subsurface geo -materials on the intensity and frequency content of earthquake ground motions, commonly known as site effects, is a crucial factor that controls the seismic hazard for many engineering project s. For this purpose, numerical wave propagation simulations , also known as ground response analyses (GRAs) , are commonly performed to capture the effects of the local site conditions on ground shaking. O ne-, two -, and three- dimensional (1D, 2D, and 3D, respectively) GRAs have been developed to simulate t he complexity of the seismic waves’ propagation patterns , while linear, equivalent -linear , and fully nonlinear methods have been used to model the soil behavior (Zalachoris and Rathje, 2015) . In practice, however, the vast majority of GRAs are performed using only 1D simulation of vertically propagating, horizontally polarized shear (SH) waves . The 1D GRA assumes a medium composed of homogenous , horizontally stacked layers of infinite lateral extent subject to a planar wavefront of SH waves. In reality , the subsurface conditions at many sites are much more complex , with heterogeneities of various length scales in different directions (de la Torre et al., 2019; Kramer, 1996) . These complexities engendered the need to utilize direct observations of site ef fects from earthquake ground motions recorded a t multiple depths in downhole arrays to better understand true site response . While several studies have shown that 1D GRAs can predict key features of the recorded response at some downhole array sites (Kim and Hashash, 2013; Tao and Rathje, 2019) , all researchers agree that the predicted and recorded amplification amplitudes rarely match (Afshari and Stewart, 2019; Boore, 2004; Kaklamanos et al., 2013; Thompson et al., 2009) . These discrepancies have been attributed to multiple sources of potential uncertainties , including: (1) inaccurate soil properties, (2) unrealistic constitutive model s, and (3) violations of the 1D site assumptions (Rathje et al., 2010; Thompson et al., 2012) . The latter includes non- vertical wave incidence, spatial variability of the subsurface conditions, and topographic and basin effect s, among others (Kaklamanos et al., 2013) . Many research efforts aimed at addressin g the uncertainties affecting 1D site response predictions have been undertaken. Arguably, of the potential sources of uncertaint y, inaccurate soil properties are the easiest to identify and mitigate through reliable measurements with the application of quality control measures . Additionally, significant research has been performed to determine constitutive models that best capture nonlinear soil behavior, although this remain s an active area of research (Groholski et al., 2016; Phillips and Hashash, 2009; Thompson et al., 2012) . Several studies have also investigated breakdowns in 1D GRA assumptions, specifically the spatial variability of subsurface conditions , through both observations at downhole array sites (Assimaki and Steidl, 2007; Rathje et al., 2010; Teague et al., 2018; Thompson et al., 2009, 2012) and numerical simulations (de la Torre et al., 2019; Field and Jacob, 1993; Makra and Chávez -García, 2016) . For example, Fie ld and Jacob (1993) , through Monte Carlo simulations of the shear wave velocity (Vs) and damping profiles, concluded that incorporating variability in the input soil properties significantly affects the predicted amplifications. Assimaki and Steidl ( 2007) , in a study at five downhole array sites, observed that the back- calculated damping values are substantially higher than those estimated from lab measurements, which was attributed to wave scattering in the field. Thompson et al. ( 2009) used 3D spatially correlated random fields to predict the site response 3 at a downhole array from the Kiban -Kyoshin network (KiK -net) in Japan. They found a more satisfactory fit to the observed response by incorporating subsurface variability and argued that this variability diminishes the downgoing wave effect. Rathje et al. ( 2010) assessed the influence of Vs and nonlinear soil properties through Monte Carlo simulations and found that modeling their variability generally reduces the predicted surface motions. Makra and Chávez- García (2016) showed through a set of 2D and 3D simulations of the Euroseistest network that lateral heterogeneity may cause amplifications at frequencies different from the resonance frequencies predicted from 1D V s profiles. Researchers have thus hypothesized that subsurface heterogeneity causes scattering and attenuation of the seismic waves, which lead to an apparent loss of energy that is not capt ured by 1D GRAs (Afshari and Stewart, 2019; Tao and Rathje, 2019; Thom pson et al., 2009, 2012; Zalachoris and Rathje, 2015) . Recognizing the indisputable effect of variability in soil properties, several methods with varying degree s of complexity have therefore been proposed in recent years to better model the variability of subsurface conditions. For example, stochastic randomization of Vs (Rathje et al., 2010; Toro, 1995) , utilizing suites of V s profiles from surface wave testing (Griffiths et al., 2016a, 2016b; Teague et al., 2018; Teague and Cox, 2016) , and increas ing the small- strain damping ratio (Dmin) through damping multiplicative factors or other means (Afshari and Stewart, 2019; Tao and Rathje, 2019) are all simplified methods that attempt to incorporate complex 3D variability of subsurface conditions through modifications to 1D analyses . The blind application of randomized Vs profiles to account for aleatory variability has been recently called into question by a number of studies (Griffiths et al., 2016a, 2016b; Teague et al., 2018; Teague and Cox, 2016) , which have shown that Vs profiles developed in this manner often do not represent realistic site conditions . While a new Vs randomization approach proposed by Passeri et al. ( 2020) looks to be more promising, it has not yet been widely used in practice. Teague et al. ( 2018) found that utilizing suites of Vs profiles from the inversion of surface wave data was more successful at predicting the observed site response than stochastic randomization at the Garner Valley Downhole Array. Tao and Rathje (2019) , in a study on four downhole array sites, found that either increas ing D min or using Vs randomization provided comparable results when site-specific inputs that reflect the degree of spatial variability are used. Afshari and Stewart (2019) , by investigating three different damping models at 21 California downhole arrays, found that the damping model inferred from the site -specific decay parameter (κ 0) generally provide d better results than laboratory- and V s- based D min values . Notwithstanding potential improvements gained by attempting to incorporate some degree of spatial variability into 1D GRAs , adopti ng these methods in mainstream p ractice is still frustrated by a number of limitations. For example, Vs randomization procedure s require site -specific input parameters for stochastic modeling that are not typically available, necessitating the assumption of generic parameters derived for other sites (Stewart et al., 2014; Toro, 1995) , bringing additional uncertainty into the site response predictions (Passeri et al., 2019; Teague et al., 2018) . Additionally, incorporating non-unique V s profiles from the inversion of surface wave data requires significant expertise and use of arrays of great apertures to infer information at depths governing site response . However, no guidance exists on the footprint that needs to be covered by surface wave arrays in order to appropriately capture spatial variability , which might lead to the inclusion of subsurface conditions that are not actually influencing the site response , or vice versa. Also, while damping multiplicative factors are believed to be dependent on the level of spatial 4 variability , the appropriate factor cannot be known a priori and no current method is available to link information that can be measured at the site to an appropriate factor. Furthermore , although the application of D min multiplicative factors reduces the predicted a mplifications near the fundamental site frequency, this also tends to result in underpredicted amplifications at higher frequencies. For reason s noted above , it is clear that incorporating site-specific spatial variability into GRAs is important, especially using approaches that are not fully stochastic or require inputs that cannot be known a priori . Furthermore, investigating an approach that could be used within the 1D framework is desir able, as 1D GRAs are easier to perform and within the proverbial comfort zone of most engineers. In a companion paper (Hallal and Cox, 2020) , a new framework called the ‘H/V geostatistical approach’ was proposed for developing pseudo- 3D Vs model s that capture key aspects of site- specific spatial variability using a single measured Vs profile and estimates of fundamental site frequency (f 0; a key parameter governing site effects) from horizontal -to-vertical spectral ratio (H/V) noise measurements (f0,H/V). This paper incorporat es this 3D site -specific spatial variability in a meaningful and practical way into 1D GRAs . Besides incorporating a realistic pseudo -3D Vs model, a key advantage of implementing the H/V geostatistical approach in 1D GRAs is the ability to predic t the site response as a function of the spatial variability across different footprints , and consequently, provide insights on the spatial area likely influencing site response. First, information on the case study sites, the Treasure Island Downhole Array (TIDA) and Delaney Park Downhole Array (DPDA) , and the ir recorded ground motions are presented . Second , the theoretical linear -viscoelastic transfer functions obtained from a single, invasive V s profile at each site are compared with the empirical transfer functions calculated from small- strain ground motions recorded at the downhole arrays . Third, the spatial variability of f 0 across each site measured using widely distributed, non- invasive H/V measurements is quantified and the discrepancies between the empirical and theoretical transfer functions are re- examined in light of the site -specific spatial variability revealed by the H/V measurements . Fourth, w e describe the application of the H/V geostatistical approach to 1D GRAs and demonstrate it by present ing its results for the two case study sites. Finally, insights are provided on the potential spatial area affecting site response and recommendations are given for future improvements . Site and Ground Motion Information The t wo downhole array sites for w hich pseudo- 3D Vs models were developed in the companion paper (Hallal and Cox, 2020) , the TIDA and DPDA, are investigated in this study regarding their recorded versus predicted linear -viscoelastic site response . Importantly, these two sites are distinct in that TIDA is believed to be an archetype of 1D conditions , where 1D GRAs are expected to provide a satisfactory match to empirical/recorded ground motions , whereas DPDA is the ant ithesis of such conditions , where even highly modified 1D GRAs fail to accurately represent empirical ground motions (Tao and Rathje, 2019) . For each site, the geology and subsurface properties are briefly described in what follows, along with the downhole array and ground motion information. A more detailed discussion on the geology and subsurface conditions is presented in the companion paper . 5 Treasure Island Downhole Array The T IDA site is located in northern California on an artificial hydraulic fill island constructed on a natural sandspit northwest of Yerba Buena Island (YBI ), an outcropping bedrock island in San Francisco Bay (Rollins et al., 1994) . A Vs profile measured using PS suspension logging (Graizer and Shakal, 2004) and simplified subsurface stratigraphy for TIDA are shown in Figure 1a. Existing data indicates that the subsurface near the downhole array site consists of the following sequence of layers: loose sandy hydraulic fill, Young Bay M ud (YBM), and Old Bay M ud (OBM) overlying bedrock. The bedrock is located at approximately 90 m and is composed of interbedded sandstone, siltstone, and shale of the Franciscan Formation (Gibbs et al., 1992; Pass, 1994) . The TIDA site is instrumented with triaxial accelerometers at the surf ace and throughout the soil profile in to the bedrock, as indicated by the black triangle symbols in Figure 1a. Only the surface and deepest accelerometers, which are numbered 00 and 06 , respectively, are considered further hereafter. Downhole accelerometer number 06 will simply be referred to as the “rock” accelerometer. A suite of 52 low amplitude ground motions recorded at the site ( both horizontal components from 26 unique events) was selected to calculate the “true”/empirical site response between the r ock and surface accelerometers . By including only low amplitude motions, we ensure that the challenges associated with modeling nonlinear soil behavior do not interfere with attempts to understand spatial variability . The peak ground accelerations (PGA s) of the surface motions were on the order of 0.00 2 to 0.04 g, with Richter Local Magnitude (M L) ranging from 3.0 to 5.6 Figure 1. Measured Vs profiles, triaxial accelerometer sensor depths (black triangle symbols), and simplified stratigraphy at: (a) the Treasure Island Downhole Array (TIDA), and (b) Delaney Park Downhole Array (DPDA) sites. The measured Vs profile at TIDA is from de ep PS suspension logging (Graizer and Shakal, 2004) and the stratigraphy is based on the seismic downhole borehole presented in Gibbs et al. (1992) . The measured Vs profile at DPDA is from seismic downhole (DH) (Thornley et al., 2019) and the stratigraphy is based on geologic data from Combellick (1999) . 0 50010001500200025003000Shear Wave Velocity, Vs (m/s) 0 20 40 60 80 100 120 140Depth (m)TIDA 00 01 02 03 04 05 06FILL YBM ?? OBM ?? ROCK ?? (a)PS (Graizer and Shakal, 2004)0 500 1000 1500 2000Shear Wave Velocity, Vs (m/s) 0 10 20 30 40 50 60 70Depth (m)DPDA D0 D1 D2 D3 D4 D5 D6GOW CLAY ?? TILL ?? (b)DH (Thornley et al., 2019) 6 and distance ranging from 11 to 120 km. The se ground motions were selected from a larger candidate set based on a signal -to-noise ratio (SNR) criterion of greater than 3 dB between 0.5 and 10 Hz. This frequency range encompasses those frequencies that are of most interest from an engineer ing perspective. Additionally, the expected fundamental mode resonance frequency at TIDA is greater than the minimum screening frequency of 0.5 Hz. Delaney Park Downhole Array The D PDA site is located in downtown Anchorage, Alaska in Delaney Park. A Vs pro file measured using seismic downhole (Thornley et al., 2019) and a simplified subsurface stratigraphy for DPDA are shown in Figure 1b. Geologic data indicates that the subsurface near the downhole array site consists of the following sequence of layers: glacial outwash (GOW) , followed by clay from the Bootlegger Cove Formation overlying glacial till (Combellick, 1999) . The glacial till is located at approximately 47 m and consists of very dense sand and gravel . The DPDA site is instrumented with triaxial accelerometers at the surface and through out the soil profile in to the till, as indicated by the black triangle symbols in Figure 1b. Only the surface and deepest accelerometers, which are denoted D0 and D6, respectively, are considered further hereafter. Downhole accelerometer number D6 will simply be referred to as the “rock” accelerometer. A suite of 56 low amplitude ground motions recorded at the site ( both horizontal components from 28 unique events ) with PGAs at the surface on the order of 0.001 to 0.01 g, M L ranging from 3.0 to 5.1, and distance ranging from 10 to 96 km were selected to calculate the empirical site response between the rock and surface accelerometers . The same SNR criterion of more than 3 dB in the 0.5 – 10 Hz frequency range was applied during the selection of these ground motions , noting that the expected fundamental mo de resonance frequency at DPDA is greater than the minimum screening frequency of 0.5 Hz. Linear -Viscoelastic GRAs General Framework The “true” site response at each downhole array site is represented in this study using the empirical transfer function (E TF). The ETF is a frequency -dependent measure of site response meant to represent the SH wave transfer function, and is calculated as the ratio of the Fourier amplitude spectra (FAS) of the horizontal surface accelerations divided by the horizontal rock acceleration s. ETFs between the rock and surface accelerometers were computed using procedures similar to those documented in Teague et al. ( 2018) . The acceleration time histories were baseline- corrected and then filtered using a fifth -order Butterworth acausal filter with a passband of 0.5 – 10 Hz ; the pass band is determined based on frequencies with SNR > 3 dB . The FAS were smoothed using a log- scale rectangular window , as described by Goulet et al. ( 2014) , and the ETF for a single ground motion was then computed as the ratio of the smoothed FAS. The representative ETF for the site is finally taken as the geometric mean of all individual ETFs f rom the two horizontal components of all earthquake events . For a lognormal distribution, which is commonly assumed for the population of ETFs (Thompson et al., 2012) , the geometric mean is equivalent to the median of the sample. We also compute the natural log arithmic standard deviation of the ETFs (σlnETF) at each frequency as a means to capture the variability in the ETFs . As the ETF s are 7 calculated from weak ground motions with PGA s at the surface much less than 0.1 g, which is reported by several researchers as the threshold below which soil nonlinearity is imperceptible (e.g., Kaklamanos et al., 2013) , they are believed to be a representative of the linear -viscoelastic site response. To evaluate the modeling accuracy of GRAs , the theoretical transfer function (TTF) is computed and compared with the median ETF. Linear -viscoelastic, theoretical SH wave transfer functions were computed for the V s profile s using the closed- form solution for damped, uniform soil over elastic bedrock (Kramer, 1996) . The required input soil properties include V s, mass densit y, and D min for each layer. The damping profiles were assigned using the laboratory- based Dmin from Darendeli (2001) . It should be noted that since the ETFs represent the surface -to-rock ratio in which the rock time histor ies include both upgoing and downgoing waves, the TTF is modeled to include those effect s using what is referred to as the within boundary condition at the rock accelerometer (Thompson et al., 2012; Zalachoris and Rathje, 2015) . The presence of downgoing waves in the recorded rock motions w as verified by deconvolving the rock acceleration time history with that of the surface, similar to procedures described by Wen and Kalkan (2017) . The modeling a ccuracy of GRAs is qua ntitatively assessed in this study using two goodness - of-fit parameters , similar to Teague et al. (2018) : the Pearson correlation coefficient (r) and the transfer function misfit (m TF). Several parameters have been used in previous studies to judge the goodness -of-fit between empirical and theoretical site response, all of which have positive a nd negative attributes. The use of multiple parameters is therefore recommended to offset the limitations of a single parameter (Tao and Rathje , 2019; Teague et al., 2018) . The Pearson correlation coefficient quantifies how well the shape of the ETF and TTF align, including the location of resonance frequencies , with r > 0.6 proposed by Thompson et al. ( 2012) as an indicator of good alignment . However, r is oversensitive to extreme amplitude values (Legates and McCabe, 1999) , which can lead to a bias even when the peaks are well -aligned. In addition to quantifying the alignment of the peaks , we quantify the residuals between the TTF and ETF using a transfer function misfit valu e. This value indicates, on average, how many standard deviations the TTF is away fro m the ETF (Teague et al., 2018) . Thompson et al. (2012) computed goodness -of-fit parameters over logarithmically spaced frequencies between t he first (i.e., fundamental mode ) and fourth (i.e, 3rd-higher mode) peaks in the TTF. However, such criterion might fail to capture the fit to the entire fundamental mode of the ETF , especially if the fundamental mode of the ETF occurs at a frequency lower than that of the TTF. Thus, we propose that the frequency range should be determined by the ETF. Also, to reflect the fit to the width of the ETF ’s modes , we propose that the range should extend beyond the frequencies at which the peaks occur. As such, we compute r and m TF at logarithmically spaced frequencies between the lower frequency of the half -amplitude bandwidth of the first peak to the higher frequency of the half -amplitude bandwidth of the four th peak in the median ETF. Such calculation is only possible if the ETF has clear peaks. In this study, w e limit the higher frequency to 10 Hz regardless of the location of the fourth peak in the ETF, as the ground motions only have SNR greater than 3 dB u p to that frequency. Assessment of 1D GRA s at Case Study Sites The TTF s obtained from the measured V s profiles at the TIDA and DPDA site s are compared to the ir respective median ETF ± σ lnETF in Figure 2. While the locations of fundamental and higher 8 mode resonan ce frequencies in the TTF well -match those of the ETF at both sites , the amplitudes of the TTFs are significantly over -estimated and the widths of the TTFs at the fundamental resonance frequency are generally underestimated (particularly at TIDA) . The discrepancies in the transfer function amplitudes are more pronounced at the DPDA site (Figure 2b). Additionally, the ETF at the TIDA site shows the presence of a double amplification peak near the fundamental mode resonan ce frequency, as labeled in Figure 2a. However, this secondary peak is not captured by the TTF. The values of r and m TF at both sites are indicated in brackets in the figure legend. Based on the r > 0.6 criterion set by Thompson et al. ( 2012) , the TIDA site is considered to be a well-modeled/ good- fit site whereas DPDA falls short of such a classification. While the mismatch near the fourth peak at the DPDA site might be contributing to this lower r value, it is likely that the more pronounced overprediction of transfer function amplitudes is driving it. This is more evident in the misfit values, which show that the TTF at the DPDA site is on average 2.34 standard deviations away from the ETF , compared to 1.63 standard deviations at the TIDA site. As noted in the Introduction, s everal researchers have observed s imilar discrepancies between the amplitudes of linear -viscoelastic ETFs and TTF s at downhole array sites and attributed them to inaccurate soil properties and/or violation of the 1D GRA assumptions (Kaklamanos et al., 2013; Rathje et al., 2010; Tao and Rathje, 2019; Teague et al., 2018; Thompson et al., 2009, 2012) . Specifically, research suggests that the subsurface spatial variability may cause attenuation of seismic waves through mechanisms not captured in laboratory damping measurements, amplification at frequencies different from the expected resonance frequencies, and redistribution of the seismic energy over a wider frequency band (de la Torre et al., 2019; Makra and Chávez - García, 2016) . Because Figure 2 shows similar observations , we believe that subsurface variability, or more generally breakdowns in the 1D GRA assumptions, are a highly plausible reason for the mismatch. Similar to Tao and Rathje ( 2019) , we also hypothesize that subsurface conditions at the DPDA site are more variable than those at the TIDA site, causing more scatter ing in the seismic Figure 2. Comparison of the empirical and theoretical transfer functions (ETF and TTF, respectively) at the: (a) TIDA, and (b) DPDA sites. The Pearson correlati on coefficient (r) and transfer function misfit (m TF) values between the experimental and theoretical data at each site are indicated in brackets in the legend. 0.5 1 3 10 Frequency (Hz)100101102TF AmplitudeTIDA Secondary Peak (a)Med. ETFlnETF TTF PS [r=0.63 ; mTF=1.63] 0.5 1 3 10 Frequency (Hz)DPDA (b)Med. ETFlnETF TTF DH [r=0.54 ; mTF=2.34] 9 waves and consequently poorer model predictions. To test these hypotheses, in- situ, surface -based measurements aimed at quantifying the variability in subsurface conditions around the downhole array sites were conducted. H/V Measurements To efficiently characterize the spatial variability over a wide area around the downhole array sites, surface -based H/V ambient noise measurements were conducted. Background information on the H/V method along with the framework implemented in this study t o process the H/V measurements are discussed in the Hallal and Cox ( 2020) compani on paper . The variation in f 0,H/V obtained from H/V measurements at 97 and 96 locations around the TIDA and DPDA site s, respectively, are shown in Figure 3. The f 0,H/V values at TIDA ( Figure 3a) are quite consistent across most of the island with values of approximately 0.75 – 0.80 Hz. However, they increase abruptly in the south where bedrock shallows near the causeway to YBI. Unlike TIDA, the distribution of the f 0,H/V values at DPDA ( Figure 3b) shows considerable heterogeneity in almost all lateral direc tions. The heterogeneity is minimal in close vicinity to the downhole array, but becomes more significant in as little as 300 m away. The f 0,H/V values decrease in the north direction but increase in the south- east quadrant. Closer scrutiny of the mapped spatial variability in f 0,H/V values could help us test the relation between the complexity of the subsurface conditions and the mismatch in the 1D GRA predictions . In order to quantitatively represent the extent of the spatia l variability at each site, we choose to use the semivariogram of f 0,H/V, which is a measure of the degree of variability between observations as a function of their separation distance. We normalize the measurements by the maximum f0,H/V at each site to ensure that the difference in magnitude of f 0,H/V between the sites does not affect the semivariance. We also compute the omnidirectional semi variogram since we are interested in the overall variability rather than that along a certain azimuth. Figure 3. Fundamental site frequency map inferred from H/V measurements (f 0,H/V) across the: (a) TIDA, and (b) DPDA sites. The location of the downhole array at each site is denoted with a black triangle. 10 The semivariograms of the normalized f0,H/V values at both sites are compared in Figure 4. While t he sites have comparable variability at small distances, DPDA has an increased semivariance beyond a lag distance of 200 m. While Tao and Rathje ( 2019) suggested that the DPDA site is more variable than TIDA, this was based on a qualitative assessment of the different geology and depositional environments. The authors highlighte d the need for further investigation through detailed site characterizations. Our results corroborate their hypothesis that DPDA has more variability in the subsurface properties , which is likely responsible for the poorer 1D TTF results shown in Figure 2b. Thus , the f0,H/V semivariograms ( Figure 4) provide quantitative evidence that supports the presence of more complex subsurface conditions at the DPDA site than at the TIDA site . This is also in line with findings presented in the companion paper (refer to Figure 11a in Hallal and Cox 2020) , which show that the standard deviation in f0,H/V is higher at DPDA beyond a specific area around the downhole arrays. Additionall y, if the degree of spatial variability determines the ext ent of wave scattering, then Figure 4 suggests that the variability at distances relatively far from the site may be of interest. Not only do the mapped f0,H/V values show that the DPDA site is more variable, but they also potentially explain the observed double amplification peak in the ETF near the fundamental mode resonan ce frequency at the TIDA site . The observed secondary peak in the ETF at a frequency just higher than the fundamental mode occurs at 0.93 Hz (Figure 2a). Interestingly, this frequency closely coincides with the average f0,H/V near the causeway to YBI (Figure 3a). Thus, it might be possible that the site response at the TIDA site is influenced by the subsurface conditions at great distances, as much as 1 k m away . This possibility is examined in greater detail later in the paper. The concept that the site response at a s pecific location is influenced by some non- finite radial area around the site seem s intuitive . Additionally, it is well known that low frequency earthquake waves have long wavelengths, which sample large portions of the subsurface. However, no experimental study we are aware of has examined the lateral extent that influences site response , or more importantly, provided quantitative evidence indicating that variable subsurface conditions at significant distances could influence the recorded site response. We would like to emphasize that we are not referring to topographic or basin effects. While these effects can also extend over significant distances, we do not believe they are influencing the small -strain ground motions recorded at either downhole array site. Figure 4. Comparison of the experimental omnidirectional semivariograms of the normalized f 0,H/V values from the TIDA and DPDA sites. 0 100 200 300 400 500 600 700 800 Lag Distance, h (m)10-510-410-310-210-1Semivariance, TIDA DPDA 11 As highlighted in Hallal and Cox ( 2020) and this paper, t he mapped f 0,H/V values at each downhole array site underscore the importance of incorporating site-specific spatial variability to captu re true site response. Measurements of site -specific vari ability should reflect both its magnitude and skewness; we refer to skewness as variability being favored in a certain direction, such as having mostly stiffer or softer conditions in a certain area or direction around the site. We have observed that the DPDA site has variability that i s moderately skewed to higher f 0,H/V toward the south- east, whereas the TIDA site has relatively lower variability that is skewed to higher f 0,H/V values toward the south, near the causeway to YBI . The observations suggest that the magnitude of the spatial variability is related to the amplitude of the peaks and the skewness can be related to the amplification at frequencies different from the expected resonance frequencies in the ETFs (Figure 2). Such variability could not be represented by purely stochast ic models based on generic parameters derived for other sites and, consequently, these methods may not be able to account for site-specific spatial variability in a meaning ful way. Thus , to attempt to model true site effects, the H/V geostatistical approach, which account s for the site -specific spatial variability in a rigorous and methodical way , is implemented herein . Implementation of the H/V Geostatistical Approach in 1D GRAs The H/V geostatistical approach proposed in the companion paper (Hallal and Cox, 2020) uses site-specific, simple measurements to develop reasonable and representative pseudo- 3D Vs model s that can be used to incorporate spatial variability into 1D, 2D , and 3D GRAs . In this paper we focus exclusively on incorporating these pseudo- 3D Vs models into 1D GRAs , with the ultimate goal of modeling 3D response in the future, at least for sites where the 1D approach is shown to not be adequate . A schematic summary of the H/V ge ostatistical approach and its implementation in 1D GRAs is shown in Figure 5. First, the irregularly spaced f 0,H/V measurements are kriged to obtain a uniformly estimated map of f 0,H/V values. Second , using the f0,H/V value closest to the measured Vs profi le (f0,H/V,Vs ) and the estimated f0,H/V map, a uniformly estimated scaling factor (S f0) map is obtained. Third, a spatially variable pseu do-3D Vs model is built on the same uniform grid by multiplying the layer thicknesses of the measured V s profile by S f0. Fourth, 1D GRAs are performed for each mapped Vs profile and their results are statistically combined as a means to incorporate spatial variability into 1D analyses . Finally, t he procedur e can be repeated to incorporate variability across different spatial areas. Note that steps 1 through 3 ( Figure 5) allow for building the pseudo- 3D Vs model s and do not depend on the chosen GRA simulation (i.e., 1D, 2D, or 3D), whereas steps 4 and 5 a re related to the implementation in 1D GRAs . Steps 1 through 3 are explained in greater detail in Hallal and Cox ( 2020) and we only discuss the implementation in 1D GRAs hereafter . Several researchers have shown that by combining 1D TTFs from a suite of varying Vs profiles, an averaging effect of the site response similar to that caused by wave scattering is achieved (de la Torre et al., 2019; Tao and Rathje, 2019; Teague et al., 2018) ; we follow the same approach herein. TTFs are computed for all scaled Vs profile s in the grid following the same procedure s as those implemented to compute the TTF of the measured V s profile . Then, the TTFs are statistically combined and the result s are compared to the ET F. We choose to combine TTFs 12 using both the geometric and arithmetic means. While g eometric means (i.e., lognormal medians) are common in earthquake applications, they are most appropriate when trying to average samples that are dependent on each other . For example , when combining ETFs from different ground motions at a downhole array , the sampled ETFs are dependent , as they are influenced by the same subsurface conditions . Howe ver, by dividing the spatial area into a grid and independently performing a GRA for each cell , we are assuming that the overall site response can be represented by combining results from the independently evaluated cells . In this case , where we are combining independent samples , the arithmetic mean might be more suitable (Cooper, 1996) . For simplicity, we hereafter refer to the geo metric mean as median, which is true for the commonly assumed lognormal distribution, and to the arithmetic mean as only mean. Figure 5. Schematic representation of the H/V geostatistical approach used to build pseudo- 3D Vs models and its implementation to incorporate spatial variability into 1D ground response analyses (GRAs). 13 The mean/median TTF calculations can be repeated by incorporat ing Vs profiles from different spatial areas surrounding the site (Figure 5). This is particularly valuable if we are interested in determining how large of an area influences site response. Additionally, considering different spatial areas may be important when modeling site response under different shaking intensities . For example, Assimaki et al. ( 2003) concluded that the contribution of spatial variability in soil properties strongly depends on the motion intensity. A s the intensity increases, the induced soil nonlinearity tends to eliminate the effects of the variability in soil properties. On the contrary, Thompson et al. ( 2009) noted that larger earthquakes sample a larger volume of the subsurface than sm aller earthquakes (due to generation of longer wavelength/lower frequency waves) , resulting in a larger degree of wave scattering and thus more profound effect of spatial variability on the site response. Implementation at Case Study Sites The pseudo- 3D Vs models developed in Hallal and Cox ( 2020) for the TIDA and DPDA sites using the H/V geostatistical approach are investigated herein using the 1D GRA framework discussed earlier. First, we present results from incorporating the full pseudo- 3D Vs model at each site. Second, we vary the incorporated area to assess the possibility of determining the area influencing site response at each site. Results from Incorporating the Entire Investigated Areas Treasure Island Downhole Array The pseudo- 3D Vs model for the TIDA , along with t he results obtained from statistically combining the TTFs of the scaled 1D V s profiles from the entire investigated area , are illustrated in Figure 6. A total of 646 50-m cells were considered in the analysis , resulting in a pseudo- 3D model wi th approximate surface dimensions of 1.6 km by 1.0 km. It is clear that most of the scaled Vs profiles are comparable to the measured one, implying minimal variability across the island . Additionally, the variability is skewed toward a shallower soil -rock interface (i.e., higher f0,H/V values) near the causeway to YBI in the south (Figure 6a) , which implies stiffer subsurface conditions, consistent with the observations made earlier in regards to the outcropping bedrock on YBI (refer to Figure 3a). Figure 6c compares t he 1D TTFs obtained from the H/V geostatistical approach at the TIDA site to the median ETF ± σ lnETF. In scrutiniz ing Figure 6c, several important observations can be made, which demonstrate that incorporating spatial variability using the H/V geostatistical approach is in fact yielding accurate site response estimates . It is clear, by visual examination , that the mean and median TTF s associated with the H/V geostatistica l approach better match the median ETF than t he single TTF of the 1D homogeneous invasive Vs profile (Figure 2a). Relative to the 1D site response of the measured V s profile , a clear reduction in the amplitude of the peaks and broadening in their widths are visible , suggesting an effect similar to wave scattering and redistribution of the seismic energy when incorporating 3D spatial variability . It should also be noted that the bandwidth broadening effect is more evident in the mean TTF , where a sm all secondary peak near the fundamental mode is evident at the exact same frequency as observed in the median ETF. 14 The quantitative goodness -of-fit parameters (r and m TF) listed between brackets in the legend of Figure 6c also reflect the generally excellent agreement between the statistically combined 1D TTFs and the ETF , with a 30% to 50% improvement in both parameters compared to the 1D homogeneous site response (Figure 2a). The r values for both the median and mean TTFs are higher than 0.9, with a value of 1.0 indicating a perfect relation between the ETF and TTF. It should also be noted that these r values are higher than those achieved by Tao and Rathje ( 2019) Figure 6. The pseudo- 3D Vs model developed for the TIDA using the H/V geostatistical approach and the result of its implementation in 1D GRAs. The uniformly estimated f 0,H/V map is shown in (a) and the spatially variable Vs profiles are shown in (b) using a 1D representation. The results from implementing the pseudo- 3D Vs model in 1D GRAs are shown in (c), in which the median ETF and its +/ - uncertainty bounds, the individual TTFs of the scaled Vs profiles, and their median and mean TTFs are compared. The r and m TF values between the experimental and theoretical transfer functions are indicated in brackets in the legend of (c). The scaled Vs profiles and the individual TTFs ar e colored by the f 0,H/V values of their corresponding cells shown in (a). TIDA (a)-500 0 500 1000 X(m)-200 0 200 400 600Y (m) 0.6 0.8 1.0 1.2f0,H/V(Hz)0 1000 2000 3000Shear Wave Velocity, Vs (m/s) 0 20 40 60 80 100 120Depth (m) (b) 0.5 1 3 10 Frequency (Hz)10-1100101102TF Amplitude (c)Med. ETFlnETF Individual TTFs H/V Geostat. Med. TTF H/V Geostat. [r=0.92 ; mTF=1.13] Mean TTF H/V Geostat. [r=0.91 ; mTF=0.80]Individual Scaled Measured 15 and Afshari and Stewart (2019) , who both obtained a maximum r value of 0.82 when attempting to account for spatial variability using modifi cations to small- strain damping. Additionally, the misfit value s are close to 1.0, imply ing that, on average, the combined TTFs are within approximately 1 standard deviation from the median ETF. These results are very promising, considering the simplistic approach adopted to incorporate spatial variability into 1D GRAs . Delaney Park Downhole Array The pseudo- 3D Vs model shown in 1D along with t he results obtained from statistically combining the TTFs of the scaled 1D V s profiles from the entire investigated area at the DPDA site are presented in Figure 7. A total of 1,024 50-m cells were considered in the analysis , resulting in a pseudo- 3D model with a footprint of 1.6 km by 1.6 km . The scaled Vs profiles show that the variability is relatively significant, with many Vs profiles having shallower or deeper glacial till. However, this variability in the scaled Vs profiles is somewhat symmetric about the measured Vs profile , with till being located, at most, approximately 10 m shallower or deeper across the extent of the considered area. Figure 7c compares the 1D TTFs obtained from the H/V geostatistical approach at the DPDA site to the median ETF ± σ lnETF. Similar to observations at the TIDA site, the geostatistical approach provides a superior fit to the ETF than that computed from the 1D homogeneous site response based on a single V s profile (Figure 2b). This is supported by the goodness -of-fit parameters indicated in Figure 7b. However, while the r value for both the mean and median TTFs increased by about 40% to values above the threshold of good predictions set by Thompson et al. (2012) (i.e., r = 0.6) , only the mTF value of the median improved by about 30% compared to the homogeneous response (Figure 2b), while that of the mean is almost unchanged. Although the amplitude s of the peaks of the mean TTF are reduced, th e amplitudes of the shoulders in the vicinity of the fundamental and higher modes increase d due to broadening of the peaks , resulting in the relatively unchanged misfit . We also note that s imilar to results at TIDA, the r values obtained from the geostatistical approach are higher than those achieved by Tao and Rathje (2019), who obtained a val ue of 0.67 using the D min multiplier approach . Spatial Area Influencing Site Response The initial implementation of the H/V geostatistical approach presented above did not make any explicit attempts to account for a specific spatial area, or to incorporate spatial variability over different scales. Rather, the incorporated area at each site w as selected somewhat arbitrarily to be symmetric around the downhole array sites, when possible, and to extend as far as the available f0,H/V measurements. In this section, we vary the incorporated spatial area and assess the possibility of determining an optimal area that gives the best match to the median ETF. In the companion paper (Hallal and Cox, 2020) , we showed that the descriptive statistics of the spatial variabilities (i.e., standard deviation and skew ness) change with incorporated area, and hence, there is a need to determine the area of interest likely influencing site response. We first start with the point measurement (i.e., the homogeneous 1D response of the measured V s profile ). Then , we increment ally increase the area surrounding the Vs profile by one cell ( i.e., by 50 m ) in each of the four directions (north, south, east, and west). This can be best implement ed at sites which have a symmetric estimated f 0,H/V map, such as that o f the DPDA. For sites with certain limitations that cause the map to be asymmetric, such as the TIDA, we keep on increasing the incorporated area in the directions that can be extended, until covering the entire estimated f 0,H/V map. 16 Deciding on the optimal area is not a trivial task. While the quantitative goodness -of-fit parameter s (r and m TF) aid in this judgement , it can be difficult to evaluate the overall agreement between ETF and TTF based on these parameter s only. Teague et al. ( 2018) mentioned that even with similar goodness -of-fit parameter s, visual examination can arguably show that one TTF is a superior fit. Thus , the incorporated area that gives the best match to the observations will be selected based on a collective assessment of the qual itative shape of the TTF and the quantitative Figure 7. The pseudo -3D Vs model developed for the DPDA using the H/V geostatistical approach and the result of its implementation in 1D GRAs. The uniformly estimated f 0,H/V map is shown in (a) and the spatially variable Vs profiles are shown in (b) using a 1D representation. The results from implementing the pseudo- 3D Vs model in 1D GRAs are shown in (c), in which the median ETF and its +/ - uncertainty bounds, the individual TTFs of the scaled Vs profiles, and their median and mean TTFs are compared. The r and m TF values between the ETF and TTFs are indicated in brackets in the legend of (c). The scaled Vs profiles and the individual TTFs are colored by the f 0,H/V values of their corresponding cells shown in (a). DPDAN (a)-800 -400 0 400 800 X (m)-800 -400 0 400 800Y (m) 1.0 1.2 1.4 1.6 1.8f0,H/V(Hz)0 500 1000 1500Shear Wave Velocity, Vs (m/s) 0 10 20 30 40 50 60Depth (m) (b) 0.5 1 3 10 Frequency (Hz)10-1100101102TF Amplitude (c)Med. ETFlnETF Individual TTFs H/V Geostat. Med. TTF H/V Geostat. [r=0.76 ; mTF=1.62] Mean TTF H/V Geostat. [r=0.74 ; mTF=2.26]Individual Scaled Measured 17 r and m TF values. Additionally, we aim to choose the smallest area that gives a reasonable fit ; meaning, if increasing the incorporated area only results in a modest improvement to the shape of the TTF and the r and m TF values, we select the smaller area as being optimal. Below, w e only present results from the mean TTF since it is believed to be a better summary statistic for data description when combining independent events, which is inherently assumed herein. Treasure Island Downhole Array The variation of the r and m TF values as a function of the incorporated area at the TIDA site is presented in Figure 8a and 8b, respectively. Additionally, for clarity, in Figure 8c we only show the mean TTF at a few selected areas; th ose areas are highlighted in Figure 8a and 8b and outlined in the uniformly estimated f0,H/V map in Figure 8d. Furthermore, the TTFs in Figure 8c are only shown over the frequency range for which the goodness -of-fit parameters were computed (refer to the end of the section titled General Framework of Linear -Viscoelastic GRAs ). It is evident from Figure 8 that incorporating spatial variability at e ven moderate scales provides considerable improvement in the quantitative goodness -of-fit parameters . For example, when considering an area of 400 m x 400 m (i.e., the area highlighted in teal ; 0.16 km2) the r and m TF values are significantly improved. For larger incorporated areas , the quantitative parameter s remain relatively unchanged until an area of approximately 1 km2 is incorporated, after which the r values increase and the m TF values decrease, indicating better goodness -of-fit. The best results are achieved when incorporating the full area of approximately 1,600 m x 1,000 m ( i.e., the area highlighted in magenta ; 1.6 km2). However, it should be noted that the better goodness -of-fit parameters are predominantly driven by the increased widt h of the fundamental mode peak, as indicated in Figure 8c, which shows that the mean TT F of the full area better predict s the width of the fundamental mode and the secondary peak near it (refer to Figure 2a). Thus, while it is evident that shallowing bedrock near YBI is influencing the site response recorded approximately 1 km away at the TIDA, most of the improvement in the mean TTF can be achieved by incorporating a smaller area (e.g., the area highlighted in teal , which is 400 m x 400 m) . In Figure 8a, w e also show r values obtained by researchers who have attempted to incorporate spatial variability into 1D GRAs at the TIDA . The r values shown include results from the D min multiplier approach (Tao and Rathje, 2019) and three different small- strain damping models (Afshari and Stewart, 2019) : Vs-based Dmin model, lab -based D min model, and κ -informe d Dmin model . We note that while all studies use d the same V s profile to evaluate TTFs , some differences in the r values are expected simply due to the different approaches used by each research group to calculate ETFs and the frequency range over which the r calculations were made. Nonetheless, t he results show that the H/V geostatistical approach provides superior r values in comparison to the damping- based methods when even moderate spatial variability is incorporated using the pseudo- 3D Vs model . Thus, by incorporating site -specific spatial variability based on f 0,H/V into simple 1D GRAs , we are capable of better predicting the ground motions observed by the TIDA . Delaney Park Downhole Array The results of the spatial area assessment at the DPDA site are presented in Figure 9. From Figure 9a and 9b, it is evident that the goodness -of-fit parameter s improve up to an area of 800 m x 800 m (i.e., the area highlighted in light pink; 0.64 km2). Beyond this area, the r value show s minimal change , whereas the m TF increases, indicating poorer agreement with increasing incorporated area. Visual examination of the TTFs ( Figure 9c) also indicates that increasing the 18 area beyond 800 m x 800 m insignificantly affects the mean TTF. We note that the estimated f 0,H/V map has high er uncertainty beyond around 500 m from the downhole array site due to the limited measurements along some azimuths , and thus the results could plausibly change with increased future measurements . In Figure 9a, we also show the r value obtained fr om the D min multiplier approach used by Tao and Rathje ( 2019) . Note that Afshari and Stewart ( 2019) only studied California downhole arrays and thus have no results for DPDA. Similar to TIDA, a better match to the ETF is achieved by using the H/V geostatistical approach to account for spatial variabil ity in 1D GRAs at the DPDA site, provided an area greater than approximately 400 m x 400 m is used (i.e., the area highlighted in teal, which is 0.16 km2). We note that Ta o and Rathje ( 2019) used a slightly modified version of the Vs profile by Nath et al. ( 1997) , which is different than the one used in this study Figure 8. Effect of varying the incorporated spatial area in the pseudo -3D Vs model at the TIDA site on the: (a) r value, (b) m TF value, and (c) mean TTFs for selected areas plotted only over the frequency range used in the goodness -of-fit parameters’ calculation. (d) The estimated f 0,H/V map showing the selected areas whose results are highlighted in (a) and (b) and displayed in (c ) using the same colors. In (a), the r values obtained from the D min multiplier approach by Tao and Rathje (2019) and three different damping models by Afshari and Stewart (2019) are shown for comparison. 0.6 0.7 0.8 0.9 1.0r -informedDmin3 x Dminmultiplier & -informedDmin VS-based DminLab-based Dmin Vs-based Dmin(a) 0 0.5 1 1.5 2 2.5 Area (km2)0.6 0.8 1.0 1.2 1.4 1.6 1.8mTF(b) 1 3 Frequency (Hz)10-1100101102TF Amplitude(c) Med. ETFlnETF TTF PS [r=0.63 ; mTF=1.63] Mean TTF Geostatistics at selected areas -500 0 500 1000 X (m)-200 0 200 400 600Y (m) TIDA(d) 0.6 0.8 1.0 1.2f0,H/V(Hz)Med. ETFlnETF TTF PS [r=0.63 ; mTF=1.63] Mean TTFs H/V Geostat. at selected areas 19 (comparison of Vs profiles not shown in this paper; refer to Hallal and Cox 2020) . However, by examining the TTF reported by Tao and Rathje ( 2019) , it is evident that it overpredicts the amplitude of the fundamental mode and underpredicts higher modes. Significant underprediction of higher -mode amplitudes tends to occur when large increases in small -strain damping are used as a means to better match the amplitude of the fundamental mode. While the H/V geostatistical approach also overpredicts the amplitude of the fundamental mode, the amplitudes of higher modes are comparable to those of the median ETF (Figure 9c). Figure 9. Effect of varying the incorporated spatial area in the pseudo -3D Vs model at the DPDA site on the : (a) r value, (b) m TF value, and (c) mean TTFs for selected areas plotted only over the frequency range used in the goodness -of-fit parameters’ calculation. (d) The estimated f 0,H/V map showing the selected areas whose results are highlighted in (a) and (b) and displayed in (c) using the same colors. In (a), the r value obtained from the D min multiplier approach by Tao and Rathje (2019) is shown for comparison. 0.5 0.6 0.7 0.8r6 x DminMultiplier(a) 0 0.5 1 1.5 2 2.5 Area (km2)2.0 2.1 2.2 2.3 2.4 2.5mTF(b) 3 10 Frequency (Hz)10-1100101102TF Amplitude(c) Med. ETFlnETF TTF DH [r=0.54 ; mTF=2.34] Mean TTF Geostatistics at selected areas -800 -400 0 400 800 X (m)-800 -400 0 400 800Y (m) DPDA(d) 1.0 1.2 1.4 1.6 1.8f0,H/V(Hz)Med. ETFlnETF TTF DH [r=0.54 ; mTF=2.34] Mean TTFs H/V Geostat. at selected areas 20 Comparison of R esults at Case Study Sites The case study sites clearly represent areas with different degrees of spatial variability. For example, the 1D homogenous GRA based on a single invasive Vs profile yields better results at TIDA than at DPDA ( Figure 2). This was attributed to the fact that TIDA is an embodiment of local 1D site conditions, whereas DPDA is not. However, the results of this paper demonstrate that because of the significantly large area influencing site response, non- 1D aspects are inevitable, even at the TIDA site. Even so , depending on how variable the site conditions are, basic 1D GRAs can provide acceptable predictions if spatial variability is incorporated in a logical way . While it seems that the H/V geostatistical approach yields better predictions at sites closer to 1D conditions (i.e., the TIDA site), the goodness -of-fit parameter s are relatively equally improved by around 30% – 50% for both case study sites , which distinctly address simple and complex subsurface conditions at TIDA and DPDA, respectively . The assessment of incorporated spatial area showed that TIDA requires an area greater than 1.6 km2 (i.e., the area highlighted in magenta) to achieve the absolute best fit to the median ETF (Figure 8). However, very good fits could be achieved with an incorporated area less than approximately 0.16 km2 (i.e., the area highlighted in teal). On the other hand, DPDA requires an incorporated area of approximately 0.64 km2 (i.e., the area highlighted in light pink) to achieve an optimal ETF fit when considering both r and m TF values (Figure 9). Yet, even when the incorporated area is optimized at DPDA, the goodness -of-fit parameter s at DPDA are still poorer than those obtained from a 4- times smaller, non -optimized area of 0.16 km2 at TIDA. As such, DPDA will likely require more advanced 2D/3D G RAs in order to better match the recorded small - strain site response. We believe that the spatial area influencing site response at any given site depends on several factors , including: the depth of the soil -rock interface that governs f 0, the degree of th e spatial variability in material properties , and the frequency/ wavelength of the earthquake waves . We find it intuitive that sites with deeper soil- rock interface (i.e., lower f 0 values) are likely influenced by a larger area. Hence, TIDA would be expected to be influenced by a larger area than DPDA , which it appears to be . However, sites with less spatial variability in material properties will tend to behave more like 1D sites . As such, 1D GRAs would be expected to perform better at TIDA than at DPDA, wh ich they do. Finally, we believe that shorter wavelengths, which correspond to higher frequencies and softer subsurface conditions, sample a smaller portion of the subsurface and consequently should be influenced by spatial variability over a smaller area , while longer wavelengths, which correspond to lower frequencies and stiffer subsurface conditions, sample a larger portion of the subsurface and consequently should be influenced by spatial variabilit y over a larger area. This means that the effects of spatial variability are most likely frequency dependent. With only two case study sites, it is difficult to assess exactly how these factors combine to affect th e area influencing site response . However, at both sites the statistical results obtained from 1D GRAs are significantly improved when an area of at least 400 m x 400 m (i.e., the area highlighted in teal in Figure s 8 and 9; 0.16 km2) is incorporated, and even larger incorporated areas produce bet ter results at both sites. Thus, this size of an area might be considered as a minimum over which to account for spatial variability in GRAs. Further studies will be required to refine this approximation, to investigate how ground motion intensity affects this approximation, and to determine if a relatively larger or smaller area is needed when performing 2D/3D GRAs instead of 1D GRAs. 21 Interestingly, we observe that the change in the goodness -of-fit parameter s with incorporated area follows similar trends as that of the lognormal standard deviation in f 0,H/V (σlnf0,H/V) (results not shown in this paper, refer to Figure 11a in Hallal and Cox 2020) ). For example, the r value at TIDA rapidly increases at first, f ollowed by a nearly flat response, beyond which it increases again. The exact same trend is observed for σ lnf0,H/V. Similarly, we observe the same trend in r and σ lnf0,H/V at DPDA. This is expected , as changes in the TTF, and hence the goodness -of-fit parameter s, are driven by spatial variability, which is measured herein using f 0,H/V. When the uniformly estimated f0,H/V shows minimal changes (e.g., going from the teal to the light pink area at TIDA in Figure 8d) the Vs profiles also change insignificantly, and consequently their combined TTF remains unchanged. In engineering practice, where empirical recordings of ground motions are not available, engineers might gain insights on the area likely influencing site re sponse at their site by assessing the vari ation of σ lnf0,H/V as a function of incorporated area. Owing to the inherent geologic spatial autocorrelation, larger areas tend to have dissimilar properties. Thus, it is likely that σ lnf0,H/V will increase at first with incorporated area. If σlnf0,H/V then shows minimal chang e beyond a certain area, then th e additional incorporated area might not be of interest for GRAs, as the site response is also likely to remain unchanged. On the contrary, if σlnf0,H/V keeps on increasing, or increases /decreas es beyond a flat region, then determining whether this additional incorporated area is of interest to site response will require an under standing of the factors discussed earlier (i.e., the depth of the soil -rock interface governing f0, the degree of spati al variability in material properties, and the frequency/wavelength of the earthquake waves ). However, these hypothes es require further investigation by studying more downhole arrays , such that we might provide better guidance on how large an area should be incorporated into GRAs when only experimental site characterization measurements are available. Conclusions The H/V geostatistical approach proposed for building pseudo- 3D Vs models using simple, site-specific measurements has been implemented into a 1D GRA framework at two downhole array case study sites with recorded surface and rock ground motions, enabling comparisons between predicted and recorded site response . While the goal of the H/V geostatistical methodology and its application to 1D GRAs is not necessarily to “ perfectly ” predict the site response for a site, it may be reasonably successful at doing so for sites without significant spatial variability . However, the primary advantage of using the proposed approach is to incorporate V s profiles t hat realistically represent the spatial variability at the site and avoid the inclusion of profiles that are unrealistic due to assumed stochastic parameters that are not site- specific. Additionally, the application of the H/V geostatistical approach to 1D GRAs allows engineers to still effectively use simple 1D analyses at many sites . Another key advantage of the implementation of this approach to 1D GRAs is the ability to predict site response as a function of the spatial variability across different footprints. At both downhole array sites, GRAs that incorporated the site -specific variability from the H/V geostatistical approach provided a superior fit to the recorded ground motions compared to other common approaches that also attempt to account for spa tial variability in 1D GRAs . While th e 22 proposed approach requires the consideration of the variable subsurface conditions from a certain spatial footprint around the site, with no current strict guidelines on the extent of such area, our results show that 1D GRAs are significantly improved when an area of at least 400 m x 400 m (i.e., 0.16 km2) is incorporated, and even larger incorporated areas produce d better results at the two studied downhole array sites . Though a simplified pseudo- 3D imaging technique is used to model th e subsurface conditions on the basis of a single V s profile and a uniformly estimated f0,H/V map from geostatistical kriging, the proposed approach is capable of reproducing phenomena observed in the ETF, such as reduced amplificati on due to attenuation of seismic waves through scattering , increased amplification at frequencies different from the expected resonance frequencies, and redistribution of the seismic energy over a wider frequency band. Without considering site-specific inf ormation about the lateral variability in the subsurface conditions, it is generally not possible to know whether the predicted site response is adversely influenced by the assumed spatial variability model parameterizatio n. The H/V geostatistical approach proposed in the companion paper provides a systematic means to ensure the incorporation of site -specific spatial variability that governs site response in a meaningful and practical way. Furthermore, it has been shown that implementing the approach within the 1D GRA framework can provide insights on the potential spatial area influencing site response. More studies are needed at other downhole array sites to determine if the H/V geostatistical approach will work reasonably w ell for more complicated site conditions that contain, for example, multiple strong impedance contrasts. Additionally , the pseudo- 3D Vs models should be investigated in 2D/3D GRAs , at least for sites with poorer model predictions, such as DPDA, to achieve the ultimate goal of improving the modeling of site response in engineering practice. T he proposed implementation of the geostatistical approach within 1D GRAs provides a promising framework to account for site -specific spatial variability in V s using simp le 1D analyses without the need to perform advanced site response analyses at some site s, although further validation is of interest . Acknowledgements The H/V data from DPDA and TIDA were collected with funding from Pacific Gas and Electric (PG&E). However , any opinions, findings, conclusions, or recommendations expressed in this paper are those of the authors and do not necessarily reflect the views of PG&E . We would like to thank Mr. Joseph Vantassel and Mr. Michael Yust for their help in planning and col lecting the H/V data at DPDA and TIDA, Dr. Albert Kottke, Dr. David Teague, and Dr. Krishna Kumar for their help in collecting the H /V data at TIDA , and Ms. Jodie Crocker for her help in accessing and processing the ground motions .