The fully dynamic (FD) model based on Biot’s theory and the soil elastic constitutive relation are used to investigate seabed oscillatory responses. This model is frequently employed to simulate the dynamic behavior of soil, accounting for both the acceleration of the soil skeleton and the pore fluid. The governing equation of the FD model can be expressed as follows (Zienkiewicz et al., 1980):

The variables are defined as follows: σ_ij represents the total soil stress (Pa), ρ stands for the average density of the porous medium (kg/m³), ρ_f denotes the density of pore fluid (kg/m³), g_i represents the gravitational acceleration in the i-direction (m/s²), u_i is the soil displacement in the i-direction (m), w_i signifies the relative displacement of the fluid respect to the solid skeleton in the i-direction (m), k_i represents the soil permeability in the i-direction (m/s), and n stands for the porosity of the solid phase. β in Equation (3) is the compressibility of the pore fluid. The expression of β is:

Where K _f is the bulk modulus of the pore fluid (N/m²), S _r is the saturation of soil. The shear stress of the soil under the action of waves can be obtained by solving the oscillatory response of the seabed, which is adopted as a driving force for the seabed residual response (Seed and Rahman, 1978). More details about the oscillatory model which involves the stress-strain relation, discretization scheme et al., interested readers can refer to the authors’ previous publication of Sui et al. (2019a).

According to Jeng (2012), the governing equations for the seabed residual response in a 3D space can be written as (Equation 5)

where p_r designates the residual pore pressure (Pa), C_v is the coefficient of soil consolidation (m²/s), and f denotes the source term for the accumulation of pore pressure (Pa/s). Here, C_v and f are defined as (Equations 6 and 7)

where E represents the Young’s modulus of soil (Pa), γ_w is the unit weight of water (N/m³), T is the period of wave loading (s), while α _r and β _r are empirical parameters that can be determined using the following expressions (Sumer et al., 2012) (Equations 8-10):

Here, τ_xz , τ_yz , and τ_xy are the shear stresses of the soil in space (Pa) (Equation 12), which are obtained from the oscillatory module of the present model. It should be noted that the saturation degree (S_r ) can have an influential role in pore pressure accumulation. This impact is included in the oscillatory model in section 2.1, as seen in Equations 1-4. The oscillatory model provides the shear stress for the source (f) of the residual model.

To address the governing equations mentioned above, it is essential to impose appropriate boundary conditions. In this model, the lateral boundary and the bottom boundary of the seabed are treated as rigid and impermeable. Consequently, the seabed displacements and the normal gradients of pore pressure are constrained to zero (Equation 13).

At the interface between the seabed and water, the residual pore pressure is maintained at zero, while the oscillatory pore pressure equals the dynamic wave pressure. Furthermore, the vertical effective normal and shear stresses become negligible at this juncture (Equation 14):

According to the linear wave theory, the dynamic wave pressure exerted on the seabed surface can be expressed as (Equation 15):

where H represents the wave height (m), d denotes the water depth (m), α stands for the wave number (m^-1), and ω is the wave frequency (s^-1).

To ensure the model’s reliability, it is essential to conduct verification before its practical application. Three validation cases are performed, involving the comparison of model results with previous analytical solutions and experimental data. The wave and seabed parameters for the validations are listed in Table 1 .

For the first validation, the proposed model is subjected to comparison with an analytical solution for the oscillatory seabed responses as previously presented by Jeng (2012). In this case, the wave is obliquely incident at an angle of 45°. The distributions of oscillatory pore pressure and soil stresses under short-crest wave loading are shown in Figure 1 . The present numerical result and the previous analytical solutions agree relatively well. Note that, the present numerical model adopted the same governing equations, boundary conditions, wave loading etc. as the analytical model. The good comparison shown in Figure 1 proves that the discretization method and numerical scheme used in the numerical model were appropriate and reliable for the present simulation.

