<p>Wavenumber integration (WI) theory is typically employed in numerical simulations of underwater acoustic propagation, for which the solution of depth-dependent kernel functions is the core issue of integral transform modeling. In the current widely used WI models, the depth equation is discretized using the finite element method or a spectral method. The finite element model is characterized by slow error convergence but high computational efficiency, whereas the spectral model is characterized by low computational efficiency but high accuracy. We introduce a spectral element approach for evaluating kernel functions by discretizing the depth-separated wave equation in weak form. The domain is partitioned into elements, and within each, Lagrange interpolants constructed at Gauss–Legendre–Lobatto nodes serve as basis functions. This yields a global matrix that is both block-diagonal and symmetric. The results of the numerical scheme and simulation experiments show that the newly developed spectral element algorithm and model not only combine the advantages of both methods but can also flexibly highlight the strengths of either method via optimal parameter setting, demonstrating high application value.</p>