A method for numerically simulating hydrophone vertical seismic profiles (VSP) and crosswell data measured in a fluid-filled borehole (either open or cased) embedded in stratified media is presented. The method makes use of both the borehole coupling theory and the global matrix formulation for computing synthetic seismograms in a stratified medium. The global matrix formulation is used to calculate the stress field at the borehole location. Borehole coupling theory is then employed to obtain the pressure in the borehole fluid. Comparisons with exact solutions for an open borehole in a homogeneous and unbounded formation show that this method is accurate for frequencies below 2 kHz. This method is used to model the Kent Cliffs hydrophone VSP data, where good agreement between the numerical simulations and the field measurements has been found, in both traveltimes and rms amplitudes of the direct P-wave. Examples show that this method is efficient and accurate, and can be applied to model VSP and crosswell experiments using an array of hydrophones.