For the second validation, the reliability of the proposed model in forecasting pore pressure accumulation is assessed using residual pore pressure data obtained from a wave flume test conducted by Sumer et al. (2012). In this test, pressure transducers were positioned at varying depths within a silt pit. Figure 2 presents a comparative analysis of the vertical distribution of residual pore pressure with respect to seabed depth, contrasting the findings from the present numerical simulation with both the numerical simulations and experimental data from Sumer et al. (2012). The proposed model predicted the accumulation of pore pressure relatively well. It is noted that, due to the mismatch of timescales between wave propagation and soil consolidation, scale effects exist when expanding the standard wave flume experimental results from laboratory scaling to field scaling. The present model is based on the theory of Seed and Rahman (1978) where the important coefficients α and β were obtained from the large-scale cyclic shear test by Peacock and Seed (1968), and Alba et al. (1976). Therefore, the present model is capable of simulating the seabed residual response on a large scale. Figure 2 shows that the present model can also simulate the residual response in a standard flume test on a small scale as well. This is because that the alterations in soil behavior with depth appear to have limited significance in relatively shallow soil layers, precisely where the primary occurrence of wave-induced pore pressure buildup and liquefaction takes place (Sumer, 2014).

The third validation pertains to the assessment of the wave-induced response in a layered seabed. Hsu et al. (1995) established a semi-analytical solution for pore pressure and soil stresses in a layered seabed subjected to dynamic wave loading. Figure 3 illustrates a comparison between the simulated results (represented by solid lines) and the analytical results (labeled data) with respect to the vertical distribution of pore pressure and soil stress in a layered seabed. As shown in Figure 3 , the simulation results and analytical solution agreed relatively well. The reason for the error may be that the present model retains the spatial gradient terms of the soil nonhomogeneous parameters in the governing equations, which were not considered by Hsu et al. (1995).

The spatial distribution of pore pressure for a “standard” case (using all the default values in Table 1 ) at (a) t = 40T and (b) t = 100T are shown in Figure 5 . The interface between the upper and lower layers is delineated by the black dashed line. The figure illustrates that the pore pressure in the upper layer (located above the dashed line) exceeded that in the lower layer (situated below the dashed line). In addition, comparing Figures 5A, B reveals that, when the time lasts from t = 40T to t = 100T, the pore pressure in the upper layer increases significantly, varying from approximately 1 × 10⁴ to 2.5 × 10⁴ Pa. In the lower sand layer, the pore pressure exhibits almost no change, with a value of approximately 0.5 × 10⁴ Pa. This indicates that the pore pressure increased significantly in the upper silt layer, whereas it hardly increased in the lower sand layer. This is because the upper silt part has partial drainage conditions, which makes it more capable of pore pressure buildup. However, the lower sand layer has relatively high permeability, which means that pore pressure dissipation is easier than pore pressure accumulation. Another reason is that the shear stresses in the soil of the upper silt layer are much larger than those in the lower part, which provides more “power” for the accumulation of the pore pressure (which will be discussed in detail in the next section).

Figure 6 presents the time series of pore pressure accumulation at a depth of z = -15 m within the lower sand layer for various upper silt layer depths. It is evident that pore pressure gradually increases with time due to the generation of excess pore pressure. In Figure 6 , the influence of the upper silt layer’s coverage thickness on the pore pressure in the lower sand layer is depicted. Figure 6 shows that, during the period from 0 to 800 s, the pore pressure increased slightly if there was no silt layer coverage (h = 0 m) (black solid line). When the coverage thickness of the upper silt layer was h = 5 m, the pore pressure increased dramatically from zero to approximately 1.8 × 10⁴ Pa during the same period (red solid line). However, the maximum pore pressure then decreased significantly when the coverage thickness increased (h = 10 m for the red dashed line and h = 15 m for the blue solid line). This indicates an “increase–decrease” pattern with the increase in depth of the upper silt layer, which is more clearly shown in the nested graph of Figure 6 with quantitative plotting. Here, p_{r , max} indicates the pore pressure at t = 800 s, which was the end of the simulation. The nested graph of Figure 6 shows that there was a peak value (p_{r , max} = 1.8 × 10⁴ Pa) when the coverage thickness of the upper silt layer h was approximately 5 m.

