## Numerical Method  
### Formulation of the Boundary Element Method  
Consider the following 3D Helmholtz wave equation:  
\[\rho\
abla\left(\frac{1}{\rho}\
abla p(\mathbf{x})\right)-k^{2}p(\mathbf{x})=0,\ \ \mathbf{x}\in\mathbf{V}, \tag{1}\]  
where \(\mathbf{x}=(x\),\(y\),\(z\)), \(p(x)\) is acoustics pressure, and \(\rho\), \(\mathrm{V}\), \(\mathrm{k}\) is the density, volume and wavenumber of the medium, respectively. At the interface between the medium and atmosphere, pressure release boundary condition is satisfied as:  
\[p(\mathbf{x})=0,\mathbf{x}\in S_{F}, \tag{2}\]  
where \(S_{F}\) is the free surface. The Sommerfeld radiation condition is defined at \(S_{\infty}\):  
\[\frac{\partial p(\mathbf{x})}{\partial n}-ikp(\mathbf{x})=0,\ \ \mathbf{x}\to\infty, \tag{3}\]  
where \(n\) is the normal direction on the boundaries. For shallow water environments, a multilayer model is adopted to account for the inhomogeneous medium properties. In this model, we assume the density \(\rho\) and sound speed \(c\) in each layer to be constant at every instant of time. A sketch of the multilayer domain is shown in Figure 1.  
Then, at every instant of time, the homogeneous Helmholtz equation governing the sound field in each layer can be expressed as  
\[\
abla^{2}p^{i}(\mathbf{x})+k^{2}p^{i}(\mathbf{x})=0,\ \ \mathbf{x}\in v^{i}, \tag{4}\]  
where \(i\) is the number of layers with \(i=1\), \(\Delta\), \(I\), and \(v^{i}\) is the volume of layer \(i\). Continuity of pressure and particle velocity are required across the interface \(\hat{S}^{i+1}(\hat{S}^{i})\) of neighboring layers \(v^{i+1}\) and \(v^{i}\) (including at the penetrable bottoms). We define the bottom and top boundaries of the layer \(i\) as the symbols '\(\mathcal{I}\)' and \(\mathcal{V}\)', respectively. Continuity boundary condition across the interface can be expressed as  
\[p^{i}=p^{i+1},\ \ \mathbf{x}\in\hat{S}^{i+1}, \tag{5}\]  
\[\frac{1}{\rho_{i}}\frac{\partial p^{i}}{\partial n}=\frac{1}{\rho_{i+1}}\frac {\partial p^{i+1}}{\partial n},\ \ \mathbf{x}\in\hat{S}^{i+1}. \tag{6}\]  
Figure 1: Schematic of the multilayer model with a total number of layers of \(\mathrm{I}=3\). From [PERSON]'s second theorem, the boundary integral equation (BIE) for the pressure \(p^{i}(x)\) in each layer is derived from Equation (1) as  
\[4\pi p^{i}-\iint\limits_{S^{i}+S^{ki}}p^{i}_{n}(\xi)G(\mathbf{x};\mathbf{\xi})dS^{i}(\bm {\xi})+\iint\limits_{S^{i}+S^{ki}}p^{i}(\xi)G_{n}(\mathbf{x};\mathbf{\xi})dS^{i}(\xi)=0, \tag{7}\]  
where \(S^{i}\) represents the top and bottom surface of each layer, \(S^{Ri}\) is the radiation boundary in the far field, \(G(\mathbf{x};\mathbf{\xi})=\frac{1}{|\mathbf{x}-\xi|}e^{ik|\mathbf{x}-\mathbf{\xi}|}\) and \(\mathbf{\xi}\) is a dummy variable representing the point on the layer boundary \(S\). Because the Green function \(G(\mathbf{x};\mathbf{\xi})\) satisfies the radiation condition Equation (3) itself, the bottom boundary \(S^{\infty}\) does not need to be meshed in order to satisfy the far-field radiation condition. Similarly, due to the attenuation effect in the bottom layer, the radiation condition Equation (3) is also automatically satisfied at the side boundary \(S^{RI}\) in the bottom layer. However, due to the finite truncation of the side domain boundary, a proper radiation condition needs to be applied on \(S^{Ri}\) in the fluid layer (i.e., \(i\
eq 1\)). This is further discussed in Section 2.2. By taking the limit of \(\mathbf{x}\to S^{i}\) and using the far-field boundary conditions, pressure release boundary condition, and continuity interface boundary conditions, Equation (7) for each layer can be coupled into one large set of integral equations  
\[