Wave-Induced Instantaneous Liquefaction of a Non-Cohesive Seabed around Buried Pipelines: A Liquefaction-Associated Non-Darcy Flow Model Approach

: In complex marine environments, the wave-induced instantaneous liquefaction of the seabed is a key issue for the long-term safety control of marine structures. Existing computational frameworks for instantaneous liquefaction result in unreasonable tensile stresses in a non-cohesive seabed. To address this issue, a liquefaction-associated non-Darcy flow model has been proposed, but it has only been applied to the scenario of a pure seabed without a structure. In this study, we applied the previously proposed non-Darcy flow model to investigate the mechanism of wave–seabed–structure interactions under extreme wave loading considering a pipeline fully buried in a non-cohesive seabed. By comparing the liquefaction depths in the presence and absence of structures, it was found that the existence of structures weakens the attenuation of the pore pressure amplitude and influences the overall pore pressure distribution. Parametric studies were conducted. It was found that the liquefaction depth from the non-Darcy model is approximately 0.73 times that from the traditional Darcy model, regardless of whether or not a pipeline is involved. A quantitative relationship between the wave loading and structural size was established. The liquefied zone above the buried pipeline was found to be smaller than that in a pure seabed without a structure. A tentative explanation is provided for this phenomenon.


Introduction
In the wave-dominated coastal zone, the wave-induced liquefaction of the seabed often occurs due to the propagation of waves inducing excess pore pressure within the seabed, leading to a decrease in effective stress.Seabed liquefaction has a close relation with the erosion or resuspension of marine sediment [1][2][3][4][5] and can affect the stability of coastal structures [6,7].Therefore, it has become a key concern within the scope of marine engineering.Based on laboratory experiments and field measurements [8,9], wave-induced seabed liquefaction can be categorized into instantaneous liquefaction and residual liquefaction, induced by instantaneous pore pressure and residual pore pressure, respectively.This study focuses on the instantaneous liquefaction occurring in a non-cohesive seabed.
However, the existing computational methods and analytical models associated with instantaneous liquefaction can result in unreasonable tensile stresses in a non-cohesive seabed [50,51].This nonphysical phenomenon further influences the overall distribution of pore pressure [50,51].In fact, when instantaneous liquefaction occurs in a non-cohesive seabed, the physical and mechanical parameters of the seabed should change accordingly.However, most existing models employ an invariant poro-elastic assumption [14,35].
In order to eliminate unreasonable tensile behavior, Zhou et al. [52] was inspired by the literature [53,54] and introduced a dynamic permeability model.This model uses soil permeability k s that is determined by pore pressure p instead of the constant permeability assumed in conventional Darcy's law.This treatment mitigates the occurrence of unrealistic tensile stresses in the liquefaction zone.Compared with the previous dynamic permeability models [53,54], the new model built by Zhou et al. [52] agrees with the concept of an increase in permeability during soil liquefaction [55][56][57][58][59][60].However, this model [52] faces challenges in nonlinear convergence.When using larger model parameters or finer computational grids, numerical divergence may even occur, limiting its applications.Therefore, Zhou et al. [61] further converted the instantaneous liquefaction problem into a mathematical problem, i.e., a nonlinear complementarity problem (NCP).They employed specific Karush-Kuhn-Tucker (KKT) conditions for instantaneous liquefaction and a Lagrange multiplier method to solve the problem, while they optimized it using the primal-dual active set strategy (PDASS) [62].Additionally, they improved computational efficiency by employing direct delta function interpolation to condense the multipliers.This method addresses the challenge of nonlinear convergence caused by the previous dynamic permeability model [52].However, solving the NCP problem falls within the scope of a constrained variational theory framework and is not a standard module in most finite element codes (such as PORO-FSSI-FOAM [20,21]).Applying the NCP approach to broader scenarios, such as pipelines that are partially buried in the seabed [16,23], gravity foundations for offshore wind turbines [32], and monopile foundations [33,35], is also not convenient.
Recently, Zhou et al. [63] proposed a new KKT condition for wave-induced instantaneous liquefaction through the variational equivalence principle; thus, a non-Darcy flow model was proposed for instantaneous liquefaction.This model not only eliminates the phenomenon of tensile stress in the liquefied zone but also provides significant convenience for numerical implementations, because only the formulation of a seepage constitutive equation is required to be incorporated into finite element codes.Furthermore, it provides a unified explanation for various existing phenomena, such as non-cohesive granular soils with an upward seepage flow [64]; experimental results on piping in sandy gravels [65]; field trials [3]; and DEM-PFV simulations [66].However, these works were still limited to the scenario of wave-seabed interactions without a structure.
In this work, we further took into account a pipeline fully buried in a non-cohesive seabed under wave loadings in order to investigate wave-seabed-structure interactions with the use of our non-Darcy flow model [63].First, we extended our previous finite element code to a scenario that can consider a structure fully buried in the seabed.This extension was validated by comparing our numerical results with the experimental results of Turcotte [67].Then, a computational example under extreme wave loading was conducted to compare the liquefaction depths in the presence and absence of structures.Parametric studies were further accomplished considering wave, seabed, and structural parameters.Non-dimensional parametric studies were also carried out.
The remainder of this paper is as follows.Section 2 provides a brief overview of the fundamental theory and governing equations of the non-Darcy flow model considering a fully buried structure.In Section 3, the wave-seabed-structure model is validated based on classical experimental results without liquefaction.Section 4 presents the numerical results under extreme wave loading conditions, which induce seabed liquefaction.A quantitative relationship between the wave loading and structural response is exploited.Section 5 tries to provide an explanation for the liquefaction phenomenon in a seabed including a buried structure.Finally, Section 6 draws several conclusions.

