A Vector Wavenumber Integration Model of Underwater Acoustic Propagation Based on the Matched Interface and Boundary Method

: Acoustic particle velocities can provide additional energy ﬂow information of the sound ﬁeld; thus, the vector acoustic model is attracting increasing attention. In the current study, a vector wavenumber integration (VWI) model was established to provide benchmark solutions of ocean acoustic propagation. The depth-separated wave equation was solved using ﬁnite difference (FD) methods with second- and fourth-order accuracy, and the sound source singularity in this equation was treated using the matched interface and boundary method. Moreover, the particle velocity was calculated using the wavenumber integration method, consistent with the calculation of the sound pressure. Furthermore, the VWI model was veriﬁed using acoustic test cases of the free acoustic ﬁeld, the ideal ﬂuid waveguide, the Bucker waveguide, and the Munk waveguide by comparing the solutions of the VWI model, the analytical formula, and the image method. In the free acoustic ﬁeld case, the errors of the second- and fourth-order FD schemes for solving the depth-separated equation were calculated, and the actual orders of accuracy of the FD schemes were tested. Moreover, the time-averaged sound intensity (TASI) was calculated using the pressure and particle velocity, and the TASI streamlines were traced to visualize the time-independent energy ﬂow in the acoustic ﬁeld and better understand the distribution of the acoustic transmission loss.


Introduction
Ocean acoustic waves propagate over long distances because their energy attenuates less in water than does that of electromagnetic waves; consequently, acoustic waves have been widely used in target detection, environmental monitoring, and underwater communication [1]. Underwater acoustic wave propagation excited by a time harmonic source can be described by the Helmholtz equation in the frequency domain, which can be solved using the wavenumber integration technique [2,3], the fast field program (FFP) [4][5][6], normal modes [7], ray theory [8], the parabolic equation [9], the finite difference method (FDM) [10], and the spectral method [11]. As these models for predicting the scalar sound level have reached a high level of development, our interest herein is in developing a vector model to accurately predict the particle velocity.
As a result of the growing interest in observing and investigating the underwater vector acoustic field [12], research is increasingly being conducted on the structure of the local vortices of the energy flux vector in acoustic waveguides [13]. Compared with experiments, vector acoustic simulations can more rapidly and economically provide more acoustic fields; accordingly, computer simulation and visualization techniques are highly important for studying the vector acoustic field [14][15][16][17]. Although the particle velocity can be obtained by calculating the pressure gradient using the FDM [18], the grid sizes in all coordinate directions must be sufficiently small to reduce numerical errors, resulting in a large number of grid points. Hence, other methods that can be used to analytically calculate some of the particle velocity components need to be developed, such as the analytical Lie-algebraic approach [19,20] and acoustic vector models that compute the particle velocity directly. Smith [21] introduced general numerical algorithms of the acoustic particle velocity field from the parabolic equation and normal mode models, and test cases of a Pekeris waveguide and range-dependent benchmark wedge were presented to verify the results of both models. Gulin [22] developed an analytical numerical acoustic model in terms of a scalar-vector description; the Helmholtz equation was replaced by two first-order partial differential equations that were used to govern the sound pressure and the velocity vector of fluid particles, and the normal mode method was utilized to simulate the acoustic scalar-vector field in layered stratified media. The FFP, which has become the most common wavenumber integration method in underwater acoustics, uses the far-field approximation and fast Fourier transform (FFT) to accelerate the computation. For instance, Zhu [23] developed a vector acoustic FFP model for an ocean-like environment with a sound speed profile of water and a multilayered elastic sediment bottom and analyzed the correlation between the sound pressure and the particle velocity.
The objective of this study was to obtain the most accurate vector acoustic model possible that can be used to provide benchmark solutions in range-independent environments for verifying other numerical acoustic models, such as FDM models [10] and spectral method models [11]. However, although the normal mode method, parabolic equation, and FFP mentioned above have the advantages of a small computational load and good numerical stability, they have difficulty effectively satisfying the required accuracy due to approximation errors. For instance, parabolic and far-field approximations should be used in the parabolic equation [24]; the branch cut integral is commonly neglected in normal modes [25]; and the far-field approximation needs to be applied to the FFP [26,27]. In contrast, the wavenumber integration method is attractive because it is basically an implementation of the integral transform technique for general range-independent media and directly evaluates the integrals of the solutions to the depth-separated wave equation (depth equation for short) by numerical quadrature; consequently, the field solution can be accurately calculated to any desired degree [2] (page 233). However, the wavenumber integration method incurs a relatively large computational load; nevertheless, for the model to generate benchmark solutions, guaranteeing the accuracy of the model is most often the greatest challenge, with reducing the computation time usually being a secondary priority [28].
To address the abovementioned challenges, a vector wavenumber integration (VWI) model was established in this study and verified using benchmark cases of underwater acoustic propagation. The primary contributions of this study are as follows: (1) An FDM to solve the depth equation is proposed. Several methods exist for solving the depth equation, such as the propagator matrix approach [4] and the direct global matrix approach [29]; however, these approaches require that the density be approximated by the piecewise constant in each element along the depth direction, meaning that the rangeindependent environments are approximated as layered media. To further improve the solving accuracy of the depth equation solution, FD schemes with second-order (2nd-order) and fourth-order (4th-order) accuracy that can treat the density variation are proposed in this paper, and the recently developed matched interface and boundary (MIB) method is used to address the sound source singularity [30,31].
(2) A VWI model is established. In the present model, the particle velocity at any position in the sound field can be directly calculated using wavenumber integration (consistent with the calculation of the sound pressure sharing the zero-order Bessel function of the first kind) instead of calculating the spatial gradient of the pressure (two-dimensional grids are necessary, and the grid sizes in both directions need to be sufficiently small [18]).
(3) A sound field analysis method that displays both propagation loss contours and streamlines of the time-averaged sound intensity (TASI) is proposed. The TASI can be computed using the pressure and particle velocity and displayed using streamlines, thereby bettering our understanding of the transmission loss (TL) distribution and the net transport of sound energy. For instance, TASI streamlines in the test cases show that the regions through which the energy flow is concentrated form zones of sound convergence.

