Quantitative modeling of bottom‐interacting ocean acoustic waves is complicated by the long propagation ranges and by the complexity of the scattering targets. We employ a two‐dimensional (2‐D) hybrid technique combining Gaussian beam, finite difference, and Kirchhoff integral solutions of the wave equation to simulate ocean acoustic experiments within half of a convergence zone in the SOFAR channel. The 2‐D modeling approach is reasonable due to the one‐dimensional (1‐D) velocity distribution in the water column and the strong lineation of the seafloor morphology parallel to the mid‐ocean ridges. Full‐waveform modeling of ocean acoustic data requires that the topography and the material properties of the seafloor are available at scales that are several orders of magnitude smaller than typical bathymetric sampling rates. We have therefore investigated the effects on the ocean acoustic response of a stochastic interpolation scheme used to generate seafloor models. For typical grazing angles of the incident wave field (approximately 5°–20°), we found that different stochastic realizations of the same seafloor segment (sampled at 200 m) yield an intrinsic uncertainty of the order of 3–8 dB in amplitude and 0.1–0.3 s in time for individual prominent events in the reverberant acoustic field. Hybrid simulations are compared to beam‐formed ocean acoustic data collected during the Acoustic Reverberation Special Research Program (ARSRP) cruises. Side lobe noise in the observed acoustic data is simulated by adding band‐limited white noise at −30 dB relative to the maximum intensity in the synthetic data. Numerical simulations can be limited to the response of only one of the mirror azimuth beams provided that the experimental geometry is suitably chosen. For the 2‐D approximation to be valid, the cross‐range resolution of the observed data must be smaller than the characteristic scale of seafloor lineations, and the beams of interest must be approximately perpendicular to the dominant structural grain. Under these premises, the arrival time and maximum intensity of the observed backscattered acoustic events can be modeled within 0.3 s and 5–10 dB, respectively. The mismatch in arrival time is interpreted to be due predominantly to intrinsic uncertainties in the stochastic interpolation of the seafloor profiles, whereas the fact that the intensity of the backscattered events is systematically overestimated is attributed to 3‐D effects within the “footprint” of the beam and/or to underestimated noise levels.