Effects of principal stress rotation on the fluid-induced soil response in a porous seabed

Principal stress rotation (PSR) is an important feature for describing the stress status of marine sediments subject to cyclic loading. In this study, a one-way coupled numerical model that combines the fluid model (for wave–current interactions) and the soil model (including the effect of PSR) was established. Then, the proposed model was incorporated into the finite element analysis procedure DIANA-SWANDYNE II with PSR effects incorporated and further validated by the experimental data available in the literature. Finally, the impact of PSR on the pore-water pressures and the resultant seabed liquefaction were investigated using the numerical model, and it was found that PSR had a significant influence on the seabed response to combined wave and current loading.


Introduction
Recently, the physical processes of fluid-seabed interactions have attracted great attention from coastal and geotechnical engineers because of the growth in human exploration and development of offshore projects. Seabed instability due to cyclic loading, such as waves, currents, and earthquakes, is one of the main concerns of offshore geotechnical engineers involved in the design of offshore infrastructure.
It has been well known that dynamic wave pressure generated by natural hydrodynamic loading on the sea floor further induces pore-water pressure and stresses in the seabed. When the pore pressure accumulates and reaches a certain level, the effective stresses vanish and lead to soil instability as a consequence of the movement of soil particles [1,2]. Therefore, an accurate prediction of the soil response, including pore-water pressure, effective stresses, and shear stresses, is important for the design of offshore infrastructure.
On the basis of laboratory and field measurements, two mechanisms for the wave-induced soil response have been developed and reported in the literature [3][4][5], namely oscillatory and residual mechanisms. The first mechanism is the result of oscillatory excess pore-water pressures and accompanied by the attenuated amplitude and phase lag of pore-water pressure changes [6]. The second mechanism is the build-up of excess pore pressures caused by the contraction of cyclic loading [7]. As reported in Jeng and Seymour [8], the oscillatory mechanism dominates the process of liquefaction in the case of a longer wave period and small amplitude, while the residual mechanism dominates the whole process for a wave with a short wave period or large wave amplitude.
Numerous studies for the wave-induced oscillatory soil response have been carried out since the 1970s. For example, on the basis of Biot's consolidation theory [9], Yamamoto et al. [6] derived a closed-form analytical solution for the wave-induced oscillatory soil response in an infinite seabed. The scope of this framework has been further extended to a seabed of finite thickness, a layered seabed, an inhomogeneous seabed with variable permeability, and a cross-anisotropic seabed [5]. Later, several analytical solutions were proposed that incorporated dynamic soil behavior; these models include the partial dynamic (u − p) [10] and full dynamic (FD) models [11,12]. Jeng and Cha [11] investigated the applicable range of different approximations with two non-dimensional parameters and soil permeability. This applicable range was reexamined by Ulker and Rahman [12] for different soils.
In addition to analytical approximations, several numerical models for the wave-induced oscillatory soil response for more complicated cases have been developed and applied to different offshore infrastructures. For example, Jeng and Lin [13] established a finite element model (FEM) that considers the effects of variable permeability and the shear modulus. Later, FEM models were developed for cases with different offshore infrastructures, for example, breakwaters [14], pipelines [15], and mono-pile foundations [16].
In the literature, numerous investigations of the wave-induced residual soil response are available. Using the results of direct shear tests [17], Seed and Rahman [7] proposed a 1D approximation with a source term for pore-water pressure generation. Following this framework, several analytical solutions and numerical solutions for wave-induced residual liquefaction were proposed [4,5,8]. Recently, Jeng and Zhao [18] proposed a new definition of the source term by considering the instant oscillatory shear stress; then, they extended the 1D model to two dimensions and applied it to the case of a submarine pipeline [19]. The aforementioned works were based on an inelastic model with a source term. Adopting the model proposed by Sassa et al. [20] and including the dissipation of pore-water pressures in the source term, Liao et al. [21] extended the model to two-dimensional cases. In addition to the inelastic models with a source term, a poro-elastoplastic model (DIANA-SWANDYNE II) was established by Chan [22] for earthquake-induced liquefaction by adopting the Pastor-Zienkiewicz Mark-III (PZIII) model [23]. This model was modified and applied to the problem of wave-seabed interactions around marine infrastructures, such as pipelines and breakwaters [24,25].
In natural ocean environments, the co-existence of waves and currents is a common physical phenomenon, and their interaction is an important topic in coastal and ocean engineering. The presence of a current in propagating waves directly changes the flow field and causes further changes to the soil response. On the basis of the analytical solution for wave-current interactions [26], Ye and Jeng [27] were the first to investigate the wave (current)-induced oscillatory soil response in a porous seabed. Following a similar framework, Wen et al. [28] further considered the case of a submarine pipeline by using the commercial software ABAQUS. Several analytical solutions based on different soil behaviors, such as quasi-static, partial dynamic, and full dynamic models, have been developed to describe the soil response to combined wave and current loading [29]. All of these approaches are based on the third-order analytical approximation for wave-current interactions [26]. Using the numerical model for wave-current interactions [30], Zhang et al. [31] further investigated the wave (current)-induced oscillatory pore pressures in a porous seabed. Recently, Liao et al. [32] proposed coupling models for residual seabed liquefaction subject to combined wave and current loading. Although numerous theoretical studies have been carried out since 2012, only two experimental studies for the wave (current)-induced oscillatory soil response are available in the literature [33,34]. These studies were used for the validation of the present model.
None of the aforementioned investigations have considered the effects of principal stress rotation (PSR) in a marine seabed, although the continuous rotation of principal stresses is an essential feature of soil's dynamic response to cyclic loading. Unfortunately, because pure PSR is assumed, this process cannot be captured by a conventional elastoplastic model without changing the cyclic deviatoric stress amplitude of the plastic strain [35]. Several experimental results have confirmed that plastic strains are generated merely by altering the principal stress orientation in both monotonic and cyclic rotational shear tests [36,37]. On the basis of the generalized plasticity theory Zienkiewicz and Morz [38], as the first attempt, Sassa and Sekiguchi [35] developed a modified version of PZIII model by considering the effects of principal stress orientation. Their model defines a new major principal stress angle parameter (Φ) that replaces potential plastic functions, the loading functions, and the plastic modulus. However, as Zhu et al. [39] pointed out, Sassa's model [35] also has deficiencies; for example, it does not account for out-of-plane stress, which is a critical parameter in the determination of plastic flow conditions. Furthermore, the reloading effect is not considered in their model. Recently, Zhu et al. [39] proposed a modified constitutive model in which both the PSR and the out-of-plane stress are taken into consideration within the generalized plasticity theory framework. In contrast to Sassa's model [35], this model was built to consider previous events during the reload by adding a discrete memory factor, and stress invariants were added at the same time to complete the optimization of the model. However, Zhu et al. [39] only considered linear wave theory in their model. The effects of PSR on the soil response to combined wave and current loading have not been reported in the literature.
In this paper, the constitutive model proposed by [39] is adopted, and the impact of PSR is included to examine the wave (current)-induced soil response in a sandy seabed. The theoretical model for both the flow model (wave-seabed interactions) and the seabed model (with PSR) are outlined first. The validation of the present model by both wave flume tests [33,34] and centrifugal tests [40] is then described. Finally, the results of the parametric study are reported to examine the effects of PSR with combined wave and current loading.