The FDM for Solving the Depth Equation
Within a range-independent marine environment, because a simple point source is omnidirectional, the solutions of the steady-state acoustic field are symmetric about the z-axis and independent of the circumferential direction in a cylindrical coordinate system; thus, the Helmholtz equation in the frequency domain can be expressed as [2]: where z is the vertical coordinate (in m), which is positive in the direction of increasing depth and passes through the source point; z s is the source depth (z-coordinate of the source); r is the coordinate (in m) in the range direction from the source; P(r,z) is the relatively complex sound pressure field (referenced to the pressure at a distance of 1 m from the sound source) and is assumed to have a time factor of e −iωt ; ρ(z) is the medium density (in g/cm 3 ); δ is the Dirac delta function; and k(z) = 2πf /c(z) is the wavenumber (f is the sound frequency in Hz and c(z) is the speed of sound in the medium in m/s). We can apply a Hankel transform in the range direction to Equation (1) and define where φ is the wavenumber kernel function (kernel for short), k r is the horizontal wavenumber, and J 0 is the zero-order Bessel function of the first kind. Then, the sound pressure can be expressed in terms of the Sommerfeld integral: Multiplying Equation (1) by J 0 (k r r)rdr and integrating over r from 0 to ∞ ( ∞ 0 (·)J 0 (k r r)rdr), we can obtain the depth equation: where k 2 z = k 2 − k 2 r . The wavenumber integration model consists of two successive steps: one is to obtain the kernel by solving Equation (4); the other is to attain the pressure by numerically integrating Equation (3).

The Boundary and Source Interface Conditions
As shown in Figure 1, the kernel (φ) is continuous at the upper boundary, and its derivative with respect to the depth (φ ) satisfies where the subscripts "0 − " and "0 + " represent the upper and lower surfaces, respectively, of the upper boundary (grid point 0). The derivative inside the water (φ 0 + ) can be calculated using kernels on the interior grid points by the FDM. To calculate the derivative outside the water (φ 0 − ), considering the homogeneous properties of the half-space above the water column, Equation (4) can be changed to the plane wave equation away from the source: Because only upgoing waves are present in the upper half-space, the solution of Equation (6) is: where C is a constant. Then, φ 0 − can be expressed by the kernel on the upper boundary: The lower boundary conditions at grid point N can be treated as symmetric (only downgoing waves are present) to the upper boundary conditions.
The source interface conditions state that the kernel must be continuous (φ z − s = φ z + s ), and the derivative of the kernel follows the jump condition. Let z − s = z s − ε and z + s = z s + ε, where ε → 0 . Then, Equation (4) is integrated over depth z from z − s to z + s along the z-axis through the source. We can obtain the jump condition of the derivative of the kernel:

