## Numerical Method
### Formulation of the Boundary Element Method
\begin{array}{c}\alpha p^{i}-\iint_{S^{i}}p_{n}(\mathbf{\xi})G(\mathbf{x};\mathbf{\xi})d \dot{S}^{i}(\mathbf{\xi})+\iint_{S^{i}}p(\mathbf{\xi})G_{n}(\mathbf{x};\mathbf{\xi})d\dot{S}^{ i}(\mathbf{\xi})+\beta_{i}\iint_{S^{i}+1}p_{n}(\mathbf{\xi})G(\mathbf{x};\mathbf{\xi})d\dot{S}^{ i+1}(\mathbf{\xi})\\ \qquad\qquad\qquad+\iint_{S^{i+1}}p(\mathbf{\xi})G_{n}(\mathbf{x};\mathbf{\xi})d\dot{S}^{ i+1}(\mathbf{\xi})-\iint p_{n}(\mathbf{\xi})G(\mathbf{x};\mathbf{\xi})dS^{Ri}(\mathbf{\xi})\\ \qquad\qquad\qquad+\iint_{S^{Ri}}p(\mathbf{\xi})G_{n}(\mathbf{x};\mathbf{\xi})dS^{Ri}(\mathbf{ \xi})=0,\end{array}  \tag{8}\]  
where \(\alpha\) is the solid angle and \(\beta_{i}=\rho_{i}/\rho_{i+1}\). Using the constant panel method, we discretize the top and bottom surface of the \(ith\) layer into \(N_{i}\) quadrilateral elements and the side boundary, where the radiation boundary condition is applied in layer \(i\), into \(N_{Ri}\) quadrilateral elements. On each element, quantities such as the acoustic pressure \(p\) and the normal gradient of the pressure \(p_{n}\) are assumed to be constant. This transfers Equation (8) into a set of linear algebraic equations of the form  
\[2\pi p^{i}-\sum\limits_{j=1}^{N_{i}}\left(\dot{F}^{i}_{s}\right)_{j}+\sum \limits_{j=1}^{N_{i}}\left  
NRBC is applied at the boundary of the truncated computational domain in PFFT-BEM. In a radial coordinate, the first-order NRBC can be expressed as [13]:  
\[\frac{\partial p_{m}}{\partial r}+\alpha_{m}p_{m}=0,\ m=1,2,\ldots,\ M, \tag{12}\]  
and \(p_{m}\) is the acoustic pressure due to the \(mth\) propagating modes with  
\[\mathrm{p}=\sum_{m=1}^{M}p_{m}, \tag{13}\]  
where \(\alpha_{m}\) is the damping coefficient, which for 3D is defined as:  
\[\alpha_{m}=(\frac{1}{r}-ik_{m}), \tag{14}\]  
where \(k_{m}\) is the horizontal wavenumber for the propagating mode \(m\). Substituting Equation (13) into Equation (12) yields the following boundary condition at the truncated domain boundary:  
\[\frac{\partial p}{\partial r}+\sum_{m=1}^{M}\alpha_{m}p_{m}=0. \tag{15}\]  
In this study, we compute the acoustic pressure for the propagating mode in the far-field \(p_{m}\) and the corresponding wavenumber \(k_{m}(\alpha_{m})\) using the well-known normal mode method. A detailed derivation and good review of the normal mode method can be found in Chapter 5 of [14]. A similar method has been applied for 2D finite element analysis for shallow water acoustics problems [15]. However, for 3D boundary element method, the effect of the first-order NRBC is still unknown and need to be investigated. To demonstrate the efficiency of first-order NRBC, we compare the numerical results using Equation (15) at the truncated boundary with the results based on the radiation condition with \(k_{m}=k\). The latter can be expressed as Equation (16) at the truncated domain boundary:  
\[\frac{\partial p}{\partial r}+\left(\frac{1}{r}-ik\right)p=0. \tag{16}\]  
Substituting Equation (15) into Equation (8) yields  
\[ \begin{split} p^{j}-\iint_{S^{i}}\hat{p}_{n}^{i}(\xi)& G(x\xi)d\hat{S}^{i}(\xi)+\iint_{S^{j}}\hat{p}^{i}(\xi)G_{n}(x\xi)d \hat{S}^{i}(\xi)+\beta_{i}\iint_{S^{i+1}}\hat{p}_{n}^{i+1}(\xi)G(x\xi)d\hat{S} ^{i+1}(\xi)\\ &+\iint_{S^{i+1}}\hat{p}^{i+1}(\xi)G_{n}(x\xi)d\hat{S}^{i+1}(\xi )-\iint_{S^{R}}\sum_{m=1}^{M}\alpha_{m}p_{m}(\xi)G(x\xi)dS^{Ri}(\xi)\\ &+\iint_{S^{Ri}}p(\xi)G_{n}(x\xi)dS^{Ri}(\xi)=0.\end{split}  \tag{17}\]  
If we put Equation (17) in each layer into one single equation, the entire linear system of equations, Equation (17), can then be expressed as  
\[[A][x]=\{b\}, \tag{18}\]  
where \([A]\) is an \(N\times N\) dense and nonsymmetric matrix of influence coefficients, \(\{x\}\) is the vector of unknowns (\(p\) and \(p_{n}\)) on the sub-domain surfaces, and \(\{b\}\) is the known vector due to the point source. An iterative solver such as the generalized minimum residual (GMRES) method could also be applied to solve the system shown in Equation (18). GMRES normally requires \(O(N^{2})\) operations.