Theoretical Models
The present model consists of two submodels: flow and seabed submodels. A one-way coupling between the two different models is employed by the pressure continuity at the interface between the fluid and seabed domains. The fluid model is used to obtain the flow characteristics, such as wave motion, velocity field, wave pressures, etc. In the present model, the continuity of pressures is used to link the two submodels; that is, the dynamic fluid pressure is extracted and interpolated on the grid points of the solid model interface and serves as the pressure boundary condition for the seabed model. This approach has been commonly used by previous researchers [6,12,14,41,42].
The present study is based on the one-way coupling approach. Although one-way coupling has been widely used in the past and may still be effective in some cases, there are more recent approaches that effectively represent the water and bottom sediment coupled dynamics [43][44][45][46]. For example, Ran et al. [43] proposed an incompressible smoothed particle hydrodynamics (ISPH) scour model for movable bed dam break flows. Wang et al. [44] presented an ISPH simulation of scour behind seawall due to continuous tsunami overflow. Manenti et al. [45] adopted SPH model to investigate the Vajont disaster and compare their numerical simulation with 2D experiments. Wang et al. [46] further adopted their model [44] to 3D ISPH erosion model for flow passing a vertical cylinder. The technique could be further adopted to the present problem in the future.
The problem considered in this study is depicted in Figure 1. In the computation domain, the seabed thickness is h, and the water depth is d. The ocean wave propagates in the x-direction, while the z-axis is oriented upward from the seabed surface. The direction of the current can be the same as or opposite of the direction of wave propagation.