The FDM for the Depth Equaiton
As shown in Figure 1, equally-spaced grids in the depth direction are used in this paper, and h is the grid step size in meters. At the interior grid points far from the boundaries and source, the second derivatives in Equation (4) can be discretized using the standard 2nd-order accuracy FD scheme: where ρ i+1/2 = ρ(z i+1/2 ) and ρ i−1/2 = ρ(z i−1/2 ). Then, the discretized depth equation at interior points can be expressed as: At the lower surface of the upper boundary, the derivative of the kernel (φ ) can be approximated using the 2nd-order FD scheme: After inserting Equations (5) and (8), we can establish the equation: and we can attain the following relationship between the kernel on the upper boundary φ 0 and the kernels at the interior points (φ 1 and φ 2 ): Specifically, for the pressure-release upper boundary, we can set the density in the upper half-space to zero (ρ 0 − = 0, and thus, φ 0 = 0). The lower boundary can be treated as symmetrical with respect to the upper boundary, and thus, we can derive the following relationship between the kernel on the lower boundary (φ N ) and the kernels at the interior points (φ N−1 and φ N−2 ): At grid points near the source interface, the MIB method is used to form the FD schemes. The MIB is a numerical method that can treat discontinuous interfaces: the solution on each side of the source interface is smoothly extended beyond the interface by means of fictitious domains, and the fictitious values in the fictitious domains are determined simultaneously by enforcing the source interface conditions at the exact position of the interface [31,32]. Assuming that the source interface lies between the grid nodes s and s + 1, the interface splits this grid step into two pieces with sizes of θh at the upper part and (1 − θ)h at the lower part, where 0 ≤ θ < 1. Because the derivative is discontinuous at the source interface, the MIB method smoothly extends the solution on the upper side of the source interface beyond the interface by means of fictitious domains and sets the fictitious kernel F s+1 at the location of grid point s + 1; symmetrically, the fictitious kernel F s is supplied at the location of grid point s to smoothly extend the solution from the lower side of the source interface to the upper side. Using Lagrange interpolation polynomials, the kernel at the source interface can be interpolated from the upper and lower sides: Moreover, the derivatives of the kernel at the source interface can be approximated using the 2nd-order FD scheme from the upper and lower sides: Applying source interface conditions of φ z − s = φ z + s and Equation (9), the fictitious kernels in the fictitious domains can be determined simultaneously: Thus, by using the above FDM with 2nd-order accuracy, a linear system of the discretized depth equation at all grid points (excluding the boundary points 0 and N) can be obtained and solved using direct methods for solving linear systems: and To improve the numerical accuracy of the kernels, a 4th-order FD scheme is also developed, which is described in Appendix A.

VWI Approaches
In addition to the sound pressure P, the VWI model also needs to calculate the acoustic particle velocity. According to Euler's equation, the acoustic particle velocity can be calculated from the spatial gradient of pressure [33]: In the present VWI model, the particle velocity at any position in the acoustic field is directly calculated using the wavenumber integration method without a two-dimensional grid. Using Equations (3) and (29), the horizontal and vertical components of the particle velocity can be expressed as follows: where P r = ∂P/∂r, P z = ∂P/∂z, J 1 is the first-order Bessel function of the first kind (∂J 0 (k r r)/∂r = −k r J 1 (k r r)), and w = φ /(iρω) is the kernel of the vertical component of the particle velocity and can also be calculated using the FDM. For the 2nd-order FD scheme: For the 4th-order FD scheme: In the cycles of wavenumber integration, the zero-order Bessel function of the first kind J 0 (k r r) can be shared by the integrals of the pressure P in Equation (3) and the vertical component of the particle velocity W in Equation (31).

