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 of 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 describes 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 (f0,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 f0,H/V value, thereby producing a pseudo-3D Vs model. This approach is demonstrated at the Treasure Island and Delaney Park Downhole Array sites (TIDA and DPDA, respectively). While the pseudo-3D Vs models can be used to incorporate spatial variability into 1D, 2D, or 3D GRAs, their implementation in 1D GRAs at TIDA and DPDA is discussed in a companion paper.