Flow Model
Recently, the open-source code OpenFOAM has become widely used for the simulation of various coastal/ocean engineering problems; for example, Waves2FOAM and IHFOAM have been used to study wave generation [47,48], wave-structure interactions, and other coastal engineering processes [49]. In this study, IHFOAM was adopted to describe the wave-current interactions. Basically, IHFOAM solves three-dimensional Volume-Averaged Reynolds-Averaged Navier-Stokes (VARANS) equations for two incompressible phases (water and air) using a finite volume discretization and volume of fluid (VOF) method. The governing equations, including mass conservation and momentum conservation equations, can be expressed as ∂ρ u f i ∂t where and f are Darcy's volume-averaging operator and the intrinsic averaging operator, respectively; ρ is the density, computed by ρ = αρ water + (1 − α)ρ air , in which α is the indicator function defined in (4); u f i is the velocity vector; n is the porosity; p * is the pseudo-dynamic pressures; g i is the gravitational acceleration; µ e f f is the efficient dynamic viscosity, defined as µ e f f = µ + ρν turb , in which µ is the molecular dynamic viscosity and ν turb is the turbulent kinetic viscosity, given by the chosen turbulence model. In this study, the k − turbulence model is used. The last term in (2) represents the resistance of porous media and can be expressed as where the factor C is less significant than factors A and B (refer to [50] for the values). A value of C = 0.34 [kg/m 3 ] is often applied by default [51]. Each cell in the computational domain is considered a mixture of a two-phase fluid (air and water). The indicator function α varies from 0 (air) to 1 (water); α is defined as the quantity of water per unit of volume for each cell and calculated as follows: Any variation of fluid properties, such as density and viscosity, can be represented using the indicator function α considering the mixture properties: where Φ water and Φ air are water and air properties, respectively, such as the density of the fluid. The fluid's movement can be tracked by solving the following advection equation [52]: ∂α ∂t where |u f c | = min c α |u f |, max(|u f |) , in which the default value of c α is 1.
The wave generation and active wave absorption in the fluid domain were implemented within IHFOAM. Several boundary conditions were introduced: (i) the inlet boundary condition allows for generating a wave according to different wave theories as well as adding different steady current flows; (ii) the outlet boundary condition applies an active wave absorption theory to prevent the re-reflection of an incoming wave; (iii) the slip boundary condition (zero-gradient) is applied on the bottom of the fluid domain and the lateral boundary of the numerical wave flume; (iv) the top boundary condition is set as the atmospheric pressure. For the details of IHFOAM, the readers can refer to Higuera et al. [48].