Algorithm Parameters
To numerically evaluate the wavenumber integrals in Equations (3), (30) and (31), the kernel needs to be evaluated at discrete horizontal wavenumbers. In this paper, k r ∈ [0, k max ] is discretized equidistantly to N k cells, and the kernel at the midpoint of each cell is used as the average value of that cell for wavenumber integration. Thus, discretized horizontal wavenumbers used to solve the depth equation can be expressed as: where ∆k r is the horizontal wavenumber step size. To achieve more accurate integrals, ∆k r should be small enough to avoid aliasing and wrap-around effects [2], which means that a sufficient number of sampling points is needed at the maximum range R max . In this paper, the horizontal wavenumber step size is estimated using: where n R is the number of samples of the horizontal wavenumber in each 2π period (the Bessel functions exhibit periodic oscillations similar to the cosine function) at r = R max . Then, the total number of samples of the horizontal wavenumber can be calculated by N k = k max /∆k r , where k max is the finite maximum value of the horizontal wavenumber used to truncate the range of the infinite Sommerfeld integrals.
In addition, to avoid the vertical wavenumber k z = 0, the integration path is commonly moved into the complex plane [34], which means that: where ε k = 3∆k r / 2π log 10 e and e is the natural constant. Then, the discrete wavenumber integrals of Equations (3) and (29) can be expressed as:

Vector Sound Intensity
After the particle velocity field is calculated using the VWI method, the TASI streamlines can be traced to visualize the time-independent energy flow in the acoustic field. Although sound intensity streamlines have been applied in several studies [35][36][37], these streamlines are rarely shown together with the TL contours to better understand the TL distribution.
For a steady-state (time-harmonic) sound field, the TASI I can be obtained by temporally averaging the instantaneous intensity over one time period (T = 2π/ω), which can be expressed as [35]: Additionally, the two components of I can be expressed as: The TASI governs the net transport of energy, whose direction is proportional to the spatial gradient of the pressure phase and perpendicular to the resultant wave fronts (surfaces of the constant-pressure phase). The TASI streamline curves are tangent to the direction of the local TASI vector at every point [35,38], and these streamlines can describe the way in which sound energy is transferred, physically coming only from the sound source and propagating along the streamlines to infinity [39].

Benchmarking for Vector Acoustic Fields
Four benchmark cases were used to test the present VWI model: the free acoustic field, the ideal fluid waveguide, the Bucker waveguide, and the Munk waveguide. Moreover, the TASI vectors and streamlines are calculated and traced to visualize the time-averaged energy flow path, and the TL distribution is interpreted from the perspective of streamlines.

The Free Acoustic Field
If the acoustic environment is a homogeneous and unbounded medium (free space), the solution of the acoustic Helmholtz equation is analytical and can be expressed using the Sommerfeld-Weyl integral [2] (page 88): where R = r 2 + (z − z s ) 2 is the distance from the source point (0, z s ) to the position (r, z). Moreover, the exact kernel is: which can be used to verify the numerical solution of the kernel solved by using the FDM. Herein, the thickness of the uniform medium (water) is 5000 m, the sound speed and density of the medium are 1500 m/s and 1 g/cm 3 , respectively, no attenuation occurs within the medium, and the properties of the medium in the upper and lower half-space are the same as those of the water. The depths of the sound source and the receiver are First, the FDM is verified using Equation (42). Setting the grid step size in the depth direction to h = 20 m, comparisons of the kernel (k r = k r,1 ) obtained using the exact formula and the 2nd-and 4th-order FD schemes are shown in Figure 2a,b, demonstrating that the three results agree well and that the solutions of the 4th-order FD scheme are more accurate than those of the 2nd-order FD scheme. Furthermore, the actual orders of the formal 2nd-and 4th-order FD schemes were tested. The relative error of the kernel at the ith grid point can be expressed as: where φ FDM i is the numerical solution of the kernel on grid point i and φ Exact i is the exact solution. In this paper, two grid size steps (in the depth direction) of h 1 = 20 m and h 2 = 10 m are used to solve the depth equation. Then, at the depth of z = ih 1 = 2ih 2 , the actual orders of the FD schemes can be calculated using: Figure 3a,b show the relative errors of the kernel and actual orders with respect to the formal 2nd-and 4th-order FD schemes (k r = k r,1 ), respectively, demonstrating that the relative errors can be decreased by reducing the grid step size; the relative errors near the source depth (2500 m) are smaller than those near the boundaries; the 4th-order FD scheme has smaller relative errors than the 2nd-order FD scheme; and the actual orders of the 2ndand 4th-order FD schemes are close to their formal orders. It is worth pointing out that the MIB method exhibits good performance at the discontinuous source interface, meaning that the actual orders of the FD schemes do not decrease near the source interface.  Figure 4 shows the TL curves along the receiver line z = 20 m, where TL = −20 log 10 |P| (in dB) and |P| denotes the magnitude of the complex pressure P. The solutions of the 2nd-and 4th-order FD schemes (h = 20 m) match the exact analytical formula well in general, and the curve of the 4th-order FD scheme is closer to that of the exact formula than that of the 2nd-order FD scheme.  Figure 5 plots the TL contours and TASI streamlines obtained from the VWI model using the 2nd-order FD scheme (a) and the 4th-order FD scheme (b). In both figures, the TL contour lines are approximately concentric circles (the 4th-order scheme is more accurate than the 2nd-order scheme), whereas the TASI streamlines are straight lines from the source point to infinity, which is consistent with the characteristics of spherical waves in free space.

