Subsurface spatial variability is known to significantly influence the frequency content and amplitude of seismic ground shaking. A significant amount of seismic site response research over the past decade has focused on our abilities to replicate recorded ground motions at borehole array sites, where both the input (rock) and output (surface) ground motions are known. When viewed in aggregate, these studies have found that approximately 50% of borehole array sites are poorly modeled using one-dimensional (1D) ground response analyses (GRAs) based on a single shear wave velocity (Vs) profile, with individual studies reporting values between approximately 30-80%. When 1D GRAs fail to accurately predict recorded site response, the site is often considered too complex to be effectively modeled as 1D. While three-dimensional (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 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 guessed, input parameters, introducing significant uncertainties into the site response predictions. This research 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. The proposed approach distinguishes itself from previous studies in three key ways: (1) it requires only a single, accurately measured Vs profile down to engineering bedrock, (2) it relies majorly on estimates of fundamental site frequency (f₀; a key parameter governing site effects) obtained from simple horizontal-to-vertical spectral ratio (H/V) noise measurements (f₀,[subscript 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. At the two sites investigated in this research, the H/V geostatistical approach is capable of generating pseudo-3D Vs models that reliably capture important subsurface features present in geologic cross-sections. Furthermore, the 1D GRA predictions associated with the H/V geostatistical approach were more accurate than those associated with common and recently proposed strategies of accounting for Vs variability. One of the most significant contributions of this research is providing insights on the lateral area influencing seismic site response. The H/V geostatistical approach enables predicting site response as a function of the spatial variability across different footprints. The results show that 1D GRAs are significantly improved when an area of at least 400 m x 400 m (i.e., 0.16 km²) is incorporated, and even larger incorporated areas could produce better results. Thus, this size of an area might be considered as a minimum over which to account for spatial variability in GRAs. These results are supported by two-dimensional (2D) GRAs, which show that incorporating variability along at least 600 m was needed to appropriately model decreased amplification at the fundamental mode caused by wave scattering, while a lateral extent of 1700 m was needed to more accurately model other observed complex phenomena. These results and insights work toward achieving more accurate and reliable seismic hazard assessment and risk mitigation.