Executive Summary Along with source and path effects, site response is an essential component of ground motion prediction. Widely used ground motion models (GMMs), also known as ground motion prediction equations, provide an ergodic representation of each component in the sense that observations from global databases during the observation period (generally the last few decades) are taken to apply for a particular site and tectonic setting of interest, following conditioning on relevant parameters (magnitude, distance, time-averaged shear wave velocity in the upper 30 m, VS30). Such models inherently average across effects that may exhibit location-to-location variability, increasing model dispersion. The use of non-ergodic site response has gained increasing attention in recent years as a means by which to increase model accuracy and reduce model dispersion, both of which affect the outcomes of seismic hazard analysis. The analysis of non-ergodic site response can, in general, be undertaken through analysis of recordings at the site of interest, or (in the absence of such data) through the use of geotechnical simulations. The most common simulation approach, known as ground response analyses, simplifies the actual site response problem by assuming horizontal soil layers and vertically propagating waves. The objective of this research was to compile and analyze data from vertical arrays in California for the purpose of evaluating ground response analysis as a method of predicting non-ergodic site response and to estimate epistemic uncertainties associated with its application. More specifically, we investigated three questions: (1) how effective is ground response analysis at predicting observed small strain, essentially visco-elastic, site response as observed at California vertical arrays?; (2) which models for small strain damping are most effective for use in ground response analyses?; and (3) recognizing the imperfect ability of ground response analysis to capture observed site response effects, how should epistemic uncertainty in site response be represented when it is estimated using ground response analysis procedures? We consider a database of 21 California vertical arrays operated by the California Strong Motion Instrumentation Program (CSMIP) and the University of California Santa Barbara. Each of the considered arrays has ≥ four surface and downhole ground motion recordings, and cumulatively our database contains 287 ground motion pairs from 207 earthquakes. Uncorrected (version 1) acceleration time series were processed using standard procedures developed for the Next Generation Attenuation projects. Although this database is considerably smaller than the KiK-net database that has been widely used in prior research, it has two notable benefits: (1) it represents site response for a distinct region (California) with a different geologic history and (2) the available velocity profile data is of higher resolution and quality, and is mostly accompanied by geotechnical logs with detailed information on soil conditions. The processed data were plotted as surface-to-downhole transfer functions and ratios of 5% damped pseudo acceleration response spectra, each of which represents in different ways the frequency-dependent site response for essentially visco-elastic conditions over the depth range of the arrays. The site response is considered visco-elastic because relatively weak ground motions were selected for analysis. Ground response analyses were performed using the measured shear wave velocities, various damping models, and the recorded base motion as input. We find a higher percentage of California sites, as compared to KiK-net sites from Japan, to have a reasonable match of empirical and theoretical transfer function shapes. The empirical transfer functions also have a greater degree of event-to-event consistency than has been found previously in Japan. We were unsuccessful at diagnosing conditions that would indicate, a priori, whether ground response analyses are or are not effective for a particulate site. Three damping models were considered in the ground response analyses – geotechnical models, models for quality factor (Q) based on seismological inversion, and models derived from the site-specific site diminutive parameter (κ0). These models represent, to varying degrees, the attenuation of ground motions from two physical mechanisms – soil intrinsic (hysteretic) damping and wave scattering, both of which would be expected to be present to varying degrees at a given site. Despite the different mechanisms, the principle means by which to incorporate damping in ground response analysis is through the soil hysteretic damping considered in the analysis (D), which was the approach taken here. As expected, the effects of using different damping models are concentrated at high frequencies, specifically those higher that the frequency of the modelled soil column. Ground response analyses based on geotechnical models underestimate site attenuation, which has been observed previously and is expected because scattering effects are neglected. The models based on seismological inversion tend to overestimate site attenuation; this conclusion is likely not fully general, but applies to the considered data inventory. We describe a means by which to adjust geotechnical models for D using observations of κ0, more specifically the change of κ0 across the depth range of vertical arrays. This approach yielded intermediate levels of site attenuation that modestly improved prediction accuracy and reduced the dispersion of residuals relative to the other damping models. We use the residuals of ground response analysis predictions of site response, relative to observation, to quantify epistemic uncertainties. Our proposed methodology partitions prediction residuals into between- and within-site components, and takes the between-site standard deviation as a quantification of epistemic uncertainty. Our results suggest values ranging from 0.35-0.5 in natural log units, which is surprisingly consistent with related prior results from other investigators using Japanese data. We also find levels of event-to-event variability for a given site that are consistent with observations elsewhere, including Japan and Taiwan.