The effect of the coverage thickness of the upper silt layer (h) on the vertical pore pressure distribution at t = 100T is illustrated in Figure 7 . In this plot, the coverage thickness h ranges from 0 to 9 m. The figure illustrates a rapid initial increase in pore pressure, followed by a substantial decrease in the vertical direction. Notably, in the absence of silt layer coverage (h = 0 m), minimal pore pressure accumulation was observed in the vertical direction. Moreover, the coverage thickness was found to exert significant effects on the vertical distribution of pore pressure. For the various coverage thicknesses (h), the vertical locations for the maximum pore pressures were roughly constant at approximately z = −0.15 h_total . However, the maximum values exhibited a significant discrepancy at different coverage thicknesses. The maximum pore pressure increased with an increase in coverage thickness for h varying from 0 to 5 m. Subsequently, the maximum pore pressure had a large decrease with increasing h. For instance, the maximum pore pressure p _max was 16 kPa for h = 3 m, 18.3 kPa for h = 4 m, and 25 kPa for h = 5 m, but it decreased to 16 kPa when h = 6 m.

The impact of the permeability (k) and Young’s modulus (E) of the upper silt layer on the vertical distribution of pore pressure is depicted in Figures 8A, B . The dashed line demarcates the interface between the silt and sand layers. Figures 8A, B demonstrate that pore pressure accumulation predominantly transpires above the dashed line, signifying that pore pressure accumulates more readily in the silt layer, while it is more challenging to generate in the lower sand layer. In addition, the maximum pore pressure (p_{r , max}) increases with decreasing permeability. It is reasonable to assume that decreasing permeability indicates a partial drainage condition of the seabed, which would promote an increase in pore pressure. Regarding the effects of Young’s modulus, as shown in Figure 8B , the maximum pore pressure basically increases with an increase in Young’s modulus, except for one case (E = 8 × 10⁶ Pa). For instance, the pore pressure at E = 8 × 10⁶ Pa is approximately 17 kPa, which is larger than that of 15 kPa at E = 6 × 10⁶ Pa.

In addition to soil parameters, water depth are also important for pore pressure accumulation. Figure 8C illustrates the vertical distribution of (a) pore pressures and (b) shear stresses in the soil at different water depths. As indicated in Figure 8C , the predominant pore pressure accumulation was observed in the upper silt layer. In addition, for the various water depths, the locations for the maximum pore pressure are all approximately z/h_total = −0.1. The maximum pore pressure increases significantly with decreasing water depth.

Pore pressure within the seabed arises from two primary physical mechanisms: pore pressure generation and dissipation (Seed and Rahman, 1978). The source term (f) in Equation (4) is accountable for pore pressure generation and encompasses crucial variables such as shear stresses τ_xz (for the 2D cases) and the initial effective soil stress (σ ₀ ^’) defined by Equation (11). Greater shear stress τ_xz indicates a higher potential for pore pressure generation (Jeng et al., 2007; Jeng and Zhao, 2015). Sassa and Sekiguchi (1999) firstly defined the “shear stress ratio (χ)” rather than the shear stress to serve as an index of the severity of wave loading, on which the pore pressure generation depends on. The shear stress ratio is defined as the ratio of shear stress and the soil initial effective stress (Equation 16):

The calculation of shear stress (τ) is detailed in the oscillatory mode of the model (see Section 2.1).

On the other hand, pore pressure dissipation is regulated by the normalized spreading parameter S, which characterizes the ability of pore pressure to propagate through the soil skeleton. Parameter S is expressed as (Equation 17) (Sumer, 2014).

Here, in the given equation, L represents the wave length (m), and T represents the wave period (s). This parameter is first defined by Sassa and Sekiguchi (1999) and Sassa et al. (2001) (named the partial drainage factor in their paper). As outlined by Sumer (2014), the impact of the spreading parameter on pore pressure transmission can be understood from a physical perspective: a higher value of S implies enhanced spreading and dissipation of excess pore pressure.