Seabed Model
In the literature, three different models of fluid-seabed interactions have been established on the basis of different soil behaviors: quasi-static (QS, i.e., the conventional consolidation model), partial dynamic (u − p), and full dynamic (FD) models. All three are based on Biot's porous theory [9,53]. Zienkiewicz et al. [54] proposed the u − p approximation and examined the applicable range between u − p and QS for earthquake loading. The framework has been further extended to the problem of wave-seabed interactions [10][11][12]. The applicable range of different models has been clarified for various soil types [11,12].
This paper establishes a two-dimensional model that considers the rotation of the principal stress axis to analyze the seabed response to combined dynamic loading due to wave-current interactions. The dynamic Biot's equation proposed by Zienkiewicz et al. [54], the u − p approximation, was adopted. ∂σ x ∂x ∂τ xz ∂x where σ x and σ z are the effective normal stresses in the xand z-directions, respectively; τ xz is the shear stress; p e is the pore-water pressure; u s and w s represent the soil displacement in the xand z-directions; b g is the gravitational acceleration; K s is the soil permeability; ∇ 2 is the Laplace operator; n is soil porosity; ρ is the average density of a porous seabed and defined by ρ = ρ f n + ρ s (1 − n), in which ρ f is the fluid density while ρ s is the solid density. In (9), the compressibility of the pore fluid β s is defined as [55] where K w is the true bulk modulus of the elasticity of water (which can be taken as 1.95 × 10 9 N/m 2 [6]), S r is the degree of saturation, and P w0 is the absolute water pressure. When the soil is fully saturated (i.e., it is completely air-free), then β s = 1/K w since S r = 1.
The anisotropic elastic constitutive model cannot account for the directional effect of the principal stress or the dilatancy of sand. Compared with the elastic constitutive model, plastic constitutive models can more realistically simulate the stress-strain relationship of soil under dynamic load conditions and the accumulation of pore-water pressure. Therefore, Zhu et al. [39] proposed a plastic constitutive model, which was implemented in the DIANA-SWANDYNE II [22] finite element code. This code is used to analyze the seabed response of the principal stress axis to waves and ocean currents. In Zhu's plastic constitutive model [39], the loading direction vector n ij = ∂ f /∂σ ij , the plastic flow direction vector m ij = ∂g/∂σ ij , and plastic modulus H L (p , q, θ, ψ) are defined. For the theory of generalized plasticity, the loading function f (p , q, θ, ψ) and plastic potential function g(p , q, θ, ψ) do not need to be explicitly defined.
The elastic-plastic constitutive matrix can be expressed in tensorial notation as where D e ijkl is the elastic stiffness tensor. The loading direction vector n ij and the plastic flow direction vector m ij can be defined as and where χ is he control parameter to account for the effect of PSR; ψ is the major principal stress angle; α 0 , a, and c are the principal stress orientation model parameters; and M f c and M gc are model parameters related to the stress ratio. In addition, taking into account the effects of PSR, the plastic modulus H L is given as where H 0 , β 0 , and β 1 are model parameters, and where ξ is the accumulated deviatoric plastic strain.
To consider the history of loading events throughout the reloading process, the discrete memory factor H DM is introduced in the following: and γ d is the coefficient for the discrete memory factor. The plastic modulus H U for unloading is where H U0 and γ u are original model parameters.

Model Verification
To validate the present model, two comparisons with previous experimental data are presented here. First, we compare the present model with the hollow cylinder apparatus (HCA) element tests [56] for the present constitutive model. Second, we compare the present model with the wave flume test for wave (current)-induced oscillatory pore-water pressures [33,34]. Third, we compare the present model with centrifugal tests [40] for the development of pore-water pressure build-up.

Comparison with Hollow Cylinder Apparatus (HCA) Element Tests
In the validation of the constitutive model with regard to PSR, the HCA elementary test is urgent in need due to its particular ability in simulation of the PSR through altering its control parameters such as axial load, torque, inner cell pressure, and outer cell pressure. Towhata and Ishihara [56] such a test with pure rotation of principal stress axis, in which the major principal stress orientation angle ranged from −π/4 to π/4 with a constant deviatoric stress of 76.7 kPa. The main model parameters are shown in Table 1. The predicted results of the present model with PSR are plotted in Figure 2, in which the test data of Towhata and Ishihara [56] are also illustrated for comparison. It can be clearly seen that the adopted constitutive model present behaves well in capturing the effect of PSR on the development of volumetric strains with the constant amplitude of deviatoric stress. Moreover, the feasibility of the constitutive model with PSR is well validated by the good agreement between the predicted results and the measured data as shown in Figure 2.  [56].
• HCA test results The present model  Comparison of the present model with the HCA tests [56] under continuous rotation of the major principal stress axis.