Wave-Seabed-Structure Model and Liquefaction Criteria
Figure 1 depicts the wave-seabed-structure model and the process of seabed liquefaction, where d t represents the thickness of the seabed, h represents the water depth, and H represents the wave height.The coordinate value on the seabed surface is 0, and positive values exist below the seabed.Ω L denotes the zones in the seabed where instantaneous liquefaction occurs.During the process of wave loading, when the soil is located at or near a wave trough, the wave-induced seepage force is upward.When the upward seepage force equals the initial vertical effective stress of the seabed soil, the effective stress of the soil becomes zero, leading to an instantaneous liquefaction state.It should be noted that the wave in Figure 1 is linear.If considering nonlinear waves, Stokes waves [17] or cnoidal waves [68] should be employed.To simplify the model, this study applies linear wave theory [69]: where P b , whose amplitude is represented by p 0 , represents the wave pressure acting on the seabed surface (z = 0).L denotes the wavelength, while T denotes the wave period.The values of p 0 and L can be calculated by the following equations: where γ w denotes the unit weight of water, while g represents the acceleration of gravity.
The governing equations of the wave-seabed-structure model are as follows: and β represents the pore-fluid compressibility, which is determined by [71] where 0 w K represents the true bulk modulus of pore water, which can be taken as 2.0 × Γ represents the Neumann boundary.In Equation (3d), u represents the vector of displacement, while û denotes the displacement that is constrained.In Equation (3e), σ denotes the tensor of total stress, σ n denotes the unit normal vector pointing outward from σ Γ , and t denotes the surface force per unit area acting on the boundary.In Equa- tion (3f), p represents the pore pressure.In Equation (3g), w v denotes the vector of pore-fluid velocity, v n denotes the unit normal vector pointing outward from v Γ , and ˆn w v represents the boundary Darcy velocity.
It should be noted that the research presented in this paper aims to extend the non-Darcy flow model from pure seabed scenarios without structures to scenarios with structures.The presence of a structure in the seabed implies increased boundary conditions compared to the pure seabed.In order to simplify the computational model, the boundary The first two equations are derived from the theory of poro-elasticity theory [70].Equation (3a) denotes the equilibrium of the solid-fluid mixture, in which σ ′ represents the effective stress, I 2×2 denotes the second-order unit tensor, and b denotes the body force per unit volume.Equation (3b) denotes the conservation of mass, in which ε v represents the volumetric strain, t represents time, k s is Darcy's coefficient of permeability, and β represents the pore-fluid compressibility, which is determined by [71] where K w0 represents the true bulk modulus of pore water, which can be taken as 2.0 × 10 9 Pa [14].S r denotes the degree of saturation, while P abs represents the hydrodynamic pressure of water.The remaining expressions in Equation (3c-f) denote boundary conditions.Γ u and Γ σ denote two boundaries of the solid phase, where Γ u represents the Dirichlet boundary and Γ σ represents the Neumann boundary, while Γ p and Γ v denote two boundaries of the fluid phase, where Γ p represents the Dirichlet boundary and Γ v represents the Neumann boundary.In Equation (3c), u represents the vector of displacement, while û denotes the displacement that is constrained.In Equation (3d), σ denotes the tensor of total stress, n σ denotes the unit normal vector pointing outward from Γ σ , and t denotes the surface force per unit area acting on the boundary.In Equation (3e), p represents the pore pressure.
In Equation (3f), v w denotes the vector of pore-fluid velocity, n v denotes the unit normal vector pointing outward from Γ v , and vn w represents the boundary Darcy velocity.It should be noted that the research presented in this paper aims to extend the non-Darcy flow model from pure seabed scenarios without structures to scenarios with structures.The presence of a structure in the seabed implies increased boundary conditions compared to the pure seabed.In order to simplify the computational model, the boundary condition around the structure is considered an impermeable boundary, corresponding to Equation (3), where vn w is taken as zero.When the seabed experiences instantaneous liquefaction, there are two widely applied liquefaction criteria for determining the liquefaction zone (Ω L ), which use the tensile region of the sandy soil as the liquefaction zone: where p is the excessive pore pressure caused by waves.Its gradient along the vertical direction is expressed as j z = ∂p/∂z.γ ′ denotes the buoyant unit weight of the seabed and obeys the following expression: where G s represents the specific gravity of sand particles, and n represents the porosity of the sand, which is related to the void ratio e: n = e/(1 + e).

Numerical Implementation of the Non-Darcy Flow Model
In general, the liquefaction zone may vary based on the aforementioned two liquefaction criteria.Qi and Gao's [51] research revealed that, in Equations ( 5) and ( 6), the parts where p -P b > γ ′ z and j z > γ ′ represent the presence of tensile stress in the non-cohesive seabed.This is obviously nonphysical and leads to the inconsistency in liquefaction zones predicted by the two criteria.Therefore, when liquefaction occurs, even if the wave height continues to increase, the inequality signs in Equations ( 5) and ( 6) should always be equal signs.Thus, we revised Equations ( 5) and ( 6) as follows, respectively: By introducing Equation ( 8) as an additional constraint to Equation (3), Equation (3) can be rewritten as follows: where the pore-fluid velocity v w is determined by where v wx and v wz represent the horizontal and vertical components of v w , respectively.κ is the penalty factor and makes Darcy's law hold in the horizontal direction, while the significant nonlinearity caused by seabed liquefaction appears in the vertical direction.
Because the pressure gradient is related to the hydraulic gradient i (∆p = γ w i), v wz can be written as where i cr equals γ ′ /γ w and represents the critical value that determines the occurrence of liquefaction.i z represents the vertical component of the hydraulic gradient.κ ∞ represents the penalty value and is generally taken as a large value, e.g., 1000.The equations described above are shown in Figure 2 in the v wz − i z space, which appears as the non-Darcy flow model.
where cr i equals / w γ γ ′ and represents the critical value that determines the occurrence of liquefaction.z i represents the vertical component of the hydraulic gradient.κ ∞ represents the penalty value and is generally taken as a large value, e.g., 1000.The equations described above are shown in Figure 2 in the wz z v i space, which ap- pears as the non-Darcy flow model.After discretization, linearization, and other derivation steps (the specific derivation process has been provided by Zhou [63] and is omitted here), the finite element formulation using the non-Darcy flow model is obtained as follows:

