A three-dimensional (3D) finite difference (FD) model with formal fourth-order accuracy has been developed for the ocean acoustic Helmholtz equation (HE), which can be used to address arbitrary bathymetry and provide more accurate benchmark solutions for other 3D underwater acoustic approximate models. The derivatives in the acoustic HE are numerically discretized based on regular grids, and the perfectly matched layer is introduced to absorb unphysical reflections from the boundaries where Sommerfeld radiation conditions are deployed. The system of linear equations is solved using a parallel matrix-free geometric multigrid preconditioned biconjugate gradient stabilized iteration method, and the code (named COACH) is run on the Tianhe-2 supercomputer in China. Four 3D topographic benchmark acoustic cases-a wedge waveguide, Gaussian canyon, conical seamount, and corrugated seabed-are simulated to test the present FD model, and the maximum number of grid points reaches 33.15 × 109 in the wedge waveguide case, running in parallel with 988 central processing unit cores. Furthermore, the accuracy and generality of the present model have been verified by solution comparisons with other available 3D acoustic propagation models, and the two-dimensional and 3D transmission loss contours are presented to facilitate the distinguishing among the acoustic field features of these cases.