Comparison with Laboratory Experiments for the Seabed Response to Waves and Currents
The second verification of the present model compares the present model's results with experimental results. Qi and Gao [33] conducted a series of experiments to investigate the seabed response to different wave and current velocities around a single pile. In their experiments, a wave flume of 52 m × 1 m with a depth of 1.5 m was set over sandy soil. The depth of the sandy tank was 0.5 m. However, the data for the wave and opposite current cannot be adopted for the comparison with the numerical model because of the effect of the single pile, which was set in the middle of the experiments. Thus, only the cases of the wave alone and the wave with a following current (when waves and currents have the same direction of propagation) were used to verify the model. In their experiments, the pore-pressure build-up (i.e., residual mechanism) was not observed, i.e., the tests are in the range applicable to the oscillatory mechanism, for which the elastic model is used for comparison here.
Using the same conditions as those of the wave flume tests [34], the following values are assumed: the water depth is 0.5 m; the soil model is set below the wave flume with a length of 2.4 m; and the thickness of the saturated soil layer is 0.5 m. The wave profile at the free surface and the corresponding pore-water pressure beneath 0.1 m of the seabed surface are compared. The numerical results of the wave profile and pore pressures have an overall agreement with the experimental data, as shown in Figures 3 and 4.
In the above comparisons, the top subfigures show the water surface elevation and the bottom subfigures show the pore-water pressures at 10 cm beneath the seabed surface. As shown in Figures 3 and 4, the present model predicts the wave profiles (η) well, but the pore-water pressures present some differences between the numerical results and experimental data, although the general trends are in agreement. The pressent model Laboratory experiments [33] (a) Water surface elevation (η) The pressent model Laboratory experiments [33] (b) Wave-induced pore pressure in seabed (p s )  The pressent model Laboratory experiments [33] (a) Water surface elevation (η) The pressent model Laboratory experiments [33] (b) Wave-induced pore pressure in seabed (p s )

Comparison with Centrifuge Tests and Previous Numerical Model for the Seabed Response To Waves
Sassa and Sekiguchi [40] carried out a series of geotechnical centrifuge tests to investigate the process of wave-induced liquefaction in a sandy seabed. To verify the present model, we reproduced the experimental conditions and compared the results with the centrifugal experimental data [40]. Their wave tests were all performed with a centrifugal acceleration of 50 g (where g is the gravitational acceleration). The soil bed was 200 mm in width and 100 mm in depth. The submerged unit weight of soil, γ , was equal to 425 kN/m 3 , and the wave number, κ(= 2π/L 0 ), was 12.2 m −1 . The corresponding wave loading intensity p 0 and cyclic stress ratio (χ 0 = κ p 0 /γ ) were 5.0 kPa and 0.14, respectively. Other input data are listed in Table 2. In the numerical model, we converted the problem back to an environment with a gravitational acceleration of 1. As shown in Figure 5, in general, the present model has an overall agreement with the centrifuge tests [40]. By examining the comparisons closely, we observe that the present model can capture the magnitude of the maximum pore-water pressures and the time it takes to reach the maximum pore water pressures. However, there are some differences between the numerical prediction and the centrifugal tests that occur during the pore-water pressure build-up. This implies that the present model requires further improvement. However, the magnitude of the maximum pore pressures directly affects the liquefaction depth, which is more important for practical engineering design. Therefore, the present model can provide sufficient information for engineering design. In addition to the comparison with the centrifugal tests [40], the present model is also compared with Sassa's numerical model [35] in Figure 6. The phenomenon of pore pressure build-up can be observed in the first several wave cycles. After a certain wave period, p residual reaches its peak value and stabilizes because liquefaction occurs. In the present model, the calculated time it takes for the residual pore pressure to reach its peak is 1100 µs, which is less than that predicted by Sassa and Sekiguchi [35] (1200 µs). This is because Sassa and Sekiguchi [35] did not consider the effect of out-of-plane stress, which is important for determining the plastic flow direction [57][58][59].