Onset of liquefaction
, After discretization, linearization, and other derivation steps (the specific derivation process has been provided by Zhou [63] and is omitted here), the finite element formulation using the non-Darcy flow model is obtained as follows: with where d u and d p denote the discrete unknown vectors corresponding to the displacement (u) and the excessive pore pressure (p), respectively.The shape function matrices of d u and d p are N u and N p .v ' w denotes the partial derivative of v w with respect to i, and v t w,h denotes the current discrete velocity, while the subscripts k and k−1 are used to indicate the nonlinear iterations.
Based on the algebraic representations mentioned above, the iterative process for the solution (d t u , d t p ) can be derived and conveniently applied in conventional finite element codes.

Computational Model
The model without a structure and our finite element code have been validated by Zhou et al. [52,61,63], including cylinder tests under 1-D wave loading and two-dimensional (2-D) wave-seabed interactions.These validations are therefore not repeated here in this work.
This section focuses on the validation of our finite element code applied to the scenario involving a structure, which is taken as a pipeline fully buried in the non-cohesive seabed.In this scenario, a series of flume tests were conducted by Turcotte [67], wherein seabed liquefaction was not induced.The model of these tests is shown in Figure 3.In the tests, the wave period T was 0.9-2.3s and the wave height H was 3.02-14.3cm, while L was 3.53 m, h was 0.533 m, S r was 0.997, n was 0.42, k s was 0.001 m•s −1 , E was 3 MPa, and ν was 0.33.The computational model was then built according to the given conditions. in this work.
This section focuses on the validation of our finite element code applied to the scenario involving a structure, which is taken as a pipeline fully buried in the non-cohesive seabed.In this scenario, a series of flume tests were conducted by Turcotte [67], wherein seabed liquefaction was not induced.The model of these tests is shown in Figure 3.In the tests, the wave period T was 0.9-2.3s and the wave height H was 3.02-14.3cm, while L was 3.53 m, h was 0.533 m, r S was 0.997, n was 0.42, s k was 0.001 m•s −1 , E was 3 MPa, and ν was 0.33.The computational model was then built according to the given condi- tions.

Pore Pressure Response
The distribution of pore pressure around the pipeline is shown in Figure 4.The vertical axis represents the normalized pore pressure 0 / p P .Cheng and Liu [72] also con- ducted numerical studies based on Turcotte's flume experiments.For a better visual comparison, the numerical results of this study were compared with the flume experimental results of Turcotte [67] and the numerical results of Cheng and Liu [72].It can be seen that our numerical results are in basic agreement with the experimental results of Turcotte [67] and the numerical results of Cheng and Liu [72].

Pore Pressure Response
The distribution of pore pressure around the pipeline is shown in Figure 4.The vertical axis represents the normalized pore pressure p/P 0 .Cheng and Liu [72] also conducted numerical studies based on Turcotte's flume experiments.For a better visual comparison, the numerical results of this study were compared with the flume experimental results of Turcotte [67] and the numerical results of Cheng and Liu [72].It can be seen that our numerical results are in basic agreement with the experimental results of Turcotte [67] and the numerical results of Cheng and Liu [72]. Figure 5 shows the distribution of 0 / p P along the depth direction of the centerline using the constant-permeability model (which is hereafter called the CP model) and the liquefaction-associated non-Darcy model (which is hereafter called the ND model) both with and without a pipeline.Due to the absence of liquefaction in the seabed, the results of the CP and ND models are completely consistent.A comparative analysis reveals that the presence of the pipeline significantly affects the transmission of pore pressure in the seabed.Above the pipeline, the pore pressure attenuation is slower in the seabed with a pipeline, but below the pipeline, the pore pressure decreases significantly.This can also Figure 5 shows the distribution of p/P 0 along the depth direction of the centerline using the constant-permeability model (which is hereafter called the CP model) and the liquefaction-associated non-Darcy model (which is hereafter called the ND model) both with and without a pipeline.Due to the absence of liquefaction in the seabed, the results of the CP and ND models are completely consistent.A comparative analysis reveals that the presence of the pipeline significantly affects the transmission of pore pressure in the seabed.Above the pipeline, the pore pressure attenuation is slower in the seabed with a pipeline, but below the pipeline, the pore pressure decreases significantly.This can also be clearly observed in the pore pressure variations along the centerline of the seabed model (e.g., Figure 6, in which the variation in pore pressure versus time step at different depths is displayed).The specific reasons for this phenomenon will be discussed in Section 5.

Computational Model
In this subsection, our liquefaction-associated non-Darcy flow model is utilized to analyze the wave-seabed-pipeline interaction under extreme wave loadings, which induce seabed liquefaction.The computational model adopts a two-dimensional approach.The width of the seabed in the computational model is assumed to be 50 m, with a thickness of 20 m.The computational mesh consisting of quadrilateral elements is shown in Figure 7.The seabed mesh surrounding the pipeline is refined in the meshing process, with a minimum mesh size of 0.01 m.The computational parameters are listed in Table 1.The burial depth and diameter are denoted by d and D, respectively, while d is defined as the distance from the top position of the pipeline to the seabed surface.

Computational Model
In this subsection, our liquefaction-associated non-Darcy flow model is utilized to analyze the wave-seabed-pipeline interaction under extreme wave loadings, which induce seabed liquefaction.The computational model adopts a two-dimensional approach.The width of the seabed in the computational model is assumed to be 50 m, with a thickness of 20 m.The computational mesh consisting of quadrilateral elements is shown in Figure 7.The seabed mesh surrounding the pipeline is refined in the meshing process, with a minimum mesh size of 0.01 m.The computational parameters are listed in Table 1.The burial depth and diameter are denoted by d and D, respectively, while d is defined as the distance from the top position of the pipeline to the seabed surface.

