Many recent studies have shown that we are generally unable to accurately replicate recorded ground motions at most borehole array sites using available subsurface geotechnical information and one -dimensional (1D) ground response analyses ( GRAs ). When 1D GRAs fail to accurately predict recorded site response, the site is often considered too complex to be effectively modeled as 1D. While 3D numerical GRAs are possible and believed to be more accurate, there is rarely a 3D subsurface model available for these analyses. The lack o f affordable and reliable site characterization methods to quantify spatial variability in subsurface conditions, particularly regarding shear wave velocity ( Vs) measurements needed for GRAs , has pushed researchers to adopt stochastic approaches, such as Vs randomization and spatially correlated random fields. However , these stochastically generated models require the assumption of generic, or assumed , input parameters , introducing significant uncertainties into the site response predictions . This paper des cribes a new geostatistical approach that can be used for building pseudo- 3D Vs models as a means to rationally account for spatial variability in GRAs , increase model accuracy , and reduce uncertainty . Importantly, it requires only a single measured Vs profile and a number of simple, cost -effective, horizontal -to-vertical spectral ratio (H/V) noise measurements. Using Gaussian geostatistical regression , irregularly sampled estimates of fundamental site frequency from H/V measurements (f 0,H/V) are used to generate a uniform grid of f0,H/V across the site with accompanying Vs profiles that have been scaled to match each f 0,H/V value, thereby producing a pseudo- 3D Vs model. This approach is demonstrated at the Treasure I sland and D elaney Park Downhole Array sites (TIDA and DPDA, respectively) . While the pseudo- 3D Vs model s can be used to incorporate spatial variability into 1D, 2D , or 3D GRAs , their implementation in 1D GRAs at TI DA and DP DA is discussed in a companion paper. 2 Introduction Many recent studies ( e.g., Afshari and Stewart, 2019; Kaklamanos et al., 2013, 2015; Kaklamanos and Bradley, 2018; Tao and Rathje, 2019; Teague et al., 2018; Thompson et al., 2009, 2012; Zalachoris and Rathje, 2015) have shown that we are generally unable to accurately replicate recorded ground motions at most borehole array s ites using available subsurface geotechnical information and one -dimensional (1D) ground response analyses ( GRAs). When 1D GRAs fail to accurately predict recorded site response, the site is often considered too complex to be effectively modeled as 1D. Whi le 3D numerical GRAs are possible and believed to be more accurate, there is rarely a 3D subsurface model available for these analyses. T he lack of affordable and reliable site characterization methods to quantify spatial variability in subsurface conditions, particularly regarding shear wave velocity (Vs) measurements needed for GRAs, is an impediment to future progress in accurate site response modeling. Indeed, even at our most valued downhole array sites in the U.S., Japan, and elsewhere globally, we ar e generally limited to 1D representations of Vs. Common Vs characterization technique s can be subsumed under two categories: invasive /borehole methods and noninvasive /surface wave methods. Invasive methods (e.g., downhole, crosshole, PS logging) require placing the seismic source and/or receiver(s) in the subsurface material being characterized . These methods produce 1D Vs profiles at the measurement location that are generally believe d to be more accurate than those from noninvasive methods, even though they are often more time consuming and expensive. N oninvasive methods (e.g., active - and passive -surface wave measurements, full waveform inversion) use arrays of receivers placed across the ground surface , and thus sample a broad er spatial area . However, most surface wave approaches rely on the basic assumption that the medium beneath the array can be approximated by 1D geometry (Foti et al., 2011) . In the presence of variable subsurface conditions, surface wave methods result in a complex averaging of material properties beneath the array. This fact, combined with the nonunique nature of the surface wave inverse problem , results in 1D Vs profiles that are generally believed to be less accurate than those from invasive methods, even though they are often a cost and time efficient alternative for deep profiling. Regardless of whether invasive or noninvasive methods are used for dynamic site characterization , typically a single deterministic 1D Vs profile is obtained for each measurement location . In the absence of additional information, these 1D Vs profiles are used to perform 1D GRAs , which assume that all subsurface layers are isotropic , laterally homogeneous , and extend infinitely. Researchers have shown that th ese assumptions greatly limit the effectiveness of 1D GRAs due to the laterally variable conditions characterizing the shallow geology of most sites (Afshari and Stewart, 2019; Assimaki et al., 2003; Kaklamanos and Bradley, 2018; Tao and Rathje, 2019; Thompson et al., 2009, 2012) . Depending on the extent of spatial variability, 1D GRAs can be inadequate and incapable of capturing complex effects , such as wave scattering, locally generated surface waves, and basin effects, among others (Kaklamanos et al., 2013) . In an attempt to replicate the inherent spatial variability in Vs within the 1D G RA framework , researchers proposed using probabilistic statistical models (e.g., Toro, 1995) to randomly vary the properties of a reference 1D Vs profile (Kaklamanos and Bradley, 2018; Rathje et al., 2010; Tao and Rathje, 2019) . While this approach still relies on 1D GRAs with assumed laterally homogeneous subsu rface layers, the combined result of all randomly generated Vs profiles can 3 reproduce wave scattering effects (Nour et al., 2003) . However, th is method rel ies on stochastic probabilistic modeling approaches and require s site-specific input parameters that are not typically available, bringing additional uncertainty to the predictions (Passeri et al., 2 019; Teague et al., 2018) . This uncertainty can be significant , and blind Vs randomization has been shown to yield unrealistic Vs profiles that lead to poor site response predictions (Griffiths et al., 2016a, 2016b) that are different than ground motions recorded at downhole array sites (Teague et al., 2018) . Additionally, even if randomization methods can reproduce wave scattering effects, other compl ex phenomena caused by highly variable subsurface conditions cannot be captured. Recognizing the imperfect ability of modified 1D GRAs to account for 3D spatial variability, advanced 2D/3D GRAs are slowly gaining popularity. Such analyses require accurate information on the 2D/ 3D geometry and material properties of the subsurface, specifically the Vs . Within this context, s everal methods for constructing 2D/3D Vs models using in- situ measurements have been developed. The simplest approach is to utilize multiple spatially distributed 1D Vs profiles to construct pseudo -2D (e.g., Luo et al., 2008; Mi et al., 2017; Parolai et al., 2001) and 3D (e.g., Pan et al., 2018; Piatti et al., 2013; Strobbia et al., 2011) Vs models by spatial interpolation techniques . For example, this approach has been applied to construct pseudo- 2D Vs profiles from multi- channel analysis of surface waves (MASW) by interpolating a series of 1D Vs models derived from inversion of waveforms collected using different receiver arrays (e.g., Ismail et al., 2014; Mahvelati and Coe, 2017; Park, 2005) . While th is method can be used to produce pseudo- 2D/3D Vs models, the input 1D velocit y profiles repre sent spatial averages over the length of the receiver array , leading to limited resolution. Additionally, this method require s significant numbers of array/ Vs measurement locations with good azimuthal coverage around the site , which renders it extremely time consuming and cost prohibitive. Other methods have been proposed to obtain true 2D/3D Vs model s with higher resolution, with the most common being tomography and full waveform i nversion (FWI). Tomographic methods are based on cross -correlating and summing ambient noise measurements from pairs of sensors to assemble 2D frequency -dependent group or phase velocity maps, which are then inverted to obtain 3D Vs models. The method has been applied on continental and regional scales using very low frequency waves to image the Earth’s crust and upper mantle (Chen et al., 2018; Ritzwoller and Levshin, 1998; Yang et al., 2007) and has been extended to higher frequencies to image the near surface Vs at local scales (Picozzi et al., 2009; Pilz et al., 2012) . Rather than inverting dispersion information, FWI attempts to reconstruct the subsurface model that could produce theoretical waveforms m atching the entire wavefield measured by a large array of sensors. This requires a full solution of the governing wave equation in 3D. The approach has been applied to characterize the near surface Vs and detect anomalies, such as voids, among other applic ations (Kallivokas et al., 2013; Tran et al., 2013; Tran and Sperry, 2018) . Both methods , however, require solving large inverse problems, which are nonlinear, ill -posed, and nonunique (Pan et al., 2019; Pilz et al., 2012; Virieux and Operto, 2009) . Additionally, in processing and inverting data from tomography, many critical assumptions are made , including, great -circle (i.e., approximately straight) propagation path between sensor pairs , lateral smooth variation of Vs, and uniform distribution of ambient noise . These assumptions are generally violated, which has been shown to introduce bias into tomography (Chen et al., 2018; Picozzi et al., 2009; Ritzwoller and Levshin, 1998; Yang et al., 2007) . FWI also faces challenges pertaining to inconsistent source excitation, variable coupling of sensors, lack of low -frequency waves, and 4 high computational complexity in terms of both time and memory (Tran et al., 2013; Tran and Sperry, 2018; Virieux and Operto, 2009) . Finally, the need for wide -aperture surveys with numerous , simultaneously record ing sensors , which is not practical in urban areas, is another limitation for both tomography and FWI . For these reasons, these methods are seldom applied to obtain 2D/3D Vs models for ground motion studies in practice. The lack of affordable and reliable 3D characterization methods , particularly regarding Vs measurements, has pushed researchers to adopt stochastic approaches to construct 2D/3D Vs models for GRA s. Using a single 1D Vs profile , 2D/3D Vs models are developed by continuously perturbing the 1D information to mimic continuous spatial variability. Stochastic probabilistic modeling approaches have been widely used for this purpose in geotechnical earthquake engineering. The most widely used probabilistic modeling approach utilizes spatially correlated random fields to introduce randomness to a baseline 1D Vs model (e.g., de la Torre et al., 2019; Thompson et al., 2009) . While conceptually simple, generating 3D spatially correlated random fields is numerically and computationally intensive (Thompson et al., 2009) , especially considering that several randomly generated realizations of the subsurface conditions may be needed to assess the uncertainty. Also, similar to 1D Vs randomization methods, these stochastically generated 2D/3D models require the assumption of generic input parameters derived for other sites , which may not be able to account for the site -specific spatial variability in a meaningful way. This inevitably introduc es significant uncertainties into the predictions. The collective result of the issues presented above demonstrate s the need for an affordable method to construct Vs models which account for site -specific spatial variability using approaches that: (1) are not fully stochastic , but site -specific, (2) do not require inputs that cannot be known a priori , and (3) can conveniently develop Vs model s over large areas (hundreds of square meters) without significant acquisition and processing efforts . This paper presents a new framework called the ‘H/V geostatistical approach’ that can be used to develop site-specific pseudo- 3D Vs model s in a rigorous and simple manner. The proposed approach distinguishes itself from previous studies in three key ways: (1) it requires only a single, accurately measured V s profile down to engineering bedrock, (2) it relies majorly on estimates of fundamental site frequency ( f0; a key parameter governing site effects) obtained from simple horizontal -to-vertical spectral ratio (H /V) noise measurements (f 0,H/V), and (3) it creates models that can be used to ensure proper incorporation of site-specific spatial variability in 1D, 2D, and 3D GRAs . In a companion paper (Hallal and Cox, 2020) , the benefits of implementing th e pseudo -3D Vs model s generated using the proposed approach in 1D GRAs are discussed. In this paper, we first discuss the theory behind H/V measurements and their use for inf erring subsurface properties . Second, we describe the general framework of the proposed H/V geostatistical approach in terms of generat ing a uniform grid of f 0,H/V across a certain area with accompanying Vs profiles that have been scaled from the measured Vs to match each f 0,H/V value, thereby producing a pseudo- 3D Vs model . Third, we summarize information on t wo case study sites, the Treasure Island Downhole Array (TIDA) and Delaney Park Downhole Array (DPDA), along with the field studies used to measure the spatial variability of f 0 across each site using widely distributed, noninvasive H/V measurements. F ourth, we present the application of the H/V geostatistical approach for building pseudo- 3D Vs models at the two case study site s. Fifth, we discuss the descriptive measures/statistics of the site -specific spatial variability that should be properly accounted for and conclude with final remarks. 5 Review on H /V Theory The novelty of the proposed approach comes from the integration of H/V measurement s and their strong dependence on subsurface conditions controlling site response with geostatistic al modeling tools . A major component of the proposed approach is performing a number of spatially distributed H/V measurements across and around the site in order to capture trends in the variation of f0. The H /V method, which was first proposed by Nogoshi and Igarashi ( 1971) and later popularized by Nakamura (1989) , is continuously attracting more and more interest in seismic studies due to its recognized ability to estimate f 0 for sites with a strong impedance contrast in a cost-effective, time -efficient, and feasible way. As shown in Figure 1, the H/V method utilizes ambient noise recordings (generally, 15- to 30- minute time records are sufficient) from the horizontal and vertical components of a single, three- component seismometer to estimate f 0. The H/V curves are computed as the ratio of the Fourier a mplitude spectra (FAS) of the horizontal and vertical components of the noise reco rd. While each horizontal component can be considered separately, or rotated to a particular orientation based on horizontal polarization analyses, it is most common to combine the two horizontal components using one of several methods (Albarello and Lunedei, 2013; Cox et al., 2020; Fäh et al., 2001; Haghshenas et al., 2008) . Experimental and theoretical studies have shown that if the H /V spectral ratio curve of ambient noise exhibits a dominant, well -defined peak (or multiple clear peaks, if multiple strong impedance contrasts are present at different depths), then the lowest frequency peak can be used to infer the fundamental shear wave resonan ce frequency at the location of the recording seis mometer (Bard and SESAME -Team, 2004; Lachet and Bard, 1994; Lermo and Chavez -Garcia, 1993) . Owing to this simple acquisition and processing, many spatially distributed H/V measure ments can easily be conducted across a site of interest as a means to infer the variation of fundamental site Figure 1. Simplified schematic of the horizontal- to-vertical spectral ratio (H/V) method used to infer the fundamental site frequency (f 0,H/V). Soil Deposits Rock Vertical Horizontal Sensor H/V Spectral Ratio Curve FrequencyH/V Spectral AmplitudeRatio f0,H/V 6 frequency . As fundamental site frequency is being estimated from H/V measurements, we refer to it as f 0,H/V. Several competing theories exist regarding the basis for the relation between f 0,H/V and f0. One theory hypothesizes that the H/V spectra l ratio curve is a proxy to the Rayleigh wave ellipticity curve, which is simply the ratio of the amplitude of the horizontal particle motion to that of the vertical particle motion (Fäh et al., 2003; Lachet and Bard, 1994) . However, the validity of this hypothesis depends on the composition of the ambient noise wavefield. If the wavefield i s composed exclusively of fundamental mode Rayleigh waves, then the H/V spectral ratio curve is equivalent to the ellipticity curve. In reality, the wavefield i ncludes different modes of Rayleigh, Love, and body waves, which are generally not separated in time. Owing to t his complex wavefield composition, the H/V spectral ratio curve will not exactly re present the ellipticity curve. If the energy partition between the different waves is known, the H/V spectral ratio curve can be corrected to obtain ellipticity (Hobiger et al., 2013; Scherbaum et al., 2003) . However, studies have shown that the energy partition varies from site to site and, even for a given site, is frequency- dependent and changes with time, making th e H/V correction a nontrivial task (Endrun, 2011) . Directly extracting the Rayleigh wave ellipticity would avoid making any assumptions on the composition of the wavefield (Hobiger et al., 2013) . However, a further evaluation of this issue is beyond the scope of this paper. Evaluating the contribution of different wave types is cri tical to applications that attempt to invert H/V spectral ratios to retrieve the Vs profile. It is well established that the frequency - dependency of ellipticity is linked to the Vs profile, and by using its relation to the H/V curve, many researchers have used the H/V curve to obtain Vs profiles directly via inversion (e.g., Parolai et al., 2001; Yamanaka et al., 1994) . However, inaccurate identification of the proportion of fundamental mode Rayleigh waves might severely bias the inversion (Bonnefoy- Claudet et al., 2006) . Additionally , the H/V curve is dimensionless, causing the inversion to be subject to a strong tradeoff between the layer s’ thicknesses and Vs (Scherbaum et al., 2003; Yamanaka et al., 1994) . That is, by scaling the wave velocities and layer thicknesses by the same constant, the ellipticity curve remains exactly the same (Hobiger et al., 2013) . It is thus obvious that an H/V curve alone is not sufficient to obt ain the Vs profile, and additional information, such as wave velocities, layer thicknesses, or dispersion estimates, is needed. Regardless of the composition of the wavefield, researchers have shown that the position of the H/V spectral ratio peak is unaff ected (Bonnefoy- Claudet et al., 2006) . Thus, whateve r the predominant wave type is, H/V spectral ratio curves can be used to infer the fundamental shear wave resonance frequency , at least for sites with a strong impedance contrast. Therefore, we hereafter focus on common methods used to obtain f0,H/V and it s relation to subsurface conditions. Although relatively simple calculations are required for the H /V technique , there exist several deterministic and statistical methods to obtain the f 0,H/V (Cox et al., 2020) . Cheng et al. ( 2020) summarized the existing techniques and proposed a new method that directly accounts for azimuthal variability , as several studies have show n the existence of azimuthal variability in f 0,H/V and attribute it to anisotropic ambient vibration sources and/or complex subsurface conditions which do not fall under the idealized 1D assumptions . If H/V measurements are performed at multiple locations, the extent to which the f 0,H/V varies across the area of interest can be examined and mapped in 2D. The spatial variability of f 0,H/V reflects the degree of variability in the subsurface V s structure that governs site response . This 7 includes variations in the V s of soil layers, their thicknesses, and/ or in the depth to a major impedance contrast. Potential c auses for changes in f 0,H/V can be inferred by examining the equation of the fundamental site frequency of a single representative soil layer with an average shear wave velocity (V s,avg) and thickness (H) overlying a rigid half -space. For such c onditions, f0 can be approximated from the quarter -wavelength relationship , as shown in Equation 1. It is intuitive from this equation that f 0 is directly proportional to V s but inversely proportional to the depth to a major impedance contrast, which is ty pically that of the soil- rock interface. Thus , an increase in f 0 is associated with shallower rock and/or an overall stiffer soil column, and vice versa. This relation ship forms the basis for the H/V geostatistical approach proposed herein to build pseudo- 3D Vs models. 𝑓𝑓0≅𝑉𝑉𝑆𝑆,𝑎𝑎𝑎𝑎𝑎𝑎 4𝐻𝐻 (1) H/V Geostatistical Approach for Building Pseudo -3D Vs Model s The H/V geostatistical approach proposed herein uses site-specific, simple measurements to develop a reasonable and representative pseudo- 3D subsurface model that can be used to incorporate spatial variability into 1D, 2D , and 3D GRAs . The site characterizati on information required to implement this approach is: (1) at least one V s profile down to engineering bedrock (i.e., the impedance contrast governing f0 and site response), and (2) a number of H /V measurements spatially distributed across the site. The appropriate number of H /V measurements will depend on the size and complexity of the site , but should be sufficient to capture trends in the variation of f 0 across the site. Additionally, it may be necessary to conduct H /V measurement s outside the proper boundaries of your site in order to capture spatial variability over the area that is impacting site response. Insights into the area required to account for spatial variability in GRAs are discussed in Hallal and Cox ( 2020) . A schematic summary of the H/V geostatistical approach is shown in Figure 2. First, the irregularly spaced f0,H/V measurements are interpolated (for example, using kriging) to obtain a uniformly estimated map of f 0,H/V values. Second , using the f0,H/V value closest to the measured Vs profile (f0,H/V,Vs ) and the estimated f0,H/V map, a uniformly estimated scaling factor (S f0) map is obtained. Finally , 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 of each cell . The application could be easily extended to include multiple measured Vs profiles if desired . The pseudo- 3D Vs model can be visualized by stacking the scaled Vs profiles from all cells in 1D, by plotting 2D Vs cross sections along certain azimuths, or by viewing different cross sections in 3D (refer to Figure 2). The methodology is explained in greater detail in the subsequent subsections and then demonstrated for two downhole array case study sites. 8 Uniform Spatial Estimation of f 0,H/V The first step in the proposed approach aims at uniformly mapping the variation of f 0,H/V over an area of interest using irregularly sampled H/V measurements. Exhaustive measurements on a uniform grid are not practical , or even possible due to access and/or budget limitations , and thus irregular spatial sampling is always performed. The f 0,H/V values at unsampled locations are then obtained via spatial estimation methods. An in -depth assessment of different estimation techniques is beyond the scope of this paper and we only provide a systematic method rather than the “best” Figure 2. Schematic representation of the H/V geostatistical approach proposed herein to build pseudo- 3D Vs models that can be used to incorporate spatial variability into ground response analyses. 9 one for the desired application. Ultimately, individual users can judge which spatial estimation technique is bes t for their application. In this study, we use the kriging with a trend approach, also known as k riging after detrending (Hengl et al., 2004; Isaaks and Srivastava, 1989) . Compared to other methods, kriging has the advantage of incorporating valuable information regarding the site - specific spatial correlation and anisotropy, in addition to minimizing the error of the estimated values. Moreover , combining kriging with a trend has been proven to yield more accurate results, especially for geotechnical -related applications, where the data sets nearly always show some pattern of spatia l continuity related to geology (Hengl et al., 2004; Isaaks and Srivastava, 1989) . The estimated f 0,H/V map is obtained in a two- step approach, in which the trend and residuals are fitted separately and then summed. We model the trend in the f 0,H/V values using a convolution- based method, which is a moving window average of the measurements. The trend is then subtracted from the data and ordinary kriging is performed on the residuals. This step requires modeling the spatial correlation of the measur ements, or its reverse, which is commonly done using the semiv ariogram . The semivariogram is a measure of the degree of variability between observations as a function of their separation distance. Kriging on the residuals is then performed on a grid with a specified cell size. We use 50 -m cells herein , however smaller cell sizes might be needed for sites which show significant spatial variability within small distances. After performing kriging to map the residuals over a uniform area, the trend is added (Isaaks and Srivastava, 1989) and the final output is a uniformly estimated f 0,H/V map. Building a Pseudo -3D Vs Model The uniformly estimated f0,H/V map is only a means to develop a pseudo- 3D subsurface model that captures the key aspects of the site -specific spatial variability likely to influence site response . To achieve this, a single measured V s profile is repeatedly modified /scaled for each cell such that the modified profile is consistent with the f 0,H/V value of that cell . The task of perturbing the V s profile to match f 0,H/V is not trivial, as change s in f0,H/V result from variations in the layers’ Vs, their thicknesses, and /or the depth to the major impedance contrast (refer to Equation 1). If we make the simplifying assumption that the layer thicknesses and depth to the major impedance contrast are primarily c ontrolling the change in f0,H/V, then the individual layer thicknesses can be altered to achieve the desired f 0,H/V. While th is assumption inherently includes some uncertain ty, it provides a feasibl e way of economically developing representative Vs profile s over a large footprint . Additionally, p revious research ha s shown that if V s profiles match the experimental site signature, which is implemented here by targeting f0,H/V, then the uncertaint y in GRAs is greatly reduced (Griffiths et al., 2016a; Teague et al., 2018) . Since we are scaling the measured V s profile to a target f0,H/V at each grid point , we need to represent the measured V s profile with an f 0,H/V value, which we refer to as f 0,H/V,V s. The f 0,H/V,V s is simply obtained from the H /V measurement closest to the location of the measured V s profile (Figure 2). By having an initial f 0,H/V ,Vs and a target f0,H/V,i for the ith cell in the grid , we can compute a simple scaling factor as the ratio of these two values for that cell (S f0,i), as shown in Equation 2. 𝑆𝑆𝑓𝑓0,𝑖𝑖=𝑓𝑓0,𝐻𝐻/𝑉𝑉,𝑉𝑉𝑠𝑠 𝑓𝑓0,𝐻𝐻/𝑉𝑉,𝑖𝑖≅�𝑉𝑉𝑆𝑆,𝑎𝑎𝑎𝑎𝑎𝑎 4𝐻𝐻� 𝑉𝑉𝑠𝑠 �𝑉𝑉𝑆𝑆,𝑎𝑎𝑎𝑎𝑎𝑎 4𝐻𝐻� 𝑖𝑖=𝑉𝑉𝑆𝑆,𝑎𝑎𝑎𝑎𝑎𝑎 4𝐻𝐻𝑉𝑉𝑠𝑠 𝑉𝑉𝑆𝑆,𝑎𝑎𝑎𝑎𝑎𝑎 4𝐻𝐻𝑖𝑖=𝐻𝐻𝑖𝑖 𝐻𝐻𝑉𝑉𝑠𝑠 (2) 10 Vs profiles are finally developed for each cell by s imply multiplying the layer thicknesses of the measured V s profile by Sf0,i of each cell (Figure 2). This ensures that the thickness of the soil above the major impedance contrast of the measured V s profile (H Vs) will be scaled by S f0,i and consequently, as shown in Equation 2, f 0,H/V,Vs will be scaled by 1/ Sf0,i, achieving the target f 0,H/V,i . This is true because we are assuming that the measured Vs remains constant for each layer and that the thicknesses of all layers are estimated by the same scaling factor, implying that the V s,avg of the soil column above the major impedance contrast also remains constant for all cells . We acknowledge that although variability in the layers’ V s or thickness scaling can be introduced, this causes the problem of modifying the V s profile to be ill -posed with infinite possibilities. To avoid the need to account for this non- uniqueness and additional uncertainty, we decided not to include it as part of this study. One approach to avoid assuming either con stant Vs or constant thickness scaling for each layer is to invert the measured H/V curves while keeping the layer s’ thicknesses or Vs fixed, respectively . Recall that inverting the H/V curve alon e is not sufficient to uniquely determine the Vs profile , and generally either the layer s’ thicknesses or Vs must be fixed while the other parameter is inverted for. However, this still require s correcting the H/V curve to obtain ellipticity by assuming the energy partition of waves . Another limitation of inverting the H/V curves is that Vs profiles can only be inverted at H/V measurement locations rather than over a uniform grid. This is true because H/V curves at unsampled locations cannot be estimate d; we only estimate f 0,H/V values at unsampled locations rather than the full H/V curve . Thus, this approach will further require interpolating the Vs profiles inverted at each H/V measurement location , which are spatially irregular, to obtain a pseudo -3D Vs model . Future resear ch on this topic will be required. As shown in the accompanying paper, the simple scaling approach used herein provide s Vs models that work well for integrating site -specific spatial variability into GRAs . Information on Case Study Sites The proposed H/V geostatistical approach is validated at t wo downhole array sites , the TIDA and DPDA. In this paper, we only discuss the development of the pseudo- 3D Vs models a t both sites and assess them in light of available geologic information. In a companion paper (Hallal and Cox, 2020) , the pseudo- 3D Vs models are used to account for spatial variability in 1D GRAs and the results are compared to recorded ground motions for validation. The two downhole array sites were selected to represent two distinct geographic and geologic conditions. Importantly, TIDA is believed to be an archetype of 1D conditions, whereas DPDA is the antithesis of such conditions (Tao and Rathje, 2019) . For eac h site, the geology and subsurface properties are described in what follows, along with the H/V measurements . Treasure Island Downhole Array The T IDA site is located in a seismically active region of northern California, approximately 5 km northeast of San Francisco, 12 km west of the Hayward fault, and 20 km east of the San Andreas fault. The downhole array is deployed in 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 simplified stratigraphy for TIDA is shown in Figure 3a along with the depths of the triaxial ground motion accelerometers, as indicated by the black 11 triangle symbols. Existing data indicates that the subsurface near the downhole array site consists of around 15 m of loose sandy hydraulic fill, followed by approximately 15 m of Young San Francisco Bay M ud. The Young Bay M ud (YBM) is underlain by silty sand, Old Bay Mud (OBM), and a mixture of gravel, sand, and clay, all of which are sedimentary deposits of Quaternary age. The bedrock is located at approximately 90 m and is composed of interbedded sandstone, siltstone, and shale of the Franciscan Formation, which is of Jurassic -Cretaceous age (Gibbs et al., 1992; Pass, 1994) . The site has been characterized t hrough several geotechnical testing programs, including V s measurements performed using multiple invasive borehole methods. The V s profiles compared in Figure 3a include measurements from seismic downhole (Gibbs et al., 1992) , shallow crosshole (de Alba et al., 1994) , as well as deep PS suspension logging (Graizer and Shakal, 2004) . The tests were performed in close vicinity to one another, all within an area of approximately 40 m2. The seismic downhole and PS suspension logging Vs profile s indicate a transition from soil to rock at approximately 90 m. However, the PS suspension logging profile shows a significantly higher bedrock V s compared to the seismic downhole measurements. The seismic downhole report mentions that wave arrivals in rock were obscure due to poor energy transmission and wave interference (Gibbs et a l., 1992) . Because of this uncertainty and the limited depth of the crosshole testing , only the PS suspension log is used throughout this study for the TIDA site. Figure 3. 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 profiles at TIDA are from seismic downhole (DH) (Gibbs et al., 1992) , shallow crosshole (CH) (de Alba et al., 1994) , and deep PS suspension logging (Graizer and Shakal, 2004) and the stratigraphy is based on Gibbs et al. (1992) . The measured Vs profiles at DPDA are from inversion of surface wave data (SW) (Nath et al., 1997) and seismic downhole (DH) (Thornley et al., 2019) and the stratigraphy is based on geologic cross sections 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)DH (Gibbs et al., 1992) CH (de Alba et al., 1994) 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)SW (Nath et al., 1997) DH (Thornley et al., 2019) 12 Delaney Park Downhole Array The D PDA site is located in southern Alaska in one of the most seismically active regions in the U.S. Seismicity in this area is influenced by the Castle Mountain fault, Bruin Bay fault, Border Ranges fault, several faults in Cook Inlet, and the Aleutian megathrus t subduction zone (Wen and Kalkan, 2017) . The downhole array is deployed in Delaney Park in downtown Anchorage, which is a lowland area west of the Chugach Mountains. A simplified soil profile for DPDA is shown in Figure 3b (Combellick, 1999) along with the depths of the triaxial ground motion accelerometers . Geologic data indicates that the subsurfac e near the downhole array site consists of around 15 m of glacial outwash (GOW) , followed by approximately 30 m of clay from the Bootlegger Cove Formation, overlying glacial till. The glacial till is located at approximately 47 m and consists of very dense sand and gravel . All these deposits are of Quaternary age and are underlain by thick Tertiary sediment s of sandstone, siltstone, and claystone. The Tertiary deposits in Anchorage are located starting between 100 and 300 m below the surface and are a few kilometers thick (Badal et al., 2004; Combellick, 1999) . Prior to 2019, no V s measurements had been conducted in the immediate vicinity of the DPDA site. Thus, the only V s profile available for GRAs was from the inversion of surface wave data collected by Nath et al. ( 1997) at a site approximately 250 m away. More recently, Thornley et al. (2019) performed seismic downhole measurements in the borehole of the deepest DPDA accelerometer. The Vs profiles b y Nath et al. ( 1997) and Thornley et al. ( 2019) extend approximately 60 m into the Quaternary deposits , as shown in Figure 3b. While the comparison shows differences in the V s over the top 15 m, both profiles indicate a transition from soft soil to hard till ( engineering rock- like material ) at a depth of approximately 47 m. However, only the seismic downhole profile is us ed throughout this study for the DPDA site since Thornley et al. (2019) reported that it provides an improvement in the site response predictions over the previous profile and it is in closer vicinity to the borehole array. H/V Measurements To efficiently characterize the spatial variability over a wide area around the downhole array sites, surface -based seismic measurements were conducted . Specifically, ambient vibrations were recorded using broadband seismometers deployed as part of both surface wave Microtremor Array Measurements (MAM) and single -station deployments; herein w e only discuss the H/V results. Because azimuthal variability in H/V might be of interest in the present study, t he f 0,H/V were obtained following the procedures proposed by Cheng et al. ( 2020) . Using measurements from the two perpendicular horizontal components (i.e., NS and EW) , horizontal motions at six different azimuths spaced at 30 ⁰ were calculated using horizontal polarization analysis . The ambient noise records were divided into many 60- sec time windows and a frequency -domain window -rejection algorithm was used to remove contaminated time windows , as discussed in Cox et al. ( 2020) . Then, for each azimuth and time window, an H/V spectral ratio curve was computed as the ratio of the FAS of the windowed azimuthal horizontal motion to that of the vertical motion and an f 0,H/V value was obtained from the amplitude of the lowest frequency peak. The representative f0,H/V is finally taken as the lognormal -median of the f 0,H/V values of all azimuths and all time windows at the location of the H/V measurement. While this lognormal -median value that accounts for azimuthal variability will be referred to herein as f 0,H/V for simplicity, it is more rigorously notated as LM f0,AZ by Cheng et al. ( 2020) . 13 A map showing the variation in f0,H/V obtained from H /V measurements at 97 locations around the TIDA site is shown in Figure 4a. The f 0,H/V values are quite consistent across most of the island and equal to approximately 0.75 – 0.80 Hz. However, the values increase abruptly in the south where bedrock shallows near the causeway to YBI. Treasure Island was constructed in an estuarine depositional environment, which is generally characterized by a uniform deposition of sediments (Tao and Rathje, 2019) . It follows that if minimal variability exits in the f undamental frequency and in surficial soil layers, then based on Equation 1, the depth to bedrock is also expected to be generally uniform. As discussed earlier, an increase in the fundamental frequency is attributed to an increase in the average Vs and/or a decrease in the dept h to bedrock. We believe that the major contributor to the rapid increase in f 0,H/V is shallower bedrock, which is expected as the rocks of the Franciscan Formation slope upward and ultimately outcrop on YBI . The mapped spatial distribution of f0,H/V inferred from 96 H/V measurements around the DP DA site is shown in Figure 4b. Unlike TIDA, t he distribution of the f 0,H/V values shows considerable heterogeneity in almost all lateral directions. 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. This variability is expected in a glacial depositional environment, in which glacial movements lead to variable geologi c features on the order of tens to hundreds of meters (Tao and Rathje, 2019) . It is expected that lower frequencies correspond to deeper glacial till in the north direction and higher frequencies correspond to s hallower glacial till in the east direction. However, while the depth to the glacial till might be driving the change in f 0,H/V, we also expect variations in the thickness and V s of the surficial layers to play a role. Implementation at Case Study Sites This section presents the implementation of the proposed H/V geostatistical approach at the TIDA and DPDA sites. First, the f0,H/V measurements are kriged to obtain a uniform ly estimated f0,H/V map using the kriging with a trend approac h. Second, the uniformly estimated f 0,H/V map and Figure 4. Fundamenta l 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. 14 the H/V measurement nearest to the location of the measured Vs profile are used to calculate a uniform scaling factor map. Finally, the uniformly calculated scaling factor map and the measured Vs profile ar e used to build a pseudo- 3D Vs model on the same uniform grid by scaling the layer thicknesses of the measured Vs profile . Uniform ly Estimated f0,H/V Maps The first step of the H/V geostatistical approach ( Figure 2) requires uniformly mapping f 0,H/V over an area of interest using the irregularly sampled H/V measurements ; this is achieved herein using the kriging with a trend model approach (Hengl et al., 2004; Isaaks and Srivastava, 1989) . The kriging with a trend model approach is illu strated for the TIDA and DPDA site s in Figures 5 and 6, respectively . Presented are: (a) measured f 0,H/V values, (b) f 0,H/V trend map, (c) f 0,H/V residual map, and (d) the final uniformly estimated f 0,H/V map. The f0,H/V trend m ap at TIDA ( Figure 5b), which is a moving window average of the measurements ( Figure 5a), is comparable to the final uniformly estimated f 0,H/V map ( Figure 5d). Hence, adding the trend map and the residual map (Figure 5c) obtained from kriging on the residuals (i.e., the difference between the measurements and the trend) only causes minimal improvement s in the final f0,H/V map by ± 0.1 Hz at most. This is mainly due to the relatively dense and uniform H/V spatial sampling at TIDA. Figure 5. Uniformly estimated f 0,H/V map at TIDA. Shown on a 50- m grid are: (a) 97 measured f 0,H/V values, (b) f 0,H/V trend obtained from a moving window average of the measurements, (c) uniformly estimated f 0,H/V residual map obtained from kriging on the residuals (shown using circles), which are calculated by subtracting the trend from the f 0,H/V measurements, and (d) the final uniformly estimated f 0,H/V map obtained from summing the trend map (i.e., b) and the kriged residual map (i.e., c). The measured f 0,H/V values are presented in (b) and (d) for comparison. 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) (b)-500 0 500 1000 X (m)-200 0 200 400 600Y (m) (c)-500 0 500 1000 X (m)-200 0 200 400 600Y (m) -0.30 -0.15 0.00 0.15 0.30Residual (Hz) (d)-500 0 500 1000 X (m)-200 0 200 400 600Y (m) 15 On the other han d, due to the more spatially irregular f0,H/V measurements at DPDA, performing kriging on the residuals ( Figure 6c) causes more significant improvement in the estimated map. It is clear that by adding the kriged residual map to the trend, the resulting uniformly estimated f 0,H/V map ( Figure 6d) better matches the measured values compared to the trend map alone ( Figure 6b). This is most evident in the south- east quadrant where the final uniformly estimated map is improved by up to 0.3 Hz . Thus, depending on the site and the spatial distribution of the measurements, different spatial estimation techniques could be used to obtain a uniformly estimated f 0,H/V map. Simple modeling approaches, such as a moving window average of the measurements , could be potentially used when the measurements are relatively dense and Figure 6. Uniformly estimated f 0,H/V map at DPDA. Shown on a 50- m grid are: (a) 96 measured f 0,H/V values, (b) f 0,H/V trend obtained from a moving window average of the measurements, (c) uniformly estimated f 0,H/V residual map obtained from kriging on the residuals (shown using circles), which are calculated by subtracting the trend from the f 0,H/V measurements, and (d) the final uniformly estimated f 0,H/V map obtained from summing the trend map (i.e., b) and the kriged residual map (i.e., c). The measured f 0,H/V values are presented in (b) and (d) for comparison. 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) (b)-800 -400 0 400 800 X (m)-800 -400 0 400 800Y (m) (c)-800 -400 0 400 800 X (m)-800 -400 0 400 800Y (m) -0.30 -0.15 0.00 0.15 0.30Residual (Hz) (d)-800 -400 0 400 800 X (m)-800 -400 0 400 800Y (m) 16 uniform enough to capture changes in f 0,H/V (e.g., TIDA) , whereas more advanced approaches, such as kriging, could be needed when the measurements are sparse and spatially irregular (e.g., DPDA) . However, regardless of the distribution of the measurements, kriging is expected to provide the estimated f 0,H/V map that best fits the measurements since it is a generalized least - squares regression. Uniformly Calculated S f0 Maps After estimating f 0,H/V over a uniform grid, t he second step requires dividing the H/V measurement closest to the location of the measured V s profile (i.e., f0,H/V,Vs ) by the uniformly estimated f0,H/V map to obtain a scaling factor map . The scaling factor maps for the TIDA and DPDA sites are presented in Figures 7 and 8, respectively. Presented are: (a) the uniformly estimated f0,H/V maps , and (b) the uniformly calculated Sf0 maps; the f0,H/V,Vs values used in the calculations are highlighted using magenta in (a) . The f0,H/V,Vs at the TIDA is equal to 0.78 Hz, resulting in S f0 value s ranging between 0.65 and 1.10 (Figure 7b). The estimated f 0,H/V map ( Figure 7a) shows relatively minimal variability in f 0,H/V up to around 500 m away from the downhole array site in all directions , resulting in scaling factor s close to 1.0 near the TIDA boreholes . However, the distribution of S f0 at TIDA appears to be negatively skewed (i.e., skewed to low S f0) due to the higher f 0,H/V values near the causeway to YBI (refer to Figure 4a). Since the layer thicknesses of the measured Vs profile will be multiplied by the negatively skewed scaling factors, then we expect the scaled Vs profiles at TIDA to be also skewed toward shallower bedrock. As for DPDA, the f0,H/V,Vs is equal to 1.24 Hz, r esulting in S f0 values ranging between 0.75 and 1.20 (Figure 8b). The variability in the uniformly estimated f 0,H/V map ( Figure 8a) is relatively significant, resulting in S f0 that deviate to values both greater and less than 1.0 in as little as 300 m away from the DPDA boreholes . However, the distribution of S f0 at DPDA appears to be relatively symmetric , with nearly equal number of cells having scaling factors that deviate by approximately ± 0.3 from 1.0. Thus, we expect the scaled Vs profiles at DPDA to also be relatively symmetric Figure 7. Uniformly calculated S f0 map required to scale the layer thicknesses of the measured Vs profile to match the estimated f 0,H/V map at the TIDA site. Shown on a 50- m grid are: (a) the uniform ly estimated f0,H/V map obtained from the kriging with a trend model approach (i.e., Figure 5d), and (b) the uniformly calculated Sf0 map obtained from dividing the f 0,H/V measurement nearest to the TIDA (also location of Vs measurement) by the estimated f 0,H/V value of each cell (i.e., a). Note that the f 0,H/V measurement nearest to the TIDA is highlighted using magenta 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) (b)-500 0 500 1000 X (m) 0.6 0.8 1.0 1.2 1.4Scaling Factorf0,H/V,Vs 17 about the measured Vs regarding depth to glacial till. This is in contrast to TIDA, where the spatial variability in S f0 is less significant in terms of its magnitude but skewed to low values. Pseudo -3D V s Models To finally build the pseudo- 3D Vs models, Vs profile s are developed on the same uniform grid by simply multiplying the layer thicknesses of the measured Vs profile by the S f0 of each cell . The resulting pseudo- 3D Vs model at the TIDA site is shown in Figure 9. Presented are: (a) the uniformly estimated f 0,H/V map, (b) the scaled Vs profiles shown in 1D, but color coded to match values in the uniformly estimated f 0,H/V map, (c) two 2D Vs cross sections, and (d) a 3D view of the 2D Vs cross sections. The locations of the cross sections are shown in Figure 9a. A total of 646 50-m cells were considered in the analysis at TIDA, resulting in a pseudo- 3D model with approximate surface dimensions of 1.6 km by 1.0 km. As expected, the relatively uniform scaling factors of 1.0 near the TIDA boreholes (Figure 7b), result in scaled V s profiles (Figure 9b) that are mostly comparable to the measured one. It is also clear that the variability is skewed toward a shallower soil- rock interface (i.e., higher f 0,H/V values) near the causeway to YBI in the south, which implies stiff er subsurface conditions ; this is evident in cross section A -A’, which as shown in Figure 9a, extends near the causeway to YBI . The minimal spatial variability close to the TIDA site is consistent with available geologic cross sections , which show that the near subsurface in that area is horizontally stratified (Papadopulos and Eliahu, 2009) . The geologic cross sections , however, do not extend close to YBI. Nevertheless , as discussed earlier, shallower bedrock is expected in the south direction as the rocks of the Franciscan Formation slope upward and ultimately outcrop on YBI . Figure 8. Uniformly calculated S f0 map required to scale the layer thicknesses of the measured Vs profile to match the estimated f 0,H/V map at the DPDA site. Shown on a 50- m grid are: (a) the uniformly estimated f0,H/V map obtained from the kriging with a trend model approach (i.e., Figure 6d), and (b) the uniformly calculated Sf0 map obtained from dividing the f 0,H/V measurement nearest to the DPDA (also location of Vs measurement) by the estimated f 0,H/V value of each cell (i.e., a). Note that the f 0,H/V measurement nearest to the DPDA is highlighted using magenta 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) (b)-800 -400 0 400 800 X (m) 0.6 0.8 1.0 1.2 1.4Scaling Factorf0,H/V,Vs 18 The pseudo- 3D Vs model at the DP DA site is demonstrated in Figure 10. 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 ( Figure 10b) mirror many features observed in the f0,H/V map (Figure 10a) and Sf0 map ( Figure 8b). As expected, the highly variable f0,H/V and Sf0 maps near the DPDA boreholes result in many scaled V s profiles with both shallower or deeper glacial till. This variabil ity, while significant, is relatively symmetric about the measured Vs profile. The lower f0,H/V values in the north result in scaling factors above 1.0, and consequently, deeper glacial till . This is evident in the scaled Vs profiles shown in dark red (Figure 10b) and in cross section B -B’. Conversely, t he higher f0,H/V values in the south -east quadrant result in scaling factors below 1.0, and consequently, shallower glacial till, which is evident in th e scaled Vs profiles shown in bl ue (Figure 10b) and in cross section A -A’. The pattern of the Vs variability is in g ood agreement with several geologic cross sections around the site , which show deeper and shallower glacial till in the north and east directions, respectively (Combel lick, 1999) . Figure 9. Pseudo- 3D Vs model developed using the H/V geostatistical approach at the TI DA site. Shown are: (a) the uniformly estimated f 0,H/V map, (b) 1D representation of the spatially variable V S profiles calculated from the measured profile and the scaling factors (Figure 7b), (c) 2D Vs cross sections A -A’ and B-B’, and (d) 3D fence diagram of the 2D Vs cross sections. Note that the scaled 1D Vs profiles in (b) are colored by the f 0,H/V values of their corresponding cells shown in (a). The borehole array and triaxial accelerometer sensor depths are shown in black and white, respectively, in both (c) and (d). 19 Discussion We have observed that the TIDA site has low variability that is skewed to shallower bedrock toward the south, near the causeway to YBI , whereas the DPDA site has high Vs variability but which is relatively symmetric about the measured Vs profile near the downhole array site . Such variability could not be represented by purely stochastic models based on generic parameters derived for other sites or simply assumed. This underscore s the importance of building Vs models based on in-situ measurements to ensure proper incorporation of site -specific spatial variability in GRAs. The aforementioned evaluation of the site -specific variability at each site was qualitative and herein, we present quantitative measures of the spatial variability at each site. We believe that the Vs models should reflect both the magnitude and skewness of the spatial variability . The magnitude is the extent of the variability a nd can be quantified using the standard deviation. The Figure 10. Pseudo- 3D Vs model developed using the H/V geostatistical approach at the DPDA site. Shown are: (a) the uniformly estimated f 0,H/V map, (b) 1D representation of the spatially variable V S profiles calculated from the measured profile and the scaling factors (Figur e 8b), (c) 2D Vs cross sections A -A’ and B-B’, and (d) 3D fence diagram of the 2D Vs cross sections. Note that the scaled 1D Vs profiles in (b) are colored by the f 0,H/V values of their corresponding cells shown in (a). The borehole array and triaxial accelerometer sensor depths are shown in black and white, respectively, in both (c) and (d). 20 skewness should indicate whether variability is favored in a certain direction, such as having mostly stiffer or softer conditions in a certain area or direction around the site , and is quantified using the sample skewness . The descriptive statistics of the site -specific spatial variability (i.e., standard deviation and skewness) are evaluated herein for the uniformly calculated Sf0 values. We choose S f0 instead of f0,H/V because the former values are normalized ; this ensures that the difference in magnitude of f0,H/V between the sites does not affect the results. While the descriptive statistics can be calculated at each site over the entire grid, we present the results as a function of several discrete incorporated areas . This allows us to compare the variability at both sites across similar incorporated areas and eventually provide insights on the area likely influencing site response. We incrementally increase the area surrounding the downhole array site by one cell (i.e., by 50 m) in each of the four directions (north, south, east, and west) and calculate the standard deviation and skewness for each incorporated area. For TIDA, in which the geometry of the map is asymme tric, we keep on increasing the incorporated area in the directions that can be extended, until covering the full map. Magnitude of Spatial Variability We calculate the lognormal standard deviation in S f0 (σlnSf0) since lognormal distributions are common in geotechnical earthquake applications ; this also allows for a seamless transition between the standard deviation of S f0, f0,H/V, and that of depth to rock (z rock), as discussed later in this section . The σlnSf0 comparison is presented in Figure 11a. Additionally, for clarity, we highlight different areas in Figure 11a; those areas are outlined in the uniformly calculated Sf0 maps at TIDA and DPDA in Figures 11c and 11d, respectively. While the sites have comparable variability at small incorporated areas , DPDA has an increased σ lnSf0 beyond a n incorporated area of 400 m x 400 m (i.e., the area highlighted in teal ; 0.16 km2). This observation is in line with Tao and Rathje (2019) , who hypothesized that the DPDA site is more variable than TIDA based on geology and results from 1D GRAs . However, in this paper, we additionally highlight the need to identify an area of interest , as DPDA is more variable only after a certain area is incorporated. This issue is addressed in Hallal and Cox ( 2020) , in which the effects of the incorporated area on the site response predictions are examined. The evaluation of the standard deviation of S f0 in logarithmic space results in an interesting correspondence between σlnSf0 and that of f 0,H/V (σlnf0,H/V) and depth to rock ( σlnZrock ). The lognormal standard deviation depth to rock is commonly used as an input parameter in Vs randomization approaches to estimate the variation in the depth to bedrock. Following Equation 2 of the H/V geostatistical approach, and by replacing the thickness of the soil column above bedrock (H) by the depth to bedrock (zrock), the relation shown in Equation 3 can be obtained. 𝑆𝑆𝑓𝑓0=𝑧𝑧𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟 𝑧𝑧𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟 ,𝑉𝑉𝑉𝑉→𝑙𝑙𝑙𝑙(𝑧𝑧𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟)=𝑙𝑙𝑙𝑙(𝑧𝑧𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟 ,𝑉𝑉𝑉𝑉)+𝑙𝑙𝑙𝑙(𝑆𝑆𝑓𝑓0) (3) Then, by evaluating the variance of the left - and right -hand side of Equation 3, we can easily show that the variance in z rock is equal to that in S f0 in logarithmic space, as shown in Equation 4. This is true because the bedrock depth of the measured Vs profile (z rock,Vs) is a constant and does not affect the value of the variance. Hence, by taking the square root of the variance, we can prove 21 that σlnSf0 = σlnZrock . The s ame logic can be used to show that σ lnSf0 = σlnf0,H/V. Thus, without having to calculate scaling factors and just by calculating σlnf0,H/V, engineers can obtain an estimate of the site-specific σlnZrock and consequently, use better -informed parameters when applying Vs randomization for GRAs. However, this is only true under the assumptions of the H/V geostatistical approach, wherein the layer thicknesses and depth to the major impedance contrast are assumed to be primarily controlling the change in f 0,H/V. 𝑉𝑉𝑉𝑉𝑉𝑉 [𝑙𝑙𝑙𝑙(𝑧𝑧𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟)]=𝑉𝑉𝑉𝑉𝑉𝑉�𝑙𝑙𝑙𝑙�𝑧𝑧𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟 ,𝑉𝑉𝑉𝑉�+𝑙𝑙𝑙𝑙�𝑆𝑆𝑓𝑓0��=𝑉𝑉𝑉𝑉𝑉𝑉�𝑙𝑙𝑙𝑙�𝑆𝑆𝑓𝑓0��→𝜎𝜎𝑙𝑙𝑙𝑙𝑙𝑙𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟=𝜎𝜎𝑙𝑙𝑙𝑙𝑆𝑆𝑓𝑓0 (4) Figure 11. Effect of varying the incorporated spatial area in the uniformly calculated S f0 map on: (a) the lognormal standard deviation in S f0 (σlnSf0), and (b) the sample skewness (G 1) of S f0. (c & d) The uniformly calculated S f0 map at TIDA and DPDA, respectively, showing the selected areas whose results are highlighted in (a) and (b) using t he same colors. Note that σ lnSf0 is equal to the standard deviation of f 0,H/V in logarithmic space ( σlnf0,H/V). 0 0.05 0.10 0.15lnSf0 (a)-500 0 500 1000 X (m)-200 0 200 400 600Y (m) TIDA (c) 0.6 0.8 1.0 1.2 1.4Scaling Factor (b)0 0.5 1 1.5 2 2.5 Area (km2)-2 -1 0 1 2G1 -800 -400 0 400 800 X (m)-800 -400 0 400 800Y (m) DPDA (d) 0.6 0.8 1.0 1.2 1.4Scaling FactorTIDA DPDA Highly Skewed Moderately Skewed Symmetric 22 Skewness of Spatial Variability Another important descriptive statistic of the site -specific spatial variability is its skewness. While the standard deviation quantifies the extent of the variability, it does not provide any information on the direction of this varia bility , which as shown in Hallal and Cox ( 2020) , affects the site response . To quantify this directio n, we use the sample skewness (G 1) as shown in Equation 5 (Joanes and Gill, 1998) . In this equation, n is the total numbe r of cells in the incorporated area. 𝐺𝐺1=�𝑙𝑙(𝑙𝑙−1) 𝑙𝑙−21 𝑙𝑙∑�𝑆𝑆𝑓𝑓0,𝑖𝑖−𝑆𝑆𝑓𝑓0�����3𝑙𝑙 𝑖𝑖=1 ��1 𝑙𝑙∑�𝑆𝑆𝑓𝑓0,𝑖𝑖−𝑆𝑆𝑓𝑓0�����2𝑙𝑙 𝑖𝑖=1 �3 (5) The variation of the skewness as a function of the incorporated area at both sites is presented in Figure 11b. Generally, a G 1 value of zero indicates a symmetric distribution, whereas a negative value implies skewness to low values, and vice versa. We observe that the distribution of S f0 at TIDA is approximately symmetric up to an area of approximately 1,400 m x 1,000 m (i.e., the area highli ghted in dark green; approximately 1.4 km2), beyond which it starts to become highly negatively skewed, whereas DPDA is mostly negatively skewed but gradually becomes less skewed beyond an area of 800 m x 800 m (i.e., the area highlighted in light pink ; 0.64 km2). Based on the G 1 boundaries set by Bulmer (1979) for symmetric, moderately s kewed, and highly skewed distributions , which are highlighted in Figure 11b using green , yellow , and red shading, respectively, the distribution of S f0 when incorporating the full area at TIDA is considered to be highly skewed (G 1 < -1) whereas that at DPD A is at the boundary between symmetric and moderately skewed (G1 ≈ -0.5). Again, the change of the skewness with area illustrates the importance of specifying an area of interest over which the spatial variability needs to be modeled. It also highlights th e need to figure out the spatial area that impacts site response at any given location. The areas influencing recorded site response at TIDA and DPDA are investigated in the companion paper by Hallal and Cox ( 2020) . Conclusions Accounting for spatial variability in subsurface conditions is of paramount importance for seismic analysis , design, and risk assessment. While this paper is not the first to recognize the need to implement an approach that could be tailored to account for site-specific spatial variability in GRAs , no study we are aware of attempts to do so usi ng simple, cost -effective site measurements that could be feasibly obtaine d. This paper proposes a method for building a pseudo- 3D Vs model based on a geo statistical framework of H/V measurements that can be used to incorporate spatial variability into 1D, 2D, and 3D GRAs . The new framework utilizes easily attained measurements that provide valuable information about the variation of properties known to str ongly impact site - specific wave propagation, does not require the assumption of a probabilistic model to map the spatial variability with tens to hundreds of realizations , and enables the construction of a realistic pseudo- 3D Vs model . 23 The scaled Vs profil es built using the H/V geostatistical approach are expected to be a more realistic representation of the actual subsurface conditions than purely stochastic models , which generally use input parameters derived for other sites or simply assumed. The rationale behind this argument is that a unique spatial variability exists for each site, meaning that calculated parameters based on relationships derived for other sites in the literatur e may not represent the actual site properti es. Additionally, the subsurface properties are spatially correlated rather than completely random, and this spatial correlation also varies from site to site . The H/V approach addresses both of these issues by integrating site -specific H/V measurements, w hich are dependent on subsurface conditions, with geostatistical tools that can model the site -specific spatial correlation of the subsurface properties. Additionally, the approach utilizes the experimental site signature, which allows for V s uncertainty t o be accounted for in a rational manner. While the implementation of the H/V geostatistical approach at the two case study sites presented herein was each based on approximately 100 H /V measurements, the approach could certainly be implemented with fewer m easurement locations. Even so, the H/V data was collected easily at each site in less than two f ull workdays by crews of 4 -5 people using 10 seismometers. This is less than the amount of time dedicated to typical geotechnical site characterizations on even routine project s. Furthermore, this testing could have been performed faster, as m any of the H/V measurements were clustered in arrays as part of surface wave MAM testing , which are not needed for implementing the H/V approach. For future applications , the H/V measurements should be relatively uniformly distributed across the area of interest, with good azimuthal coverage where possible . The proposed methodology for building psuedo- 3D Vs models aims to achieve the ultimate goal of improving the modeling of site response in engineering practice. The approach attempts to account for realistic Vs spatial variability in a meaningful way using site -specific measurements that are relevant to site response. Furthermore, as the same pseudo -3D Vs models can be incorporated in 1D, 2D, and 3D GRAs, this allows for future research on the appropriate /necessary dimensionality of GRAs at different sites . In addition, the pseudo- 3D Vs model can provide valuable information on the subsurface structure for applications such as estimating the depth to significant Vs horizons and lognormal standard deviation on depth to rock for use in ground motion prediction equations and Vs randomization. While t he accuracy of th e generated Vs models was evaluated only in comparison to geologic information, a further validation and comparison of the predicted and observed ground motions at TIDA and DPDA using 1D GRAs that incorporate the pseudo- 3D Vs models are discussed in the ac companying paper by Hallal and Cox ( 2020) . Acknowledgements The H/V data fro m TIDA and DPDA were collected with funding from Pacific Gas and Electric (PG&E). However, any opinions, findings, conclusions, or recommendations expre ssed 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 collecting the H/V data at TIDA and DPDA, and Dr. Albert Kottke, D r. David Teague, and Dr. Krishna Kumar for their help in collecting the H/V data at TIDA. 24