Results and Discussion
In this study, two new features were incorporated into an existing model for soil response: (1) PSR effects and (2) the combined wave and current loading. In this section, the seabed is considered to be an elastoplastic medium, and the discussed results are from simulations using the generalized plasticity model PZIII and modified PZIII model with PSR in the finite element analysis program DIANA-SWANDYNE II. Nevada dense sand was adopted for a seabed with elastoplastic behavior, and the soil parameters, given by Sassa and Sekiguchi [40], were determined experimentally and are specified in Table 3. In the computations, the seabed length L s = 180 m, and seabed thickness d = 30 m. In order to ensure the numerical convergence suggested by Ye et al. [42], in the numerical model, the maximum horizontal mesh size was less than the wavelength L w /40, where L w is 88 m, and the maximum vertical mesh size was half of the horizontal mesh size. Therefore, the horizontal mesh size and vertical mesh size were 1 m and 0.5 m, respectively. Furthermore, the time step ∆t was T/40, where T is the wave period and equal to 8 s in this study. Table 3. Parameters used in dynamic constitutive model for parametric study.

Parameters Original PZIII The Present Model Unit
40,000 40,000 kPa

Seabed Liquefaction
Generally, the literature reports two different mechanisms of fluid-induced soil liquefaction [3]: momentary liquefaction and residual liquefaction. Momentary liquefaction normally occurs near the wave trough and in an unsaturated seabed when the upward seepage force is higher than the overlying pressure. However, the effect of momentary liquefaction is much smaller than that of residual liquefaction on the stability of offshore structures. As mentioned previously, the soil gradually loses stability as the wave spreads over the surface of the seabed. When the pore-water pressure increases, the effective stress between soil particles decreases. The condition for residual liquefaction occurs when the pore pressure reaches the maximum value and the effective stress approaches zero, at which point the soil loses its bearing capacity. The seabed's instability is a consequence of the horizontal or vertical movement of soil particles [1]. In such a situation, the soil acquires the behavior of a liquid. The liquefaction of the seabed has an essential impact on the safety of offshore structures.
In order to quantitatively study the liquefaction characteristics of the sandy seabed involving PSR effects, a parameter called liquefaction potential is introduced, as defined below.
where σ zd is the wave (current)-induced dynamic vertical effective stress and σ z0 is the initial vertical effective stress.
With the liquefaction criterion proposed by Okusa [60], the sandy seabed liquefies at L p = 1. In practice, however, the value of L p may not reach 1. This is because sand is a non-viscous granular material and cannot withstand any tensile stress. Wu et al. [61] suggested that the adjustment coefficient α r should be 0.78-0.99 for liquefaction in sandy soils according to different soil characteristics. In this study, it is assumed that soil liquefaction occurs when the liquefaction potential reaches 0.9. Figure 7 shows the liquefaction zone of soils at different times for the same wave condition. The numerical results for two cases-with and without PSR taken into account-are included. As shown in the figure, after 10 wave cycles, the seabed in the original PZIII model just begins to liquefy, and the liquefaction depth is less than 1 m. However, the seabed in the PSR model is markedly liquefied, and the depth is about 4 m. After 20 cycles of continuous wave action on the seabed surface, the depth of soil liquefaction changes slightly in the original PZIII model. However, in the present model, the liquefaction depth increases from 4 meters to 8 m. At the end of the simulation, the soil liquefaction depth in the original PZIII model is about 2 m. However, when considering the effect of PSR, the liquefaction depth increases to 14 m in the present model. It can be inferred that PSR plays a vital role in the stability of the seabed, and it also has an essential influence on the liquefaction depth. If the influence of PSR is neglected, the likelihood of the liquefaction of a sandy seabed is seriously underestimated, which poses a significant threat to coastal engineering.   Please note that the present model does not consider the process of progressive liquefaction (i.e., post-liquefaction). This is why the predicted liquefaction depth in Figure 7 is large (up to 14 m). As reported in the literature regarding wave-induced post-liquefaction [20,62], the maximum liquefaction depth approaches a constant value when the concept of progressive liquefaction is taken into account. However, this concept is not included in the present model. Therefore, the predicted liquefaction depth at t/T = 30 may be overestimated. This indicates that the existing model requires further improvement.

