## 4 Applications of 3D PFFT-BEM
### Model Comparisons of the PFFT-BEM with the 3D Parabolic Equation Method
#### 4.1.1 Formulation and Assumptions of the PE Method  
In order to demonstrate the comparison between PFFT-BEM and 3D PE, we briefly describe the formulation of the PE method following previous studies (e.g., [14; 19]). We begin with the Helmholtz equation in the frequency domain to determine the acoustic wave pressure \(p(\mathbf{r})\) in a cylindrical coordinate system:  
\[\
abla p+{k_{0}}^{2}n^{2}p=0, \tag{27}\]  
where \(k_{0}=\omega/c_{0}\) is the reference wavenumber, \(c_{0}\) is the reference sound speed (defined based on the sound profiles [14]), \(\omega\) is the sound frequency, and \(n=c_{0}/c=k/k_{0}\) is the index of refraction. The derivation of PE methods in a Cartesian coordinate system is similar to that one in a cylindrical coordinate system  
Figure 10: (**a**) Shallow water summer sound profile with a gradient region at 20–40 m; (**b**) transmission loss (TL) in the waveguide with summer sound profile: comparison among PFFT-BEM results with different number of layers (\(N\)\({}_{i}\)). and is therefore omitted here. By applying the separation of variables and accounting for only outgoing propagation sound waves, we can approximate \(p(r)\) by  
\[p(r)=H_{0}(k_{0}r)\varphi(r,\theta,z,\omega), \tag{28}\]  
where \(H_{0}\) is the zeroth-order Hankel function of the first kind. If we assume that there is only outgoing propagation in the range, the Helmholtz equation reduces to an initial boundary value problem governing by a parabolic equation as  
\[\frac{\partial\varphi}{\partial r}=ik_{0}\Big{(}\sqrt{I+X+Y}\Big{)}\varphi, \tag{29}\]  
with \(\varphi(r=r_{0},\theta,z)=\varphi^{0}(\theta,z)\), where, \(\varphi^{0}\) denotes the initial outgoing field at \(r=r_{0}\), \(I\) is the identity operator, \(X\) is the 2D depth operator in the \(r\)-\(z\) plane, and \(Y\) is the azimuthal operator, defined as  
\[X=\Big{(}n^{2}-1\Big{)}I+\frac{1}{k_{0}}\frac{\partial}{\partial z}\Bigg{(} \frac{1}{\rho}\frac{\partial}{\partial z}\Bigg{)}, \tag{30}\]  
\[Y=\frac{1}{\big{(}k_{0}r\big{)}^{2}}\frac{\partial^{2}}{\partial\theta^{2}}. \tag{31}\]  
Solving Equation (29) usually has a large computational cost. Several methods approximate the nonlinear term (square root operator in Equation (29)) to reduce the computational cost. The 3D PE proposed by [PERSON] [5] expands the square root operator and only retain the linear terms as  
\[\sqrt{I+X+Y}=I+\frac{1}{2}(X+Y)+O\Big{(}X^{2},Y^{2},XY\Big{)}. \tag{32}\]  
It is known that this formulation allows only for narrow angle propagation along \(\theta\). For wide angle approximation, [PERSON] [20] approximates the square root operator using the Pade approximation along \(z\) and a linear approximation long \(\theta\) as  
\[\sqrt{I+X+Y}=I+\sum_{k=1}^{n_{p}}\frac{a_{k,n_{p}}X}{I+b_{k,n_{p}}X}+\frac{1}{2 }Y+O\Big{(}X^{2},Y^{2},XY\Big{)}, \tag{33}\]  
where \(n_{p}\) is the number of terms, and \(a_{k,n_{p}}\), \(b_{k,n_{p}}\), \(1\leq k\leq n_{p}\) are the real coefficients given by  
\[a_{k,n_{p}}=\frac{2}{2n_{p}+1}sin^{2}\Bigg{(}\frac{k\pi}{2n_{p}+1}\Bigg{)}, \tag{34}\]  
\[b_{k,n_{p}}=cos^{2}\Bigg{(}\frac{k\pi}{2n_{p}+1}\Bigg{)}. \tag{35}\]  
Another approximation form, which includes part of the nonlinear term, was proposed by [PERSON] and [PERSON] [21] in a cylindrical coordinate 3D PE method as  
Figure 10:0  
The numerical results calculated based on these approximations usually have a phase error in the far field [20]. These errors are quantified in the present study by comparison between direct simulation results from PFFT-BEM and results from 3D PE. Since PE only has one range derivative in the marching direction, it is not possible to enforce both boundary conditions across the vertical boundaries of the marching directions [14]. For most 3D PEmethods, only continuity of pressure is enforced across these vertical boundaries. To conserve energy at each marching step, [PERSON] [22] introduced a variable \(\overline{p}=p/\sqrt{\rho c}\), where \(c\) is the sound speed and \(\rho\) is the medium density.