## 1 Introduction  
Fluid-saturated porous media is a widespread saturated soil, which is a combination of soil skeleton and fluid filled in the pores. The porous media theory and computational schemes have an important influence in geomechanics. Specifically, porous media analyses govern the physical behavior of cohesive soils consolidation and failure [1; 2; 3; 4] and cohesionless soils [5; 6] response when subjected to static and dynamic loading. The mechanical properties of the pore fluid and the fluid-solid couple interaction due to the relative movement between the pore fluid and the soil skeleton [7; 8] will both directly affect the overall mechanical properties of the saturated porous media; in particular, the saturated sand is prone to liquefaction under dynamic loading [9]. Previous seismic investigations all indicated that ground failure and lateral large deformation due to soil liquefaction are the main causes of damage to the structure [10; 11; 12; 13; 14; 15; 16; 17]. The macroscopic mechanical properties of saturated porous media are much more complex than single-phase media, so it is necessary to establish the equations to describe the mechanical behavior of the soil skeleton and pore water, respectively, and to consider the fluid-solid coupling interaction. In 1956, [PERSON] [18] proposed a wave propagation theory of saturated porous media based on a quasi-static consolidation model, considering the inertial effects of the soil skeleton and pore fluid [19]. On that basis, [PERSON] et al. [20] rewrote the equations by using different simplifying assumptions, such as the _u-U-p_, _u-U_, _u-p_, and _u-w_ models. Among them, the _u-p_ formulation is a mixed equation of the vector displacement \(u\) and the scalar pore pressure \(p\). Compared with the full vector formulation of the dynamic equation (e.g., _u-U_), the _u-p_ equation has fewer degree of freedom and can be directly solved for the fluid pressure variables without the need of secondary calculations. It is suitable for describing the dynamic problems of saturated porous media under low-frequency loading. Usually, the dynamic equations of saturated porous media can be solved by the analytical method only under a few special boundary conditions and load forms [21]. Therefore, numerical methods represented by finite element method (FEM) are the main methods for solving the dynamic problems of two-phase media. In the time domain, the integration methods can be categorized into implicit methods and explicit methods, such as the Newmark method [22], the Wilson-\(\theta\) method [23], the generalized-\(\alpha\) method [24; 25], and the central difference method. The \(\mathbf{u}\)-\(\mathbf{U}\) and \(\mathbf{u}\)-\(\mathbf{w}\) formulation are often solved by using the staggered explicit-explicit [26; 27], explicit-implicit [28; 29] and implicit-implicit [28] methods. However, the implicit-implicit or explicit-implicit scheme is basically used to solve the \(\mathbf{u}\)-\(p\) formulation of saturated porous media. This is because in the \(\mathbf{u}\)-\(p\) formulation, the fluid phase equation is formally a diffusion equation in terms of pore pressure \(p\) and is usually solved by an implicit integration method. For example, [PERSON] [30] and [PERSON] et al. [31] established unconditionally stable staggered implicit-implicit solution procedures for the \(\mathbf{u}\)-\(p\) dynamic equations of two-phase media. [PERSON] and [PERSON] [32] improved the Zienkiewicz method by introducing a pore pressure correction method, which obtains both displacement and pore pressure through the estimated-corrected value calculation. Since the \(\mathbf{u}\)-\(p\) formulation is a mixed form of dynamic equations, using the same order of shape functions of displacement and pore pressure would lead to instability of the results. To overcome this problem, [PERSON] et al.