Effect of Currents
Waves and currents usually coexist in the natural marine environment, and the effects of currents on the seabed cannot be ignored. Currents not only change the length and direction of wave propagation but also affect the stability of the seabed. In this section, to demonstrate the effects of currents with PSR on the seabed response, a velocity of 1 m/s for both following and opposing currents was added to the present model to compare the soil response with the principal stress axis rotation (subject to combined wave and current loading). The term "following current" means that waves and currents propagate in the same direction, and the term "opposing current" indicates the opposite direction of waves and currents. Figure 8 plots the variation in pore-water pressure and effective stress between soil particles for different current directions (following current U 0 = 1 m/s, no current U 0 = 0, and opposing current U 0 = −1 m/s). As seen from the figure, among the three cases, the soil for which the direction of current propagation is the same as that of wave propagation is liquefied first. When the current's direction is opposite of the wave propagation direction, the instability of soil liquefaction is effectively prevented. It is also noted that the pattern of the liquefaction zone is wave-like in appearance. A possible explanation is that the water particles move in the horizontal direction with combined wave and current loading, causing the wave pattern in the liquefaction zone. Unfortunately, there is no experimental evidence available to confirm these two findings. More detailed experimental works are required in the future. As shown in the figure, for the case with a following current, the seabed at z = −15 m is liquefied after 22 wave cycles, while liquefaction occurs after 28 wave cycles for the case without a current. However, the soil with reverse current experiences 30 wave cycles and does not show liquefaction.   Figure 9 shows the horizontal and vertical displacements of soil particles for different directional currents at a depth of 5 m. It can be seen from the figure that the lateral displacement of soil remains almost unchanged before soil liquefaction, but the vertical displacement increases continuously. After 13 wave cycles, the liquefaction of soils occurs for both cases (with following current and no current), and the transverse displacement increases rapidly with time. This trend is more evident in the presence of the following current. In the vertical direction, the vibration amplitude of soil displacement is gentle, but when the soil is liquefied, the vibration of vertical displacement is significant. This phenomenon shows that before soil liquefaction, soil particles are continuously compressed, and pore-water pressure gradually rises but does not dissipate. These conditions eventually lead to soil liquefaction when pore pressure is higher than the vertical effective stress. When the soil loses its stability after liquefaction, the transverse and vertical displacements change significantly.   Figure 10 illustrates the distribution of liquefaction potential when considering the impact of PSR in the vertical direction of soil for different current conditions. It clearly shows that the soil near the surface is more prone to liquefaction. Moreover, the depth of soil liquefaction gradually increases. It also shows that when the soil depth is the same, the liquefaction potential in the following current case is the most significant, so soil liquefaction occurs more easily. On the other hand, when the liquefaction potential is the same, the maximum soil liquefaction depth in the following current situation is greater than the others. In other words, this result shows that the following current accelerates the process of soil instability, making soil liquefaction more likely to occur. Reverse currents have an opposite effect, so they are conducive to soil stability.