Computational Model
In this subsection, our liquefaction-associated non-Darcy flow model is utilized to analyze the wave-seabed-pipeline interaction under extreme wave loadings, which induce seabed liquefaction.The computational model adopts a two-dimensional approach.The width of the seabed in the computational model is assumed to be 50 m, with a thickness of 20 m.The computational mesh consisting of quadrilateral elements is shown in Figure 7.The seabed mesh surrounding the pipeline is refined in the meshing process, with a minimum mesh size of 0.01 m.The computational parameters are listed in Table 1.The burial depth and diameter are denoted by d and D, respectively, while d is defined as the distance from the top position of the pipeline to the seabed surface.

Pore Pressure Response
To investigate the influence of a fully buried pipeline on the pore pressure r in a non-cohesive seabed, three values of burial depth d are considered, as pres Table 2, wherein the pipeline diameters are fixed at the same value of burial de d/D = 1.Figures 8-13 give the temporal distributions of pore pressure at the top of the pipeline and the vertical distributions of pore pressure above the top of the under three different schemes.In these figures, b P represents the wave pressure to the seabed surface.The wave pressure is instantaneously transmitted from the sea surface to the surface.However, due to the unsaturated nature of the seabed and the compress pore water flow, the wave pressure cannot be promptly transmitted within the resulting in phase lag over time (e.g., Figures 8,10 and 12) and amplitude atte along the depth direction (e.g., Figures 9,11 and 13).The amplitude attenuation an lag are inherent characteristics of the wave-induced instantaneous pore pres sponse.

Pore Pressure Response
To investigate the influence of a fully buried pipeline on the pore pressure response in a non-cohesive seabed, three values of burial depth d are considered, as presented in Table 2, wherein the pipeline diameters are fixed at the same value of burial depth, i.e., d/D = 1.Figures 8-13 give the temporal distributions of pore pressure at the top position of the pipeline and the vertical distributions of pore pressure above the top of the pipeline under three different schemes.In these figures, P b represents the wave pressure applied to the seabed surface.

Pore Pressure Response
To investigate the influence of a fully buried pipeline on the pore pressure response in a non-cohesive seabed, three values of burial depth d are considered, as presented in Table 2, wherein the pipeline diameters are fixed at the same value of burial depth, i.e., d/D = 1.Figures 8-13 give the temporal distributions of pore pressure at the top position of the pipeline and the vertical distributions of pore pressure above the top of the pipeline under three different schemes.In these figures, b P represents the wave pressure applied to the seabed surface.The wave pressure is instantaneously transmitted from the sea surface to the seabed surface.However, due to the unsaturated nature of the seabed and the compressibility of pore water flow, the wave pressure cannot be promptly transmitted within the seabed, resulting in phase lag over time (e.g., Figures 8,10 and 12) and amplitude attenuation along the depth direction (e.g., Figures 9,11 and 13).The amplitude attenuation and phase lag are inherent characteristics of the wave-induced instantaneous pore pressure response.According to Figures 8, 10 and 12, the results obtained from the CP model and ND models are close to each other.The variation in pore pressure at the top position of the pipeline follows the waveform of the wave motion.Under Scheme I, with a pipeline burial depth of 0.042 m, the wave pressure amplitude 0 p at the seabed surface is 1500 Pa, while the pore pressure amplitude at the top position of the pipeline is 1300 Pa, resulting in an amplitude attenuation of 13.3%.Compared to the amplitude of the wave pressure b P ap- plied to the seabed surface, the pore pressure amplitude at the top position of the pipeline is delayed by about two computational steps, indicating a phase lag time of about 0.05T.Under Scheme II, with a pipeline burial depth of 0.084 m, 0 p at the seabed surface is 1500 Pa, while the pore pressure amplitude at the top position of the pipeline is 1200 Pa, resulting in an amplitude attenuation of 20.0%.Compared to the amplitude of b P applied to the seabed surface, the pore pressure amplitude at the top of the pipeline is delayed by about three computational steps, indicating a phase lag time of about 0.075T.Under Scheme III, with a pipeline burial depth of 0.126 m, b P at the seabed surface is 1500 Pa, while 0 p at the top position of the pipeline is 1100 Pa, resulting in an amplitude attenu- ation of 26.7%.The pore pressure amplitude at the top position of the pipeline appears at the same time as in Scheme II and is delayed by about three computational steps relative to the amplitude of b P applied to the seabed surface, indicating a phase lag time of about 0.075T.According to the above analysis, it can be concluded that as the burial depth of the pipeline increases, the degree of amplitude attenuation and the phase lag also increase.According to Figures 9, 11 and 13, the pore pressure attenuation at a depth z = 0.04 m can be analyzed.The attenuation under Schemes I, II, and III is 13%, 25%, and 26.9%, respectively.With an increase in the structural burial depth, the degree of pore pressure attenuation also increases.In other words, a fully buried pipeline can mitigate the amplitude attenuation of pore pressure transmission within the seabed.
It can be seen that pipelines of different sizes and at different burial depths have different effects on the transfer of pore pressure in the seabed.The weakening of the seabed by the fully buried pipeline in the seabed is mainly due to the significant influence of its diameter and burial depth on the pore pressure transfer in the seabed, which is mainly manifested in the pore pressure amplitude attenuation and phase lag.The amplitude attenuation and phase lag are also the main reasons for the instantaneous liquefaction of the seabed.
Figure 14 shows the liquefaction zones from the CP and ND models in the above three schemes.It can be seen that the liquefaction depth above the pipeline is nearly close to zero due to the liquefaction shielding effect, which is discussed in Section 5.As a result, the pore pressure distributions above the pipeline from the CP and ND models are close to each other, as shown in Figures 8-13.To better illustrate the difference between the CP The wave pressure is instantaneously transmitted from the sea surface to the seabed surface.However, due to the unsaturated nature of the seabed and the compressibility of pore water flow, the wave pressure cannot be promptly transmitted within the seabed, resulting in phase lag over time (e.g., Figures 8,10 and 12) and amplitude attenuation along the depth direction (e.g., Figures 9,11 and 13).The amplitude attenuation and phase lag are inherent characteristics of the wave-induced instantaneous pore pressure response.
According to Figures 8, 10 and 12, the results obtained from the CP model and ND models are close to each other.The variation in pore pressure at the top position of the pipeline follows the waveform of the wave motion.Under Scheme I, with a pipeline burial depth of 0.042 m, the wave pressure amplitude p 0 at the seabed surface is 1500 Pa, while the pore pressure amplitude at the top position of the pipeline is 1300 Pa, resulting in an amplitude attenuation of 13.3%.Compared to the amplitude of the wave pressure P b applied to the seabed surface, the pore pressure amplitude at the top position of the pipeline is delayed by about two computational steps, indicating a phase lag time of about 0.05T.Under Scheme II, with a pipeline burial depth of 0.084 m, p 0 at the seabed surface is 1500 Pa, while the pore pressure amplitude at the top position of the pipeline is 1200 Pa, resulting in an amplitude attenuation of 20.0%.Compared to the amplitude of P b applied to the seabed surface, the pore pressure amplitude at the top of the pipeline is delayed by about three computational steps, indicating a phase lag time of about 0.075T.Under Scheme III, with a pipeline burial depth of 0.126 m, P b at the seabed surface is 1500 Pa, while p 0 at the top position of the pipeline is 1100 Pa, resulting in an amplitude attenuation of 26.7%.The pore pressure amplitude at the top position of the pipeline appears at the same time as in Scheme II and is delayed by about three computational steps relative to the amplitude of P b applied to the seabed surface, indicating a phase lag time of about 0.075T.According to the above analysis, it can be concluded that as the burial depth of the pipeline increases, the degree of amplitude attenuation and the phase lag also increase.
According to Figures 9, 11 and 13, the pore pressure attenuation at a depth z = 0.04 m can be analyzed.The attenuation under Schemes I, II, and III is 13%, 25%, and 26.9%, respectively.With an increase in the structural burial depth, the degree of pore pressure attenuation also increases.In other words, a fully buried pipeline can mitigate the amplitude attenuation of pore pressure transmission within the seabed.
It can be seen that pipelines of different sizes and at different burial depths have different effects on the transfer of pore pressure in the seabed.The weakening of the seabed by the fully buried pipeline in the seabed is mainly due to the significant influence of its diameter and burial depth on the pore pressure transfer in the seabed, which is mainly manifested in the pore pressure amplitude attenuation and phase lag.The amplitude attenuation and phase lag are also the main reasons for the instantaneous liquefaction of the seabed.
Figure 14 shows the liquefaction zones from the CP and ND models in the above three schemes.It can be seen that the liquefaction depth above the pipeline is nearly close to zero due to the liquefaction shielding effect, which is discussed in Section 5.As a result, the pore pressure distributions above the pipeline from the CP and ND models are close to each other, as shown in Figures 8-13.To better illustrate the difference between the CP and ND models, another computational scheme (d = D = 0.168 m) is considered, with the numerical results shown in Figure 15. Figure 15a,b are the liquefaction zones from the CP and ND models, respectively.It can be seen that the liquefaction zone from the ND model is apparently smaller than that from the CP model.Figure 15c is the vertical effective stress above the pipeline centerline, showing that the ND model effectively mitigates the presence of tensile stresses caused by the CP model.