The vertical distribution of shear stress of soil τ_xz, the shear stress ratio χ and the effective normalized spreading parameter S_e with various coverage thicknesses of the upper silt layer are shown in Figures 9A–C . The effective normalized spreading parameter S_e is newly defined in this study, which considers the effects of the upper-layered silt on the values of S above the point of concern. A detailed calculation of S_e can be found in the Appendix ( Figure 10 ). Figure 9A demonstrates that, akin to the vertical distribution of pore pressures, the shear stress τ_xz initially rises but subsequently declines with increasing seabed depth. The obvious disconnection in the curved line of the shear stress distribution is caused by the sudden change in the soil properties at the interface between the silt and sand layers. This disconnection was likewise observed in the pore pressure distribution, as depicted in Figure 7 , and it is more pronounced in the shear stress distribution presented in Figure 9A . This observation suggests that the potential for pore pressure generation (τ_xz ), diminishes at both the surface and bottom of the seabed, while exhibiting a peak value at the midpoint. Additionally, the maximum value of shear stress τ_xz increased as the coverage thickness followed by a decrease with an increasing h. This is also similar to the basic conclusions obtained on pore pressure variation ( Figure 7 ). Figure 9B illustrates the effects of the coverage thickness h on the maximum shear stress ratio χ. Note that, the shear stress ratio χ could be seen as the “power” for the generation of the pore pressure (Sassa and Sekiguchi, 1999). This means that a larger χ corresponds to a greater likelihood of pore pressure accumulation. Figure 9B shows that, in the vertical direction, the maximum shear stress ratio χ greatly decreases with the increase of the seabed depth. For an increasing coverage thickness h, there is a decrease in the maximum shear stress ratio. The above simulation results indicate that the pore pressure at the upper seabed layer has more power to accumulate with a smaller coverage thickness h. The effects of coverage thickness (h) on the vertical distribution of S_e are shown in Figure 9C . The effective normalized spreading parameter S_e exhibited an increase in value with the ascending seabed depth. This trend can be attributed to the significantly higher Young’s modulus E associated with the lower sand layer. This indicates that the pore pressure dissipation ability increases with seabed depth when the upper seabed has partial drainage. Moreover, Figure 9C illustrates the effective normalized spreading parameter S_e diminishes as the coverage thickness of the upper silt layer (h) increases. This relationship can be readily inferred from the definition of S_e , taking into account that the Young’s modulus E is lower for the upper silt layer. This implies that the ability of pore pressure dissipation decreases with an augmentation in coverage thickness (h).

Figure 11 shows the (A) the maximum pore pressure (obtained from Figure 7 ) and (B) the shear stress ratio χ, the normalized spreading parameter S_e with various thicknesses of the upper layer h. It shows that, the increase in the thickness of the upper layer (h) results in a rise of residual pore pressure until h reaches 5 m. After that, residual pore pressure exhibits a decreasing trend ( Figure 11A ). The reasons can be explained in Figure 11B for illustrating the changing trend of χ and S _e with various thickness of the upper layer. Figure 11B shows that an increase in the coverage thickness would cause a monotonically decreasing trend in χ and S_e . As per the definitions provided in Figure 9 , it can be inferred that an increase in the value of h results in a decrease in the shear stress ratio χ, which in turn leads to a potential decrease in the pore pressure. On the other hand, a decrease in the value of S _e corresponds to a potential increase in the pore pressure. The shear stress ratio (implying generation) and the normalized spreading parameter (implying dissipation) are the two important factors that govern the residual pore pressure. As shown in Figure 11A , the turning point for pore pressure variations occurred at a depth of 5 meters. This suggests that, for the current scenario, the normalized spreading parameter S_e has a dominant effect on the pore pressure when the depth is less than 5 meters. However, when the depth exceeds 5 meters, the dominant factor that affects the pore pressure is changed to the shear stress ratio χ.