Effect of Principal Stress Rotation with Various Wave and Soil Parameters
It is well known that when waves propagate on a porous seabed, wave parameters, including wave height and wave period, are closely related to the liquefaction of the seabed [4]. Generally, with a longer wavelength (L w ) and higher wave height (H), the seabed liquefaction depth is more obvious. Furthermore, different soil parameters, including soil permeability (K s ) and saturation rate (S r ), have a significant impact on seabed liquefaction. This section compares the effects of different wave periods, wave heights, soil permeabilities, and saturation levels on soil liquefaction with and without consideration of PSR. Four different time stages are used to illustrate the relationship between the vertical direction of soil and the liquefaction potential. Figure 11 shows the effects of wave height on the vertical distribution of wave (current)-induced liquefaction potential with and without consideration of PSR conditions. As can be seen from the figure, in both situations, the liquefaction potential becomes more significant with increasing wave height. Also, the depth of soil liquefaction gradually increases over time. When t/T = 10 and z/h = 0.2, the value of L P increases from 0.25 to 0.82 with an increase in wave height from 1 m to 3 m when considering PSR. However, the liquefaction potential value increases from 0.15 m to 0.2 m when the effects of PSR are not considered. Figure 12 illustrates the relationship between liquefaction potential and wave period with the change in soil depth. It can be seen in the figures that as the wave period increases, the value of the liquefaction potential increases and the soil is more liable to liquefy. This phenomenon is due to the increase in wavelength or wave height, both of which result in more energy. It can strengthen the interaction between the wave and seabed foundation. Also, at the same soil depth, the value of the liquefaction potential increases more significantly when PSR is taken into account. Therefore, it can be said that given different wave parameters, PSR increases the liquefaction depth and liquefaction potential of the soil. Figure 13 shows the impact of soil permeability on the liquefaction potential along the vertical direction in the sandy seabed. Similarly, there is a steady increase in the liquefaction potential with the wave propagation for 35 wave periods. Also, the results of K s = 10 −5 m/s and K s = 10 −7 m/s are almost indistinguishable. Compared with K s = 10 −2 m/s, when soil permeability is lower, the soil is more likely to liquefy. This is because the permeability of the soil is large, and the pore pressure between the soil particles dissipates rapidly during the wave propagation and does not increase cumulatively. Thus, the effective stress between soil particles is sufficient to maintain the stability of the soil, and this stability prevents the occurrence of soil liquefaction. When the permeability of the soil is relatively low, the pore-water pressure between soils does not dissipate efficiently along with wave propagation, resulting in the faster accumulation of pressure. Therefore, soil liquefaction easily occurs when soil permeability is small. When considering the existence of PSR, the liquefaction potential increases rapidly in the same situation. Therefore, it can be concluded that PSR has a considerable impact on soil liquefaction for different permeabilities.      Figure 14 presents the variation trend of the liquefaction potential with soil depth at different saturation levels at four typical time points. On the whole, the relationship between L P and different saturation levels shows the opposite trend, with a certain soil depth as the limit. Since the liquefaction potential is set near L p = 0.9, the region with low L p values effectively means that there is no liquefaction at all. Therefore, we only need to investigate the region near the seabed surface. In the surface layer of the soil, the liquefaction potential decreases with increasing saturation.

Conclusions
Principal stress rotation (PSR) is an important factor in the evaluation of wave (current)-induced seabed instability. In this study, a one-way coupled numerical model that incorporates a wave model and soil model was developed to investigate the effect of PSR on the wave (current)-induced dynamic response of an elastoplastic seabed foundation. The comparisons show that the proposed model agrees well with laboratory wave flume tests, geotechnical centrifuge tests, and previous numerical results. On the basis of the numerical examples, the effects of PSR and currents with various soil parameters were examined, and the following conclusions are drawn: (1) Principal stress rotation (PSR) has a significant effect on the soil liquefaction depth. It accelerates the growth of pore pressures and reduces the vertical effective stress, so that the soil is easier to liquefy. (2) The existence of ocean currents has an important impact on the development of the liquefaction potential of a seabed foundation. When considering the interactions between waves and currents, the soil pore pressure and effective force change significantly and have a significant impact on soil liquefaction. The following current aggravates the soil reaction and promotes soil liquefaction. On the contrary, the opposing current reduces soil instability and plays a positive role in soil stability.
With the combined action of waves and current, the seabed with porous media shows pronounced lateral expansion and vertical settlement. (4) The liquefaction potential of the elastoplastic seabed foundation increases with time and decreases with depth. This indicates that liquefaction is more likely to occur in the upper layer of the seabed foundation.
Please note that the above conclusions are based on the numerical examples presented in this manuscript, and comparable experimental data are not available in the literature. The above findings require further confirmation by experimental evidence in the future.
In this study, we adopted the model proposed by Zhu et al. [39] for the soil response to combined wave and current loading. Note that Liu et al. [63] proposed another model to modify the previous model [35], which can also be used for the present problem.