Abstract: To fully consider the horizontal coupling effect on sound propagation in a three dimensional waveguide with varying topography, a three-dimensional higher-order fluid parabolic equation model is developed in the cylindrical coordinate. In the model, the two-dimensional square operator is split into the one-dimensional square operators by Taylor expansion, and the one-dimensional square operators are replaced with the product of rational fractions of differential operators using a split-step higher-order Pad@method, as a result, the differential equation is represented with a matrix equation through the Galerkin discretization method. Besides, the energy-conserving assumption is used to deal with the ocean bottom boundaries to consider the influence of complicated ocean topographies on sound propagation. The recursion calculation of three-dimensional sound field is achieved using an alternating direction implicit format. Some sound propagation calculations of typical ocean topographies are presented, such as wedges and seamounts. Results show that the three-dimensional higher-order parabolic equation model in cylindrical coordinates can achieve higher accuracies of sound field calculations under the cases of wedges and seamounts, comparing with the existing models. The all-space field calculations are realized by the developed model.