## 1 Introduction  
In the shallow water, the vibrations and radiated noises of sailing ships propagate to the seafloor through seawater, stimulating seismic waves. Since undersea seismic waves carry useful target information as a complement to sound waves in water, they have attracted much attention in recent years in underwater target check detection [1, 2], target identification [3], and acoustic fuzing of mines [4]. However, submarine seismic wave tests require a lot of human and material resources, so the numerical prediction of the shallow submarine seismic wave field has become an effective way to study its characteristics [5, 6]. The high-precision and high-efficiency submarine seismic wave field prediction is of great significance for the application of research on submarine seismic waves. Acquiring shallow seafloor seismic wavefields requires the time-domain solution of the wavefield equation. Currently, scholars mainly focus on the time-domain solution methods for the single fluid or elastomer fluctuation equation. However, the submarine seismic wavefield equation is obtained by coupling the fluid acoustic pressure fluctuation equation and the elastomer displacement equation, which introduces a coupling term compared with the single equation. Therefore, the single equation solution method is not fully applicable to the solution of the coupled equations. For a single fluid fluctuation equation or elastomer fluctuation equation, both system equations can be expressed by the system mass matrix and stiffness matrix without considering damping, and the only difference lies in the field variables to be solved, which are generally the sound pressure for the fluid fluctuation equation and the displacement for the elastomer fluctuation equation. The classical methods for solving the system equations include the explicit center difference method and the implicit Newmark method [7], based on which scholars have conducted a series of research studies. For the fluid problem, [PERSON] [8] and others compared the performance of the explicit central difference method and the implicit Newmark method for the two-dimensional fluctuation equations with a simple example. The results show that the stability of the implicit Newmark method is better than that of the explicit method. However, the accuracy of the explicit method is higher than that of the Newmark method within the stabilization step of the explicit method, and the analysis is not convincing and representative in the simulation of ocean acoustic waves because the example is too simple. [PERSON] et al. [9] combined the explicit Newmark method with the isogeometric discretization method for acoustic wave problems with absorbing boundary conditions and analyzed the stability of the method. Although the accuracy and efficiency of the method are better than the finite difference method and the finite element method in some structural acoustics and wave propagation studies, the isogeometric discretization method is severely limited by the size of the discretized mesh and the interpolating polynomial order, which need to be simulated several times to choose the appropriate values. In addition, how to configure the parameters of the isogeometric discretization method is yet to be researched for the elastic wave problems similar to the shallow water seismic waves. For the elastomer problem, [PERSON] [10] and others proposed a weighted velocity Newmark method to suppress the pseudo-oscillation phenomenon in the dynamic structural rupture problem. However, the selection of the weights is not sufficiently theoretical supported. [PERSON] [11] and others proposed a bisecting step generalized central difference method for general dynamics problems, but it is complicated to solve the relevant parameters for multiple degrees of freedom. The time step is divided into two steps, greatly increasing the computational volume.