Numerical modeling techniques commonly used to compute the response of soil and rock media under earthquake shaking are evaluated by analyzing the observations provided by instrumented borehole arrays. The NIED Kik-Net database in Japan is selected as the main source of borehole array data for this study. The stiffness of the site and the availability of high intensity motions are the primary factors considered towards the selection of appropriate Kik-Net borehole arrays for investigation. Overall, 13 instrumented vertical arrays are investigated using over 750 recorded ground motions characterized by low (less than 0.05 g) to high (greater than 0.3 g) recorded peak ground accelerations at the downhole sensor. Based on data from the selected borehole arrays, site response predictions using 1-D linear elastic (LE) analysis, equivalent linear (EQL) analysis, equivalent linear analysis with frequency-dependent soil properties (EQL-FD), and fully nonlinear analysis (NL) are compared with the borehole observations. Initially, the low intensity motions are used to evaluate common assumptions regarding 1-D site response analysis. First, we identify the borehole wavefield best simulating the actual boundary condition at depth by comparing the theoretical linear-elastic (LE) and observed responses. Then, we identify the best-fit small-strain damping profiles that can incorporate the additional in-situ attenuation mechanisms. Finally, we assess the validity of the one-dimensional modeling assumption. Our analyses indicate that the appropriate boundary condition for analysis of a borehole array depends on the depth of the borehole sensor and that, for most of the considered vertical arrays, the one-dimensional assumption reasonably simulates the actual wave propagation pattern. In the second part of this study, we evaluate the accuracy of the EQL, EQL-FD and NL site response methods by quantifying the misfit (i.e., residual) between the simulations and observations at different levels of shaking. The evaluation of the performance of the theoretical models is made both on a site-by-site basis and in an aggregated manner. Thereafter, the variability in the predicted response from the three site response methods is assessed. Comparisons with the observed responses indicate that the misfit of simulations can be significant at short periods and large strains. Moreover, all models seem to be characterized by the same level of variability irrespectively of the level of shaking. Finally, several procedures that can be used to improve the accuracy of the one-dimensional EQL, EQL-FD and NL site response analyses, are investigated. First, an attempt to take into account the shear strength of the soil materials at large shear strains is made. Additionally, several modifications to the EQL-FD approach are proposed. The proposed modifications are evaluated against recordings from the borehole arrays. Our analyses indicate that the accuracy of the theoretical models can be, partly, increased by incorporating the proposed modifications. / text