The Ideal Fluid Waveguide
The ideal fluid waveguide is represented by a range-independent homogeneous water column with pressure-release upper and lower boundaries [2] (page 102). The density and sound speed of the water are ρ w = 1 g/cm 3 and c w = 1500 m/s, respectively, and the densities of the medium in the upper and lower half-spaces are both 0 g/cm 3 . The parameters of the sound field are as follows: the depth of the seabed is 100 m, the depths of the sound source and receiver line are 35 m and 25 m, respectively, the largest receiver range is 3 km, and the source frequency is 20 Hz. In the VWI model, n R = 10 and k max = 10 k ref . The interval step sizes in the range and depth directions are both 5 m.
The image method can be used to calculate the theoretical solutions of the ideal fluid waveguide, and the acoustic pressure can be expressed using the summation of the physical source and an infinite number of image sources as follows (in this paper, the summations of the image method are truncated at m = 500): with and H is the depth of the seabed in meters. Figure 6 shows the TL curves along the receiver line calculated by the image method and the present VWI model using the 2nd-and 4th-order FD schemes. The three curves generally match well. Moreover, at the peaks and troughs of the curves, the solutions of the 4th-order FDM are closer to those of the image method than those of the 2nd-order FDM. The TL contours and TASI streamlines obtained from the present VWI model using the 2nd-and 4th-order FD schemes are shown in Figure 7a,b, respectively, and they match well throughout the whole field. The main characteristic of the sound field is that the TASI streamlines are periodic in the range direction.

The Bucker Waveguide
The Bucker waveguide consists of a shallow water waveguide containing a bilinear sound speed (in m/s) profile in the water column [40]: The depth of the water is 240 m, and the density of the water is ρ w = 1.0 g/cm 3 . The density of the medium in the upper half-space is 0 g/cm 3 , whereas the homogeneous lower half-space of the penetrable sediment layer has a sound speed of 1505 m/s and a density of 2.1 g/cm 3 , and no attenuation occurs in this problem. The depths of the sound source and the receiver line are 30 m and 90 m, respectively, the largest receiver range is 2 km, and the source frequency is 100 Hz. In the VWI model, n R = 10 and k max = 10 k ref . The interval step sizes in the range and depth directions are both 1 m. Because the density contrast in this waveguide yields a significant number of virtual modes close to the real wavenumber axis, normal mode models ignoring the continuous spectrum cannot provide accurate solutions. Alternatively, wavenumber integration models have no restrictions on the density contrast and are therefore capable of providing exact solutions for this environment [2] (page 300). Figure 8 shows the TL curves along the receiver line calculated by the present VWI model using the 2nd-and 4th-order FD schemes, demonstrating that the two curves match well in all ranges. Figure 9 shows the TL contours obtained from the present VWI model using the 4th-order FDM, revealing that the acoustic regions with denser TASI streamlines have a smaller TL, which means that the energy distribution in the sound field is not uniform.

