## 1 Introduction  
The accurate forecasting of acoustics wave propagation and scattering in shallow water environments is important in naval operations and underwater communications. A large number of theoretical studies have been carried out to model the sound scattering effects in complex underwater environments in three dimensions (3D). The 3D coupled mode method could produce efficient and accurate prediction of the sound field for specified geometries (e.g., 3D axisymmetric geometry [1, 2], longitudinal invariant geometry [3]). The 3D ray model has been proposed for high frequency sound propagation ([4]). The 3D parabolic equation (PE) method is another powerful tool for the prediction of 3D sound propagation (e.g., [5]). Although several newly developed PE models, including two-wayscattering, have been developed, generally for two-dimensional (2D) cases [6; 7] and for 3D cases formulated in cylindrical coordinates [8], PE is mostly used for the prediction of one-way sound propagation as backward propagation is generally neglected in this method. For direct numerical simulations (i.e., discretization methods without any assumptions), despite the continued development of computing capabilities, it is a challenge to perform 3D simulations of ocean acoustics in realistic environments due to the high computational costs. Volume discretization methods, such as the finite element method (FEM), finite volume method (FVM), and finite difference method (FDM), require fine meshes over the entire volume of the computational domain. Due to the large computational domain (\(\sim\)O(10) km) and the small acoustic wavelength (\(\sim\)O(1-10) m), the 3D mesh entails a large computational cost. By only discretizing the domain, the boundary discretization method or boundary element method (BEM) has a much lower computational cost than the volume discretization methods. By discretizing the domain boundaries and interfaces into boundary elements, BEM numerically formulate the acoustics propagation problem as a linear equation system in the form of [_A][x] = [b]_. Matrix [_A_] is an \(N\times N\) dense matrix composed of influence coefficients. Normally, this linear equation system requires O(\(N^{2}\)) operations to be constructed and O(\(N^{2}\))/ O(\(N^{3}\)) to be solved by iterative/direct solvers. For applications in real shallow water acoustics, N is often in the order of O(10\({}^{6-8}\)). This induces a large computational cost and limits the use of BEM in real shallow water acoustics. Another major drawback of BEM is that classical BEM can only formulate the linear equation system for underwater medium with homogeneous properties (e.g., density and sound speed). However, complex shallow water environments with density and sound speed variations cannot be represented by a single homogeneous fluid layer. In this study, to overcome the computational cost difficulty in BEM, we used the pre-corrected fast Fourier transform (PFFT) approach to accelerate the BEM solver for the calculation of acoustics propagation in a shallow water environment. The pre-corrected fast Fourier transform (PFFT) approach has been used to increase the speed of the BEM solver from (\(N^{2}\)) to (\(N\) log \(N\)) in different areas (e.g., [9; 10; 11]). By using PFFT, we can reduce the computational cost in the order of O(10\({}^{6-8}\)). Compared with existing studies using PFFT for underwater acoustics problems (i.e., [11]), the present study made the following contributions. First, to overcome drawbacks due to homogenous assumptions, we simulated the medium inhomogeneity by coupling of multiple homogeneous sub-layers at their interfaces. Through a numerical test, we checked the convergence of the multilayer PFFT-BEM (Version 2019, Massachusetts Institute of Technology, Cambridge, MA, USA, 2019) with respect to the number of sub-layers to validate the accuracy of the multi-layer formulation in our method.