## 2 Forward Modeling of Pavement Structures
### Spectral Element Method
The potential function \(\alpha\) in Equation (6) in the frequency-wavenumber domain can be obtained by using the variable separation method:  
\[\hat{\alpha}(r,z)=A_{n}e^{-ik_{Pz}z}J_{0}(kr), \tag{8}\]  
where superscript \"\({}^{...}\)\" represents the function of the frequency-wavenumber domain, \(A_{n}\) represents an unknown constant which is determined by the boundary conditions, \(J_{0}\) is the first kind of Bessel function, and \(k\) is the wavenumber. The expression of \(k_{Pz}\) is  
\[k_{Pz}=\left(\frac{\omega_{n}^{2}}{c_{P}^{2}}-k^{2}\right)^{\frac{1}{2}}, \tag{9}\]  
where \(\omega_{n}\) is the frequency. To introduce the boundary conditions in the radial direction, the amplitude of the oscillation at the radial boundary \(r=R\) (far away from the load source) is considered to vanish. That requires the item \(J_{0}(kR)\) to be equal to 0. This condition can be satisfied at the infinite positive roots \(k_{m}R\) of the \(J_{0}\) function. Hence, the m functions,  
\[\hat{\alpha}_{mn}(r,z)=A_{mn}e^{-ik_{Pz}z}J_{0}(k_{m}r), \tag{10}\]  
In the same way, the solution of the potential function \(\beta\) in the frequency-wavenumber domain in Equation (7) is  
\[\hat{\beta}_{mn}(r,z)=B_{mn}e^{-ik_{Pz}z}J_{1}(k_{m}r), \tag{11}\]  
where \(B_{mn}\) is a constant determined by the boundary conditions and \(J_{1}\) is the first kind of Bessel function. Moreover, \(k_{Sz}\) represents the vertical shear wavenumber, and is expressed as  
\[k_{sz}=\left(\frac{\omega_{n}^{2}}{c_{SH}^{2}}-{k_{m}}^{2}\right)^{\frac{1}{2}} \tag{12}\]  
Because Equations (6) and (7) are linear homogeneous formulas, the superposition principle is applicable, i.e., the solutions of the potential functions \(\alpha\) and \(\beta\) can be obtained by the superposition of the infinitely many solutions of \(\hat{\alpha}_{mn}\) and \(\hat{\beta}_{mn}\). The inversion of the integral transformation is usually solved using infinite numerical integration [6; 26]. A high number of wavenumbers must be calculated to approximate the infinite integration. In order to obtain an accurate result, the step of the integration needs to be divided very small, which can reduce the computational efficiency significantly. In SEM, the double-sum of the discrete angular frequency \(\omega_{n}\) (\(n=1,\ 2,\ \ldots,N\)) and wavenumber \(k_{m}\) (\(m=1,2,\ldots,M\)) is used to obtain all arrays of the whole system dissipated at \(r=R\) (\(R\) is a large finite value). The summation over the wave numbers in SEM is discontinuous, and the wave number is finite. Therefore, the general solution of the wave formula of an axisymmetric system can be obtained by superimposing the special solution [10], as follows. \[\alpha(r,z,t)=\sum_{N}\sum_{M}A_{mm}e^{-ik_{Pz}z}J_{0}(k_{m}r)e^{-i\omega_{n}t} \tag{13}\]  
\[(\lambda+2\mu)\bigg{(}\frac{\partial^{2}\alpha}{\partial r^{2}}+\frac{ \partial^{2}\alpha}{\partial z^{2}}+\frac{\partial\alpha}{r\partial r}\bigg{)} =\rho\frac{\partial^{2}\alpha}{\partial t^{2}} \tag{6}\]0  
The finite values \(N\) and \(M\) in the general solution Equations (13) and (14) can be determined by the time and spatial amplitude spectra of the excitation load. The sum of \(N\) frequencies can be obtained by using the fast Fourier transform (FFT), and the sum of wavenumbers can be obtained by using the Fourier-Bessel series. The sum of wavenumbers replaces the integral from 0 to infinity in the conventional numerical integration method, which is crucial to improving the efficiency of the algorithm. For the infinite integral, when the damping ratio of the medium is zero, the integrand has singular points, which makes its integral calculation more expensive. The sum of the wavenumbers in double-sum form can effectively avoid this problem. The general solutions of the potential functions \(\alpha\) and \(\beta\) (Equations (13) and (14)) are substituted into Equation (3) to obtain the displacements of the axisymmetric problem, as follows:  
\[ \begin{cases}u_{r}(r,z,t)\\ u_{z}(r,z,t)\end{cases} =\sum_{n}\sum_{m}