Parametric Study on the Liquefaction Depth
According to the research conducted by Zhou et al. [52,61,63], when there is no liquefied zone in the non-cohesive seabed, the results obtained by the ND model are consistent with those by the model.Once liquefaction occurs in the non-cohesive seabed, the CP model results in unreasonable tensile stresses within the liquefied zone, whereas the ND model does not produce such results.Therefore, Zhou et al. [52,61,63] suggest using the CP model for calculations in practical engineering scenarios first and employing the ND model to eliminate tensile stresses if liquefaction is observed.Through the parametric study by Zhou et al. [52,61,63], it was found that the instantaneous liquefaction depth from the ND model is approximately 0.73 times that from the CP model.This relationship is obtained in the scenario of a pure seabed without structures.Here, in this subsection, we further describe a parametric study conducted for the scenario involving a pipeline fully buried in the seabed.
During the parametric study, a benchmark test was set with the computational parameters listed in Table 1, and we considered seven different parameters (Young's modulus , permeability coefficient s k , saturation degree r S , soil porosity n, wave period T, water depth h, and wave height H).Building upon the model presented in Section 3.1, d and D of the pipeline are taken as 0.083 m and 0.168 m, respectively.In Figure 16, CP represents the results of the CP model, while ND represents the results of the ND model, and the solid line represents the model with a pipeline, while the dashed line represents the model without a pipeline.Figure 16a shows that the liquefaction depth increases rapidly when the elastic modulus is between 1 and 3 MPa.As the elastic modulus exceeds 3 MPa, the rate of increase in liquefaction depth slows down and gradually approaches a constant value.In Figure 16b, it can be seen that the liquefaction depth decreases with the increase in the permeability coefficient.Figure 16c demonstrates the sensitivity of the liquefaction depth versus soil saturation.The presence of air in unsaturated soils is a key factor causing the amplitude attenuation and phase lag of pore pressure propagation in the seabed, thus directly affecting the liquefaction depth.Figure 16d shows that the liquefaction depth has a nearly linear relationship with the void ratio.In Figure 16e, it is evident that the liquefaction depth decreases linearly with the increase in water depth.Figure 16f reveals that the liquefaction depth increases and gradually approaches a constant value as the wave period increases.In Figure 16g, the change in liquefaction depth exhibits a roughly linear relationship with the wave height.The liquefaction depth increases with the increase in wave height and suddenly increases when the wave height reaches 0.8 m.This is because, during the process of increasing the wave

Parametric Study on the Liquefaction Depth
According to the research conducted by Zhou et al. [52,61,63], when there is no liquefied zone in the non-cohesive seabed, the results obtained by the ND model are consistent with those by the CP model.Once liquefaction occurs in the non-cohesive seabed, the CP model results in unreasonable tensile stresses within the liquefied zone, whereas the ND model does not produce such results.Therefore, Zhou et al. [52,61,63] suggest using the CP model for calculations in practical engineering scenarios first and employing the ND model to eliminate tensile stresses if liquefaction is observed.Through the parametric study by Zhou et al. [52,61,63], it was found that the instantaneous liquefaction depth from the ND model is approximately 0.73 times that from the CP model.This relationship is obtained in the scenario of a pure seabed without structures.Here, in this subsection, we further describe a parametric study conducted for the scenario involving a pipeline fully buried in the seabed.
During the parametric study, a benchmark test was set with the computational parameters listed in Table 1, and we considered seven different parameters (Young's modulus E, permeability coefficient k s , saturation degree S r , soil porosity n, wave period T, water depth h, and wave height H).Building upon the model presented in Section 3.1, d and D of the pipeline are taken as 0.083 m and 0.168 m, respectively.
In Figure 16, CP represents the results of the CP model, while ND represents the results of the ND model, and the solid line represents the model with a pipeline, while the dashed line represents the model without a pipeline.Figure 16a shows that the liquefaction depth increases rapidly when the elastic modulus is between 1 and 3 MPa.As the elastic modulus exceeds 3 MPa, the rate of increase in liquefaction depth slows down and gradually approaches a constant value.In Figure 16b, it can be seen that the liquefaction depth decreases with the increase in the permeability coefficient.Figure 16c demonstrates the sensitivity of the liquefaction depth versus soil saturation.The presence of air in unsaturated soils is a key factor causing the amplitude attenuation and phase lag of pore pressure propagation in the seabed, thus directly affecting the liquefaction depth.Figure 16d shows that the liquefaction depth has a nearly linear relationship with the void ratio.In Figure 16e, it is evident that the liquefaction depth decreases linearly with the increase in water depth.Figure 16f reveals that the liquefaction depth increases and gradually approaches a constant value as the wave period increases.In Figure 16g, the change in liquefaction depth exhibits a roughly linear relationship with the wave height.The liquefaction depth increases with the increase in wave height and suddenly increases when the wave height reaches 0.8 m.This is because, during the process of increasing the wave height from 0.6 m to 0.8 m, a small liquefied zone also appears at the bottom of the pipeline, as shown in Figure 17.It can be therefore concluded that E, n, T, and H positively influence the liquefaction depth, while k s , S r , and h have a negative impact on it.Additionally, the liquefaction depth in the scenario with a pipeline is generally smaller than that in the scenario without a pipeline.height from 0.6 m to 0.8 m, a small liquefied zone also appears at the bottom of the pipeline, as shown in Figure 17.It can be therefore concluded that E, n, T, and H positively influence the liquefaction depth, while s k , r S , and h have a negative impact on it.Addi- tionally, the liquefaction depth in the scenario with a pipeline generally smaller than that in the scenario without a pipeline.) in the scenario involving a pipeline.The ratio between the instantaneo liquefaction depth from the ND model and the CP model is approximately 0.73, which close to that in the scenario of a pure seabed without structures.

Study on Non-Dimensional Parameters
To establish a quantifiable relationship between wave loads and structural respon with good universality, it is necessary to develop an expression based on non-dimensio parameters.To achieve this, adopting schemes with diameters ranging from 0.042 to 0.1 m and burial depths of 1-5 times the diameters (the specific values of burial depths a diameters are given in Table 3, while r S is 0.96, ν is 0.3, and the other parameters the same as those in Table 1), we conducted a large number of numerical calculations a attempted various combinations of horizontal and vertical parameters.Eventually, found that when the horizontal coordinate is / L D (wavelength/diameter) and the v tical coordinate is ND / z d (liquefaction depth/burial depth), there is a significant lin relationship between them, as shown in Figure 18.The correlation coefficient (   According to Figure 16h (where z CP represents the liquefaction depth from the CP model and z ND represents that from the ND model, which are all from Figure 16a-g), it can be observed that there is a significant linear correlation between z CP and z ND (with R 2 = 0.9787) in the scenario involving a pipeline.The ratio between the instantaneous liquefaction depth from the ND model and the CP model is approximately 0.73, which is close to that in the scenario of a pure seabed without structures.

Study on Non-Dimensional Parameters
To establish a quantifiable relationship between wave loads and structural responses with good universality, it is necessary to develop an expression based on non-dimensional parameters.To achieve this, adopting schemes with diameters ranging from 0.042 to 0.168 m and burial depths of 1-5 times the diameters (the specific values of burial depths and diameters are given in Table 3, while S r is 0.96, ν is 0.3, and the other parameters are the same as those in Table 1), we conducted a large number of numerical calculations and attempted various combinations of horizontal and vertical parameters.Eventually, we found that when the horizontal coordinate is L/D (wavelength/diameter) and the vertical coordinate is z ND /d (liquefaction depth/burial depth), there is a significant linear relationship between them, as shown in Figure 18.The correlation coefficient (R 2 ) was as high as 0.9223, indicating a good fitting effect.In Section 4.3, it is shown that the presence of a pipeline in the seabed leads to a smaller liquefaction depth compared to the scenario without a pipeline.When a pipeline exists in the seabed, the liquefied zones above the pipeline exhibit a spatial distribution that "bypasses" the pipeline, resulting in a "liquefaction shielding phenomenon".Even if the liquefaction depth on both sides of the pipeline is below the top position of the pipeline, there will be a layer of non-liquefied soil above the pipeline, e.g., Figure 19a.As shown in Figure 19b-d, as the burial depth of the pipeline increases, keeping the diameter constant, the liquefaction depth also increases.However, when the burial depth is large enough, its further increase will not induce a further increase in the liquefaction depth.The reason is that instantaneous liquefaction generally occurs in a thin layer near the seabed surface.When the pipeline is far away from the seabed surface, the presence of the pipeline will no longer influence the liquefied zone.In contrast, when the buried depth is not large enough, the presence of the pipeline can have a significant influence on the liquefied zone.Below, a tentative explanation is for this phenomenon.
During the numerical simulation, the presence of the pipeline acts as an impermeable boundary in the seabed, causing a diversion of the flow field in the surrounding area, as shown in Figure 20.The flow does not reach the soil above the pipeline immediately, resulting in an impact on the phase lag and amplitude attenuation of the wave-induced seabed pore pressure response, and this phenomenon is referred to as the "liquefaction shielding phenomenon".In other words, when the structure is shallow, the structure will repel the liquefaction zone above it, and the area above the structure where liquefaction does not occur is referred to as the "shielding ring".It is therefore straightforward to imagine that when the pipeline is shallow enough, the liquefied zone above the pipeline can completely disappear.This prediction is studied in the following subsection.

The Shielding Effect of the Pipeline on the Liquefied Zone
In Section 4.3, it is shown that the presence of a pipeline in the seabed leads to a smaller liquefaction depth compared to the scenario without a pipeline.When a pipeline exists in the seabed, the liquefied zones above the pipeline exhibit a spatial distribution that "bypasses" the pipeline, resulting in a "liquefaction shielding phenomenon".Even if the liquefaction depth on both sides of the pipeline is below the top position of the pipeline, there will be a layer of non-liquefied soil above the pipeline, e.g., Figure 19a.As shown in Figure 19b-d, as the burial depth of the pipeline increases, keeping the diameter constant, the liquefaction depth also increases.However, when the burial depth is large enough, its further increase will not induce a further increase in the liquefaction depth.The reason is that instantaneous liquefaction generally occurs in a thin layer near the seabed surface.When the pipeline is far away from the seabed surface, the presence of the pipeline will no longer influence the liquefied zone.In contrast, when the buried depth is not large enough, the presence of the pipeline can have a significant influence on the liquefied zone.Below, a tentative explanation is provided for this phenomenon.During the numerical simulation, the presence of the pipeline acts as an impermeable boundary in the seabed, causing a diversion of the flow field in the surrounding area, as shown in Figure 20.The flow does not reach the soil above the pipeline immediately, resulting in an impact on the phase lag and amplitude attenuation of the wave-induced seabed pore pressure response, and this phenomenon is referred to as the "liquefaction shielding phenomenon".In other words, when the structure is shallow, the structure will repel the liquefaction zone above it, and the area above the structure where liquefaction does not occur is referred to as the "shielding ring".It is therefore straightforward to imagine that when the pipeline is shallow enough, the liquefied zone above the pipeline can completely disappear.This prediction is studied in the following subsection.

Onset Conditions for Liquefaction Occurrence above the Pipeline
According to Section 5.1, the shielding phenomenon of the pipeline on the liquefied zone weakens as the burial depth of the pipeline increases.This subsection will present a computational analysis of various combinations of burial depth d and diameter D to study the onset conditions for liquefaction occurrence above the pipeline.During the numerical simulation, the presence of the pipeline acts as an impermeable boundary in the seabed, causing a diversion of the flow field in the surrounding area, as shown in Figure 20.The flow does not reach the soil above the pipeline immediately, resulting in an impact on the phase lag and amplitude attenuation of the wave-induced seabed pore pressure response, and this phenomenon is referred to as the "liquefaction shielding phenomenon".In other words, when the structure is shallow, the structure will repel the liquefaction zone above it, and the area above the structure where liquefaction does not occur is referred to as the "shielding ring".It is therefore straightforward to imagine that when the pipeline is shallow enough, the liquefied zone above the pipeline can completely disappear.This prediction is studied in the following subsection.

Onset Conditions for Liquefaction Occurrence above the Pipeline
According to Section 5.1, the shielding phenomenon of the pipeline on the liquefied zone weakens as the burial depth of the pipeline increases.This subsection will present a computational analysis of various combinations of burial depth d and diameter D to study the onset conditions for liquefaction occurrence above the pipeline.

Onset Conditions for Liquefaction Occurrence above the Pipeline
According to Section 5.1, the shielding phenomenon of the pipeline on the liquefied zone weakens as the burial depth of the pipeline increases.This subsection will present a computational analysis of various combinations of burial depth d and diameter D to study the onset conditions for liquefaction occurrence above the pipeline.
Based on the numerical results in Section 4.4, we continued to calculate the liquefaction depths in the scenario involving a pipeline considering various schemes with diameters of 0.084 m, 0.126 m, and 0.168 m and burial depths ranging from 0.084 m to 0.4 m (please refer to Table 4 for specific parameters and results).The results in Section 4.4 and Table 4 are jointly plotted in Figure 21.When d/D < 2, although several schemes have liquefaction depths as high as 0.06 m, nearly half of the schemes did not show liquefied zones; i.e., the pipeline is near the seabed surface, and therefore, the aforementioned shielding effect is apparent.However, when d/D ≥ 2, liquefaction zones occur in all the numerical results, with a minimum liquefaction depth of 0.015 m; i.e., the pipeline is far away from the seabed surface and therefore has a slight influence on the liquefied zone.

Conclusions
This paper extends the application of our liquefaction-associated non-Darcy flow model from wave-seabed interactions to wave-seabed-pipeline interactions in order to investigate the influence of a fully buried pipeline on the pore pressure response and the instantaneous liquefaction of the non-cohesive seabed.Based on the numerical examples presented, the following conclusions can be drawn.
(1) Numerical simulations were conducted based on Turcotte's flume experiment, validating the reliability of our in-house code in modeling wave-seabed-structure interactions.(2) A numerical study was conducted by varying the burial depth and diameter of the pipeline, revealing that the existence of a pipeline weakens the degree of amplitude attenuation and the phase lag.Therefore, when the pipeline is shallow, the liquefied zone of the seabed with a pipeline is smaller than that in the pure seabed, which is called the "liquefaction shielding effect" in this work.(3) Under some conditions, liquefaction can even completely disappear above the pipeline, while horizontally distant areas still have liquefied zones.The onset conditions for liquefaction occurrence above the pipeline are then discussed.(4) As the burial depth of the pipeline increases, the liquefaction shielding effect weakens, resulting in an increase in the liquefaction depth above the pipeline.Once the pipeline is sufficiently far from the seabed surface, it no longer influences the liquefied zone.( 5) Based on the parametric study, it was observed that the liquefaction depth predicted by the non-Darcy model is approximately 0.73 times the value estimated by the con- It seems that a critical value of d/D can be quantitatively determined to distinguish whether liquefaction occurs above the pipeline.From Figure 21, this value can be taken as 2; i.e., when d/D ≥ 2, liquefaction always occurs above the pipeline.It should be noted that this critical value can change if other computational parameters are applied.It is also notable that even if there is no liquefied zone above the pipeline, in the horizontal direction, if a position is far from the pipeline, the shielding effect of the pipeline can also disappear, which means that liquefied zones can also occur in a shallow layer of the seabed within these areas, e.g., Figure 19b.

Conclusions
This paper extends the application of our liquefaction-associated non-Darcy flow model from wave-seabed interactions to wave-seabed-pipeline interactions in order to investigate the influence of a fully buried pipeline on the pore pressure response and the instantaneous liquefaction of the non-cohesive seabed.Based on the numerical examples presented, the following conclusions can be drawn.
(1) Numerical simulations were conducted based on Turcotte's flume experiment, validating the reliability of our in-house code in modeling wave-seabed-structure interactions.(2) A numerical study was conducted by varying the burial depth and diameter of the pipeline, revealing that the existence of a pipeline weakens the degree of amplitude attenuation and the phase lag.Therefore, when the pipeline is shallow, the liquefied zone of the seabed with a pipeline is smaller than that in the pure seabed, which is called the "liquefaction shielding effect" in this work.(3) Under some conditions, liquefaction can even completely disappear above the pipeline, while horizontally distant areas still have liquefied zones.The onset conditions for liquefaction occurrence above the pipeline are then discussed.(4) As the burial depth of the pipeline increases, the liquefaction shielding effect weakens, resulting in an increase in the liquefaction depth above the pipeline.Once the pipeline is sufficiently far from the seabed surface, it no longer influences the liquefied zone.(5) Based on the parametric study, it was observed that the liquefaction depth predicted by the non-Darcy model is approximately 0.73 times the value estimated by the conventional Darcy model, regardless of whether or not a pipeline is involved.(6) The quantitative relationship between wave loadings and structural sizes is studied.
A highly linear relationship between two sets of non-dimensional parameters, i.e., "liquefaction depth/burial depth" and "wavelength/diameter", is discovered.

Figure 1 . 2 ×Ι
Figure 1.Simplified wave-seabed-structure model with instantaneous liquefaction.The first two equations are derived from the theory of poro-elasticity theory [70].Equation (3a) denotes the equilibrium of the solid-fluid mixture, in which σ ′ represents the effective stress, 2 2 × Ι denotes the second-order unit tensor, and b denotes the body force per unit volume.Equation (3b) denotes the conservation of mass, in which v ε rep-

JFigure 4 .
Figure 4.A comparison of the wave-induced pore pressure amplitude along the periphery of the fully buried pipeline between laboratory data by Turcotte [67] and numerical simulations by Cheng and Liu [72]: (a) T = 0.9 s, H = 5.24 cm; (b) T = 1.75 s, H = 14.3 cm; (c) T = 2.3 s, H = 3.02 cm.

Figure 4 .
Figure 4.A comparison of the wave-induced pore pressure amplitude along the periphery of the fully buried pipeline between laboratory data by Turcotte [67] and numerical simulations by Cheng and Liu [72]: (a) T = 0.9 s, H = 5.24 cm; (b) T = 1.75 s, H = 14.3 cm; (c) T = 2.3 s, H = 3.02 cm.

Figure 6 .
Figure 6.Pore pressure versus time step along the x-directional centerline at different depths.

Figure 8 .
Figure 8. Pore pressure versus time step at the top position of the pipeline under Scheme I.

Figure 8 .Figure 9 .Figure 10 .Figure 11 .
Figure 8. Pore pressure versus time step at the top position of the pipeline under Scheme I.

Figure 13 .
Figure 13.Pore pressure versus depth above the top position of the pipeline under Scheme III.

JFigure 14 .Figure 14 .Figure 14 .Figure 15 .
Figure 14.Liquefaction zones from CP and ND models in Schemes I~III: (a) from CP model in Scheme I; (b) from ND model in Scheme I; (c) from CP model in Scheme II; (d) from ND model in Scheme II; (e) from CP model in Scheme III; (f) from ND model in Scheme III.

Figure 15 .
Figure 15.Liquefaction zones obtained by CP and ND models in the scheme of d = D = 0.168 m: CP model; (b) ND model; (c) vertical distribution of γ ′ z − (p − P b ).

2 R
) was high as 0.9223, indicating a good fitting effect.

Figure 17 .
Figure 17.Liquefaction zones: (a) H = 0.6 m without a pipeline; (b) H = 0.8 m without a pipeline; (c) H = 0.6 m with a pipeline; (d) H = 0.8 m with a pipeline.

Figure 18 .
Figure 18.ND / z d versus L/D in ND model.

Figure 20 .
Figure 20.A schematic diagram of flow around the pipeline.

Figure 20 .
Figure 20.A schematic diagram of flow around the pipeline.

Figure 20 .
Figure 20.A schematic diagram of flow around the pipeline.

Table 2 .
Computational schemes with various buried depths and diameters.

Table 1 .
Computational parameters for the model with seabed liquefaction.

Table 2 .
Computational schemes with various buried depths and diameters.

Table 1 .
Computational parameters for the model with seabed liquefaction.

Table 2 .
Computational schemes with various buried depths and diameters.

Table 3 .
The scheme of the burial depth and diameter.

Table 3 .
The scheme of the burial depth and diameter.

Table 4 .
Liquefaction depth under different working schemes.