The Munk Waveguide
The Munk waveguide has the following analytical sound speed (in m/s) profile in water [2] (page 356): where the scaled depth z = 2(z − 1300)/1300. The depth of the water is 5000 m, and the density of the water is ρ w = 1.0 g/cm 3 . The density of the medium in the upper half-space is 0 g/cm 3 , whereas the homogeneous lower half-space of the penetrable sediment layer has a sound speed of 1600 m/s, a density of 1.8 g/cm 3 , and an attenuation rate of 0.8 dB/λ. The depths of the sound source and the receiver line are both 100 m, and the source frequency is 50 Hz. In the VWI model, n R = 8 and k max = k ref . The interval step sizes in the range and depth directions are 300 m and 3.75 m, respectively. Figure 10 shows the TL contours obtained from the present VWI model using the 2nd-order FDM and the acoustic rays calculated using the BELLHOP model in the acoustic toolbox [41]. Figure 11 shows the TL contours and TASI streamlines obtained from the VWI model using the 4th-order FDM, indicating that the acoustic regions with denser TASI streamlines have a smaller TL, and the regions through which the energy flow is concentrated form zones of sound convergence. Moreover, the TL contours in Figures 10 and 11 match well, but the acoustic rays are significantly different from the TASI streamlines due to the following reasons. First, the TASI streamlines are calculated using the sound field solutions of the vector acoustic model, but the geometric rays are traced to calculate the solutions of the acoustic field. Second, because only one TASI vector direction may exist at any point in the field, the TASI streamlines cannot intersect; in contrast, the geometric rays can intersect among the multipath arrivals of propagating wavefronts in an ocean-like waveguide. Third, the TASI streamlines can accurately visualize the time-averaged energy flux, whereas the geometric rays intuitively visualize the approximate propagation of an acoustic field from the source to a receiver and do not include the diffracted components of the field.  Figure 12 shows the TL curves along the receiver line calculated by the present VWI model using the 2nd-and 4th-order FD schemes, and the two curves match well in all ranges. To test the effect of the density variation with depth on the acoustic field, a density profile in water is used for the Munk waveguide (other properties remain unchanged): where H = 5000 m is the seabed depth. Figure 13 shows the TL curves along the receiver line for the Munk waveguide with different density profiles in water calculated by the present VWI model using the 4th-order FDM, demonstrating that the density variation of water has an important influence on the peaks and troughs of the TL curves.

Conclusions and Discussion
The present range-independent VWI model can serve as a benchmark vector model for other acoustic models, such as FDM models and spectral method models. Moreover, by applying the particle velocity, the TASI streamlines can be computed and displayed, thereby bettering our understanding of the TL distribution, and these TASI streamlines can also be applied to other vector acoustic models. Furthermore, the particle velocity can be used to investigate the vortex-like structure and Stokes parameters [42] of the vector acoustic field.
It should be mentioned that the present VWI model is computationally intensive because it includes three Sommerfeld integrals (the pressure and two velocity components), and the linear system of the discretized depth equation using the FDM is solved directly because it is difficult to use the efficient forward elimination and backward substitution algorithm since matrix A is not a diagonal banded form as a result of the MIB method. Thus, parallel computing techniques [43] should be applied to the VWI model in future research to accelerate the computation.
In addition to the FDM combined with the MIB method in this paper, finite element methods and spectral methods [7] with high orders of accuracy can also be used to solve the depth equation, and a performance comparison among these methods may be studied in the future. Furthermore, the MIB method shows good performance in addressing the discontinuous source interface of the one-dimensional depth equation, and thus, it is very attractive to apply the MIB method to treat seabed discontinuities in FDM models [44].  Acknowledgments: The authors are grateful to M. Porter for providing the BELLHOP program [41].

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. The FD Scheme with 4th-Order Accuracy
The derivation of the FD scheme with 4th-order accuracy is similar to that of the FD scheme with 2nd-order accuracy, but the computing stencils are wider. The computational stencil coefficients of the derivatives can be calculated using Taylor series and the method of undetermined coefficients [45] (or Lagrange interpolation polynomials).

Appendix A.1. Interior Grid Points
At the interior grid points away from the boundaries and source interface, the second derivatives in Equation (4) can be discretized using the cell-centered FD scheme with 4th-order accuracy [44,46]: and the first derivatives at the half-integer points can be discretized using the same FD scheme: Then, the discretized depth equation at the interior points can be expressed as:

. Grid Points near the Boundaries
At the lower surface of the upper boundary, the derivative of φ can be approximated using the 4th-order FD scheme: After inserting Equations (5) and (8), we can obtain the following relationship: The 4th-order FD scheme of the second derivative at grid point 1 is: where [47] The stencil coefficients of the FD scheme at the lower boundary point are antisymmetric to those at the upper boundary point.

Appendix A.3. Grid Points near the Source Interface
The kernels at the source interface can be interpolated from the upper and lower sides: The derivatives of the kernel at the source interface can be approximated using the 4th-order FD scheme from the upper and lower sides: Applying the source interface conditions of Equation (9) and φ z − s = φ z + s , we can obtain: where γ = w + 1,s w − 0,s+1 − w − 1,s+1 w + 0,s . The 4th-order FD schemes of the second derivative at grid points s and s + 1 are: The other first derivatives can be approximated using Equation (A2).