To figure out the reason for the effect of permeability and Young’s modulus on pore pressure accumulation, Figures 12A–F illustrate the vertical distribution of τ_xz, χ and S_e with various permeabilities (k) and Young’s modulus E. As shown in Figures 12A, B , the shear stress and shear stress ratio did not change significantly when the permeability (k) decreased from 5 × 10⁻⁵ to 6 × 10⁻⁶ m/s. This indicates that the influence of permeability (k) on pore pressure generation is limited. Figure 12C shows that S_e decreases continuously with decreasing permeability. This inhibits pore pressure dissipation, resulting in an increase in the maximum pore pressure, as shown in Figure 8A . The effects of E on the vertical distribution of τ_xz , χ and S_e are shown in Figures 12D–F , illustrating that, in the upper silt layer, all of the τ_xz , χ and S_e have an increase with an increase in E. Detailed comments for this phenomenon will be proposed in the “Discussion”.

For the various water depths, only the factor for pore pressure generation is affected. The vertical distributions of shear stresses (τ_xz ) and shear stress ratio (χ) with various water depths are shown in Figures 12G, H . It is clear that the shear stress (τ_xz ) and shear stress ratio (χ) decrease with the increase in water depth (d). It is reasonable to assume that increasing the water depth would decrease the dynamic wave pressure at the seabed surface and weaken the soil response. This further decreases the shear stress (τ_xz ) and the shear stress ratio (χ) of the soil, which would correspondingly reduce the pore pressure generation capability. As a result, the maximum pore pressure decreases with increasing water depth, as shown in Figure 8C .

The assessment of wave-induced liquefaction was conducted with a meticulous examination to determine the stability of the seabed. It is widely acknowledged that liquefaction occurs when the excess pore pressure surpasses the initial effective stress (σ ₀ ^’), leading to a complete loss of soil resistance. In line with Jeng (1997)’s work, the criterion for the initiation of liquefaction is formulated as follows:

In this equation, γ′ represents the submerged specific weight of the soil (N/m³), and K₀ = 0.42 stands for the coefficient of lateral earth pressure.

Based on the liquefaction criterion (Equation 18), the liquefaction depths for various coverage thicknesses in the upper silt layer are shown in Figure 13 . The dashed line in Figure 13A corresponds to the overburden pressure of the soil in the vertical direction. The point of intersection between the pore and overburden pressures designates the liquefaction depth (z_L ). As depicted in Figure 13 , the liquefaction depth is zero when h = 0, as there is minimal accumulated pore pressure in the absence of an upper silt layer ( Figure 7 ). However, as the coverage thickness increases, the liquefaction depth also rises, reaching z_L = 4.5 m when h = 5 m. Subsequently, the liquefaction depth decreases significantly with increasing coverage thickness. This trend of pore pressure variation with coverage thickness is quantitatively shown in Figure 13B for clarity.

Figure 14 depicts the influence of permeability (k), Young’s modulus (E) of the upper silt layer, and water depth (d) on the liquefaction depth. Same as Figure 13A , the dashed lines in Figures 14A, C, E mark the soil initial effective stresses. It is evident that the liquefaction depth substantially diminishes with an increase in the permeability of the upper silt layer, as displayed in Figure 14A . An increase in permeability is anticipated to facilitate the dissipation of accumulated pore pressure, leading to a reduction in the liquefaction depth. In addition, Figure 14B shows that increasing the Young’s modulus first increases the liquefaction depth to the peak (z_L = 3.85 m) and then decreases it. This phenomenon is attributed to the non-monotonic impact of Young’s modulus on pore pressure, as demonstrated in Figure 8B . Furthermore, Figure 14C illustrates the effects of water depth (d) on the liquefaction depth of a layered seabed, where it is evident that the liquefaction depth decreases with increasing water depth. In practical engineering contexts, this implies that liquefaction is more likely to occur in relatively shallow water.