## 3 Validation of Multi-Layer PFFT-BEM Solutions  
In this section, we first demonstrate the correctness and efficacy of PFFT-BEM by existing benchmark solutions by an ideal waveguide shown in Figure 3. Both the top and bottom boundary satisfy pressure release boundary condition with no bottom layer. The sound source is at a depth of 36 m with frequency at 20 Hz. The other water properties of the waveguide are shown in Figure 3. The computational domain size spans \(3\) km \(\times\) 3 km \(\times\) 100 m with uniform meshing on both the top and bottom surfaces using \(\Delta l=\lambda/12\). From the comparison shown in Figure 4a, the transmission loss (TL) results of PFFT-BEM show good agreement with ones of the wavenumber integration method, which validates the accuracy of  
Figure 3: Sketch of the iso-velocity waveguide. PFFT-BEM. The transmission loss (TL) contour obtained by PFFT-BEM in the x-z plane is shown in Figure 4b. In Figure 5, we show the convergence of the numerical results towards the wavenumber integration method results with respect to the boundary element size. We compare the acoustical pressure \(p_{PFFT-BEM}\) by PFFT-BEM with \(p_{wim}\) from wavenumber integration method at z = 36 m. The errors shown in Figure 5 are defined in Equations (24) and (25). As shown in Figure 5, PFFT-BEM results converge clearly towards the wavenumber integration method solution as the element size decreases. \[Average\ error = average\left(\frac{p_{PFFT-BEM}-p_{wim}}{p_{theory}}\right), \tag{24}\] \[Max\ error = max\left(\frac{p_{PFFT-BEM}-p_{wim}}{p_{theory}}\right). \tag{25}\]  
We then check the efficiency of the first-order NRBC at the truncated numerical domain boundary. As shown in Figure 6, PFFT-BEM results clearly improve towards the wave integration method results after radiation condition at the numerical boundary is changed from \(k_{m}=k\) radiation condition  
Figure 4: (a) TL at z = 36 m obtained by PFFT-BEM and wavenumber integration method (chapter 4 of [14]); (b) TL contour obtained by PFFT-BEM. Figure 5: Numerical errors of PFFT-BEM for acoustics pressure at z = 36 m as functions of elements per wavelength. (Equation (16)) to first-order NRBC (Equation (12)) with \(m=2\). In addition, the errors of numerical results using the \(k_{rm}=k\) radiation condition become progressively larger near the truncated domain boundary. When the first order NRBC is used, there are no errors due to the boundary reflection which demonstrate the effectiveness of the first-order NRBCs used by PFFT-BEM for shallow water environments. To demonstrate the computational efficiency of the PFFT-BEM, we vary the boundary element size to increase the number of \(N\). As shown in Figure 7, both the memory and operational count of PFFT-BEM are proportional to O(\(N\)) and much smaller than the ones of the classical BEM, which is proportional to \(O\big{(}N^{2}\big{)}\). This result shows the computational advantage of the PFFT-BEM over the classical BEM for shallow water acoustic simulations. For higher sound source frequency (i.e., smaller acoustics wavelength), the computational cost, which is larger due to the element size, must be much smaller than the acoustics wavelength throughout the entire domain surface. From Figure 8, we show that the PFFT-BEM memory and computational cost is proportional to the sound frequency while it is also much smaller than classical BEM. Figure 6: Transmission Loss (TL) at z = 36m for a 20 Hz point source at a depth of 36 m in the iso-velocity waveguide: comparison of the PFFT-BEM with the \(k_{rm}=k\) radiation condition, PFFT-BEM with first-order NRBC. Figure 7: Comparison of the PFFT-BEM (slope 1) and classical BEM (slope 2) for the iso-velocity waveguide. (**a**) Memory (GB); (**b**) operational counts. In Figure 9, we show the convergence of the numerical results with respect to the number of layers used in the PFFT-BEM. The sound profiles with different gradients, which are used for the convergence tests, are shown in Figure 9a.