Assessment of Wave–Current-Induced Liquefaction under Twin Pipelines Using the Coupling Model

Abstract

Although twin pipelines in series have been used to transport hydrocarbons in engineering practice, most previous studies focused on the dynamic response of the seabed around a single pipeline. A two-way coupling model of fluid–structure–seabed interaction (FSSI) is proposed for the study of the soil response and liquefaction caused by waves and currents around twin pipelines. The present model integrates the flow model and the seabed model by introducing a boundary condition of velocity continuity in addition to the continuity of pressures at the seabed surface. Then, the inconsistency between the physical process and numerical simulation can be overcome in the one-way coupling model. Through a series of numerical simulations, the influence of different flow characteristics, soil properties, and pipeline configurations on the seabed response under the two-way coupling process were explored, and compared with the results of the single pipeline. The numerical results indicate that the twin pipeline configuration significantly alters the relevant responses compared to the single pipeline configuration, including the after-consolidation state, amplitude of velocity at the seabed surface, and distribution of pore pressure in the seabed. The parametric studies show that the amplitudes of the wave and current have significant impacts on the distribution of pore pressure in the seabed. The pore pressure in the seabed increases with the increase of forward wave current, while the results of reverse wave current are the opposite. In addition, the liquefaction range around the pipeline increases with the increase of $H_w$ and $T_w$, and increases with the decrease of $S_r$ and $k_s$. At the same time, the gaps (G) and the ratio of pipe radius ($R_1/R_2$) between the twin pipelines also significantly affect the seabed response and liquefaction distribution around the pipeline.

1. Introduction

Submarine pipelines have been widely used for gas and oil transport in marine environments. The presence of a submarine pipeline cannot only change the flow morphology around it but also aggravate the instability of the bottom of the pipeline and lead to its destruction [1]. Since the marine environment is complex, there are various factors affecting the stability of submarine pipelines. Numerous investigations have been carried out in the past. Among these, some researchers evaluated the pipeline stability, considering the seabed response and the corresponding liquefaction around pipelines under the actions of waves or combined wave–current loadings, adopting different numerical models [2,3,4]. Numerous researchers have investigated the hydrodynamics and scour processes around offshore structures in marine environments [5,6,7,8], because the hydrodynamic forces and scour holes caused by sediment transport due to complex hydrodynamics in the vicinity of the pipeline could be other vital threats to its stability.

The present work focuses on the process of wave–current–seabed–pipeline interactions and the subsequent soil response and liquefaction phenomenon near the pipeline. According to field measurements and laboratory experiments [9], the soil response caused by waves can be divided into two mechanisms. The first type is called the transient mechanism, which shows a period-averaged oscillatory soil response [10,11]. The second mechanism, known as the residual mechanism, is caused by the accumulation of pore pressure caused by the shrinkage of saturated seabed soil under cyclic loading [12,13]. In this study, we mainly investigate the transient mechanism.

There are two new contributions of the present study:

• The present model considers twin pipelines operating in parallel, whose application has significantly increased with the expansion projects of oil and gas pipelines in later stages;

• A two-way coupling model was developed, considering the interactions between the flow and solid regions in order to conduct a more realistic simulation.

Although two parallel submarine pipelines have been adopted in engineering practices to transport natural gas, most relevant research studies were limited to the study of a single pipeline. With the development of the underwater oil and gas transportation industry, it was necessary to construct a new pipeline route along the existing one or directly lay multiple pipelines in parallel on the seabed during the initial construction phase to meet the increasing demand for oil and gas transportation. When multiple pipelines exist simultaneously, the nearby hydrodynamics as well as the response of the seabed around the pipelines are unpredictable [14]. For example, the Nord Stream project is a giant gas pipeline project in Europe that crosses the Baltic Sea to transport natural gas to Western Europe, with a total length of 1224 km of subsea pipelines. The two phases of the Nord Stream project are parallel along the same route, providing 1100 billion cubic meters of natural gas to Europe annually. Such changes in the geometries of submarine pipelines require extensive studies to ensure their safety.

To date, only a few studies considered the cases of twin pipelines in tandem despite the fact that they are quite common in the offshore oil and gas industry. Among these, Zhao and Cheng [15] studied the local scour under a piggyback pipeline in a steady flow, consisting of twin pipelines with different diameters. Jo et al. [16] simulated dual pipelines in trenches, and the results indicate that both the depth and slope of the trench can affect the stability of the pipeline. Zhang et al. [17] studied the effects of burial depth and center spacing on the interaction between double-row pipelines under wave actions through a numerical simulation. The results indicate that both the burial depth of the pipeline and the relative position between the two pipelines can affect the distribution of pore pressure around the pipeline. Zhai et al. [18] examined the influences of different burial depths and spacings between two series pipelines on the response of the seabed around the pipeline by conducting experiments. Later, Chen et al. [19,20] numerically investigated the transient liquefaction and residual liquefaction around twin pipelines with identical sizes, respectively. However, none of the above works considered twin pipelines with different pipeline diameters, which could occur in the expansion projects of existing pipelines.

Regarding the second novelty of the present work, the two-way coupling algorithm adopted is more consistent with the actual physical processes. As shown in Figure 1a, the conventional one-way coupling model obtains the wave pressures from the flow model as the only external load of the seabed model. In this approach, the bottom boundary of the flow model is considered an impermeable boundary. However, the seabed is treated as a porous medium in the soil response model. This is an obvious contraction between physical phenomena and theoretical approaches. Furthermore, the influence of the seabed characteristic on the flow field cannot be reproduced in the existing one-way coupling approach. To overcome the shortcoming of the one-way coupling approach, in addition to the continuity of pressures at the seabed surface, an additional boundary condition of velocity continuity can be added at the bottom of the flow model (Figure 1b). This will reflect the impact of seabed properties on the flow model, which is named the two-way coupling model. To date, only a few researchers have used two-way coupling algorithms to study wave–seabed interactions. For example, Zhai and Jeng [21] proposed a new two-way coupling model; an impermeable seabed in the existing flow model is released, and the continuity of fluid velocity is introduced by incorporating seepage velocity and velocity resulting from soil displacements obtained from the seabed model.

This study will adopt a two-way coupling algorithm to study the soil response around twin pipelines with different diameters in tandem. Firstly, the present model was validated with existing experimental data. With the numerical model, the wave characteristics around the pipeline, seabed response, and liquefaction range of twin pipelines are compared with the single pipeline by using the two-way coupling algorithm. Finally, a series of parameter analyses are conducted to explore their influence on the liquefaction range around the twin pipelines.

2. Numerical Model

The sketch of wave–seabed–twin pipelines in tandem is shown in Figure 2. As shown in the figure, the 2D coordinate system is established based on the x-axis of the seabed surface in the direction of wave propagation, the z-axis perpendicular to the interface of the fluid and seabed, and through the center of the upstream pipeline. As suggested in the literature [22], to generate a stable current simulation, a uniform flow velocity is provided on the left-side boundary of the fluid sub-model. After the stable operation of currents, ocean waves are generated in the fluid sub-model, which can achieve the impact of the combined action of waves and currents on the seabed. Since the wave generation in the present model is based on the velocity boundary at the inlet, it does not require an additional equation for the generation of currents.

2.1. Flow Model

In the present research, the flow model (olaFlow) was developed by Higuera et al. [23] with modifications to the boundary conditions in the fluid bottom. In the flow model, the VARANS (volume-averaged Reynolds-averaged Navier–Stokes) equations are employed as the governing equations. The key equations are outlined here:

(1)$\nabla ·\mathit{U}=0,$

(2)$\frac{\partial \rho \mathit{U}}{\partial t}+\rho (\nabla ·\mathit{U})\mathit{U}+\nabla {\mathit{p}}_{\mathit{s}}=\nabla ·\tau +\rho \mathit{g},$

(3)$\frac{\partial \alpha }{\partial t}+\nabla ·\left(\alpha \mathit{U}\right)+\nabla ·\left[\alpha (1-\alpha ){\mathit{u}}_{\mathit{r}}\right]=0,$

where $\alpha$ represents the volume fraction of the fluid:

(4)$\begin{array}{c}\hfill \alpha =\left\{\begin{array}{cc}0,\hfill & \mathrm{air}\hfill \\ 1,\hfill & \mathrm{water}\hfill \\ 0<\alpha <1,\hfill & \mathrm{free} \mathrm{surface}\hfill \end{array}\right.\end{array}$

The present model, different from previous flow sub-models [4,24], integrates the flow model and the seabed model by adopting a two-way coupling algorithm, where the boundary conditions of velocity continuity and pressure continuity are adopted at the fluid and seabed interface. This overcomes the contradiction between the physical model and the actual situation. This boundary condition is obtained from the combination of seepage velocity and soil displacement-induced velocity:

(5)$\mathit{U}=-\frac{k_s}{{\gamma }_w}\nabla p_s+\frac{\partial {\mathbf{u}}_{\mathbf{s}}}{\partial t}atz=0$

In the conventional flow model, an impermeable boundary between the seabed and flow region where the normal component of the relative velocity is zero is adopted [23]. However, the present model considers this boundary as permeable by introducing the contribution of both velocities of seepage flow and solid particles from the seabed to each velocity component in this boundary (i.e., the boundary condition (5)). In (5), the left side of the formula represents the fluid velocity on the seabed surface, while the terms on the right side of the formula represent the velocity generated by seepage and the velocity generated by soil particle movement, respectively. This can be obtained from the previous time-step seabed model.

2.2. Seabed Model

The interactions between the pore fluid and soil particles are governed by the quasi-static (Q-S) Biot equation. The relationship between soil displacement and pore water pressure in the seabed can be represented by the following equation:

(6)${\nabla }^2{p}_s-\frac{{\gamma }_w{n}_s{\beta }_s}{k_s}\frac{\partial p_s}{\partial t}=\frac{{\gamma }_w}{k_s}\frac{\partial }{\partial t}\left(\frac{\partial u_s}{\partial x}+\frac{\partial w_s}{\partial z}\right),$

(7)$G_s{\nabla }^2{u}_s+\frac{G_s}{1-2{\mu }_s}\frac{\partial }{\partial x}\left(\frac{\partial u_s}{\partial x}+\frac{\partial w_s}{\partial z}\right)=\frac{\partial p_s}{\partial x},$

(8)$G_s{\nabla }^2{w}_s+\frac{G_s}{1-2{\mu }_s}\frac{\partial }{\partial z}\left(\frac{\partial u_s}{\partial x}+\frac{\partial w_s}{\partial z}\right)=\frac{\partial p_s}{\partial z},$

where ${\beta }_s$ can be represented by the following equation:

(9)${\beta }_s=\frac{1}{K_w}+\frac{1-S_r}{P_{w0}},$

The inertia terms (acceleration of pore fluid and soil particles) are neglected in the Q-S model. As reported in the literature, based on the given wave and soil parameters, the dynamic constant (${\mathrm{\Pi }}_1$ and ${\mathrm{\Pi }}_2$) can be calculated. For conditions with a relative difference of less than 3% in dynamic constants, dynamic effects can be ignored and a quasi-static (Q-S) model can be directly used [25,26].

After solving the above governing equations, the soil displacements ($u_s$ and $w_s$) can be obtained. Then, the effective normal stresses (${\sigma }_x^{\prime }$ and ${\sigma }_z^{\prime }$) and the shear stresses (${\tau }_{xz}$ and ${\tau }_{zx}$) can be determined using the generalized Hook’s law.

(10)${\sigma }_x^{\prime }=2G_s\left[\frac{\partial u_s}{\partial x}+\frac{{\mu }_s}{1-2{\mu }_s}\left(\frac{\partial u_s}{\partial x}+\frac{\partial w_s}{\partial z}\right)\right],$

(11)${\sigma }_z^{\prime }=2G_s\left[\frac{\partial w_s}{\partial z}+\frac{{\mu }_s}{1-2{\mu }_s}\left(\frac{\partial u_s}{\partial x}+\frac{\partial w_s}{\partial z}\right)\right],$

(12)${\tau }_{xz}=G_s\left[\frac{\partial u_s}{\partial z}+\frac{\partial w_s}{\partial x}\right]={\tau }_{zx},$

The relation between total stress, effective stress, and pore pressure can be expressed as follows:

(13)$\sigma ={\sigma }^{\prime }-p,$

where $\sigma$ represents total stress; p represents pore pressure. Note that the tension normal stress is taken as a positive sign in this study.

By solving the governing equation, the pore water pressure and soil displacement within the seabed can be obtained. Therefore, it is necessary to specify several boundary conditions. It is generally believed that the effective normal stress and shear stresses on the seabed surface are zero, and the pore water pressure on the seabed surface is equal to the dynamic wave pressure at the fluid bottom:

(14)$p_s=p_w,{\sigma }_z^{\prime }={\tau }_{xz}=0\mathrm{at}z=0,$

At the bottom of the seabed ($z=-h_s$), an impermeable rigid boundary condition is used, where soil displacement and normal velocity gradient are set to zero:

(15)$u_s=w_s=\frac{\partial p_s}{\partial z}=0,$

Regarding the lateral boundary, the boundary conditions of no-flow and zero soil displacement are also used:

(16)$u_s=w_s=\frac{\partial p_s}{\partial x}=0,$

Meanwhile, the submarine pipeline has an impermeable rigid boundary, which means the no-flow boundary condition is applied to its surface:

(17)$\frac{\partial p_s}{\partial n}=0,$

According to the previous investigations [27], when there are no structures or structures are present, the transverse length of the seabed model can be simulated by taking two or three times the wavelength.

2.3. Two-Way Coupling Process

In the present model, two sets of governing equations, i.e., the fluid domain (RANS equations) and soil domain (Biot’s equation) can be linked by introducing an indicator function that distinguishes between different media. Then, the system could be considered as a whole while obtaining the information, such as effective stress, soil displacement, etc., and the governing equations for the two media could be solved through iterations until the system satisfies the equilibrium condition. That is,

(18)${\alpha }_s\left[\mathrm{EQ}.\left(1\right)-\left(3\right)\right]+(1-{\alpha }_s)\left[\mathrm{EQ}.\left(6\right)-\left(8\right)\right]=0,$

where ${\alpha }_s$ is an indicator factor from the flow and seabed domains. ${\alpha }_s=1$ denotes the flow domain while ${\alpha }_s=0$ denotes the seabed domain.

There are two ways to solve the above problem with the above coupling process. First, we can directly solve both domains simultaneously. This approach was attempted by Karunarthna and Lin [28] by describing the flow outside the soil by the Reynolds averaged Navier—Stokes equations and the flow inside soil by spatially averaged Navier—Stokes equations. However, the effective stress and soil deformation cannot be taken into account in their model [28]. This approach requires the same mesh density for both the flow and seabed models. Therefore, it is time-consuming and not easy to converge in some cases. Second, we can solve the whole system through iterations with the continuity of pressures and velocity at the seabed surface. In this case, we can have fine meshes in the flow model and less dense meshes in the seabed domain. This approach is sufficient for obtaining the flow characteristics and soil behaviors while considering the influences of soil properties (i.e., soil permeability) and soil response (i.e., pore pressure and soil displacement) on the flow field. In the present model, we adopt the second approach.

Figure 3 shows the flowchart of the two-way coupling. The procedure is outlined here:

• Steps 1: The VARANS equations, as shown in shown in (1)–(3), and the VOF equation, as shown in (4), are solved by the input wave parameters and boundary conditions. Then the flow velocity and water pressure in the whole computational domain are obtained. When the entire flow field stabilizes after a certain number of iterations, the dynamic wave pressure at the bottom of the fluid is extracted and transmitted as a boundary condition to the seabed surface. In the first cycle, assuming the flow velocity at the bottom of the fluid domain is zero, a one-way coupling algorithm is used. After the second wave cycle, the boundary condition at the seabed surface will be the updated velocity obtained from the seabed model, i.e., from Step 3. The coupling process shown in Figure 4a is illustrated.

• Steps 2: The dynamic wave pressures at the fluid bottom are extracted and transported to the seabed surface as boundary conditions. Then the quasi-static (Q-S) Biot equations, as shown in (6)–(8), are solved to calculate the pore pressure and soil displacements, as shown in Figure 4a.

• Steps 3: The flow velocities at the interface of the fluid and seabed are calculated by the pore pressure and soil displacements in the seabed, as shown in (5), and transported to the fluid bottom. The flow velocities at the fluid bottom are updated, as shown in Figure 4b. Then, the flow field will be modeled with the new boundary condition at $z=0$, and the effects of the seabed characteristics will be integrated into the flow model.

• Steps 4: To date, a complete wave–seabed two-way coupling process has been finished. We then proceed to the next time step and repeat steps 1–3.

Figure 3

The conceptual flow chart of the two-way coupling model.

[Figure omitted. See PDF]

Figure 4

Schematic diagram of fluid and seabed coupling. (a) Step no. 1 & no. 2; (b) Step no. 3.

[Figure omitted. See PDF]

2.4. Mesh Convergence and Stability of the Numerical Model

To ensure the computational accuracy of the numerical model, this paper conducted mesh convergence tests on the flow model and the seabed model. As shown in Figure 5 and Figure 6, it can be seen that when $L_w/240$ is used in the x-direction, both the fluid and seabed sub-models converge, while for the z-direction, when $H_w/30$ is used in the fluid sub-model and $h_s/250$ is used in the seabed sub-model, they converge.

The following example is used to examine the stability of the two-way coupling model used. It is worth noting that the parameters used in this validation are wave height $\left(H_w\right)$ = 4 m, period $\left(T_w\right)$ = 8 s, water depth $\left(d_w\right)$ = 10 m, seabed thickness $\left(h_s\right)$ = 15 m, Poisson’s ratio $\left({\mu }_s\right)$ = 0.3, soil permeability $\left(k_s\right)$ = 1 × 10${}^{-4}$ m/s, and seabed porosity $\left(n_s\right)$ = 0.425. Figure 7 shows the time history of the wave surface variation under the two-way coupling algorithm. This study ran for 30 wave cycles. From the figure, it can be clearly seen that when the model ran for 48 s; that is, after running for the 6th cycle, the wave surface time history curve gradually stabilized and remained relatively stable until the 30th cycle. This demonstrates the stability of the present numerical model.

3. Model Validation

This paper validates the proposed model through the following experimental data, the wave and soil parameters considered in model validation are shown in Table 1:

• Validation no. 1: Comparison of the experimental results between the present model and Sun et al. [29] for a single pipe with a trench layer; in this experiment, only the pore pressures along the pipe surface were measured.

• Validation no. 2: Comparison of the experimental results between the present model and Zhai et al. [18] for twin pipes in tandem.

• Validation no. 3: Comparison of the experimental results between the present model and Chen et al. [30] for a single pipe with a trench layer; in this experiment, in addition to the pore pressures along the pipeline surface, additional measurements of pore pressures below the pipeline were taken.

Table 1

Parameters in the model validation.

${\mathit{H}}_{\mathit{w}}$ (m)	${\mathit{d}}_{\mathit{w}}$ (m)	${\mathit{T}}_{\mathit{w}}$ (s)	${\mathit{\mu }}_{\mathit{s}}$	${\mathit{k}}_{\mathit{s}}$ (m/s)	${\mathit{n}}_{\mathit{s}}$	D (m)
Sun et al. [29]
0.14	0.4	1.4	0.32	3.56 × 10${}^{-5}$	0.396	0.1
Zhai et al. [18]
0.1	0.4	1.2	0.3	3.56 × 10${}^{-5}$	0.369	0.1
Chen et al. [30]
0.08	0.4	1.2	0.3	3.57 × 10${}^{-5}$	0.369	0.08

3.1. Validation No. 1: A Single Pipeline in the Trench Layer

The first model validation is to compare the numerical results obtained by the present two-way coupling model with the experimental results of Sun et al. [29]. The experiment used a pipeline with a diameter of 0.1 m to study the pore pressure around buried pipelines in trenches under the wave action. As shown in Figure 8, eight pore pressure sensors were installed around the pipeline in this experiment to better observe the results.

Figure 9 shows a comparison of the maximal pore pressure around the pipeline in Test 10 between the numerical model and Sun et al. [29]’s experimental data. This validation selects the working condition where the pipeline is fully buried in a trench with a depth of 0.15 m. Figure 10 further compares the maximal pore pressure values at different measuring points $(z/h_s=-0.411$ and $z/h_s=-0.482$) under the pipeline. It is clear that the results obtained from the two-way coupling model in this study are basically consistent with the experimental results of Sun et al. [29]. The results indicate that the present numerical model can accurately predict the wave-induced pore pressure around the pipeline.

Figure 8

Schematic diagram of the experimental setup [29] for validation no. 1.

[Figure omitted. See PDF]

3.2. Validation No. 2: Twin Pipelines In Tandem

Recently, Zhai et al. [18] conducted a series of experimental studies to investigate the wave-induced pore pressure around twin pipelines in series. As shown in Figure 11, the length and outer diameter of the twin pipelines are 1 m and 0.12 m, respectively. Eight pore pressure sensors were uniformly installed in the center of the circumference-measuring pipeline.

Figure 12 compares the simulated excess pore pressure $(\mathsf{\Delta }u=p_b-p_s)$ under the fully $(d_t=D_m)$ buried twin pipelines at P3 (referring to Figure 11) with the measured results of [18], where $p_b$ represents the dynamic wave pressure extracted from the wave model, and $p_s$ denotes the wave-induced oscillatory pore pressure. The distance between twin pipelines in this experiment is 0.24 m. Figure 13 displays the comparison of the amplitude of the excess pore pressure around the upstream pipeline between the present model and the experimental data reported by Zhai et al. [18]. As seen from the figure, there is a slight phase difference between the numerical results and experimental results, but the overall agreement is acceptable.

3.3. Validation No. 3: A Single Pipeline in the Trench Layer

The third model validation compares the present model with the experimental results of Chen et al. [30]. As shown in Figure 14, four pore pressure sensors were installed in the seabed to measure the pore pressure around the pipeline and below the pipeline for 3 cm (P5), 8 cm (P6), and 18 cm (P7).

Figure 15 shows the changes in pore pressure around a fully buried pipeline; as seen from the figure, a slight phase difference between the numerical model and the experimental results can be found, but the overall agreement is acceptable. In summary, the present two-way coupling numerical model has the capability to simulate the pore pressures around the pipelines.

Figure 14

Schematic diagram of the experimental setup [30] for validation no. 3.

[Figure omitted. See PDF]

4. Results and Discussion

In this section, we will discuss the impacts of different wave and seabed parameters on the pore pressure around pipelines. All wave, seabed, and pipeline parameters for the study are listed in Table 2. Numerical examples involve twin pipelines with different diameters on the porous seabed. G represents the horizontal distance between twin pipelines in series.

4.1. Effect of Pipeline Configuration (Single Pipeline and Twin Pipelines)

In this section, we compare the response of seabed soil around single pipelines and a series of twin pipelines under wave actions. It is worth noting that the G between the twin pipelines is 5 m here, and only the upstream pipelines of the twin pipelines are used for the comparison with the single pipeline.

Figure 16 shows the horizontal and vertical velocities ($u_f$ and $w_f$) at the seabed surface of fully buried single and twin pipelines. The figure shows that the velocity of the seabed surface varies in both cases. For twin pipelines, the maximum $u_f$ can reach 0.3 mm/s, and the maximum $w_f$ is slightly less than 0.4 mm/s. In contrast, the maximum $u_f$ is 0.2 mm/s, and the maximum $w_f$ is slightly over 0.4 mm/s for a single pipeline. It can be seen that the horizontal velocity of the twin pipelines is greater than the single pipeline, while the vertical velocity has the opposite result.

Figure 17a shows the maximal pore pressure $\left(\right|p_s|/p_0)$ around the pipeline. It can be clearly seen that $|p_s|/p_0$ around the twin pipeline is slightly lower than that around the single pipeline under the two-way coupling algorithm. Then the variation of $|p_s|/p_0$ along the seabed depth through the center of the pipeline is displaced in Figure 17b. It is observed that $|p_s|/p_0$ below the twin pipelines is slightly lower compared with the results of the single pipeline.

4.2. Effect of Steady Current

According to previous research, the combined effect of waves and currents has significant changes in the flow field and soil response [31]. Therefore, the dynamic response and liquefaction of the seabed around twin pipelines under the combined action of waves and currents will become more complex.

The comparison of pore pressure $\left(\right|p_s|/p_0)$ below the pipeline along the seabed depth under different current velocities $\left(U_c\right)$ is represented in Figure 18. The steady current is considered in the flow model by prescribing a constant horizontal velocity on the inlet and outlet boundaries as well as the internal fields. As seen from the figure, taking the upstream pipeline as an example, when the flow velocity reaches 2 $\mathrm{m}\mathrm{/}\mathrm{s}$, the $|p_s|/p_0$ is the largest, while the $|p_s|/p_0$ obtained at the same position is the smallest when the reverse current reaches $-2$ $\mathrm{m}\mathrm{/}\mathrm{s}$. This indicates that the forward current velocity will cause an increase in the $|p_s|/p_0$ below the pipeline. On the contrary, the $|p_s|/p_0$ will decrease with the increasing reverse current velocity. This result is consistent with the effects of the current on the dynamic wave pressures [22].

4.3. Effects of Wave and Seabed Characteristics

In this section, wave and seabed parametric analyses are conducted to discuss their influences on the wave and seabed characteristics around the pipeline.

The influence of $H_w$ on pore pressure $\left(\right|p_s|/p_0)$ around the pipeline is shown in Figure 19. It is obvious that $|p_s|/p_0$ in the seabed is positively correlated with the change in $H_w$, which means that $|p_s|/p_0$ would increase with the increase of $H_w$. At the same time, the $|p_s|/p_0$ along the depth direction below the pipeline also increases with the increase of $H_w$, but the change is not insignificant.

Figure 20 shows the influence of $d_w$ on the seabed response around the pipeline. Similarly, there is a positive correlation between the pore pressure ($|p_s|/p_0$) around the pipeline and the $d_w$, which means that the $|p_s|/p_0$ around the pipeline increases with an increase in $d_w$. Figure 20c,d show the variation of $|p_s|/p_0$ along the depth below the pipeline. It can be clearly seen from the figure that when $z/h_s$ is greater than 0.3, the difference in $|p_s|/p_0$ caused by $d_w$ gradually becomes apparent.

Figure 21 displays the maximal pore pressures ($|p_s|/p_0$) around the pipeline obtained from different $T_w$. It can be seen that $|p_s|/p_0$ around the pipeline increases with the decrease in $T_w$. The difference is most pronounced below the pipeline, i.e., ($\theta$ = ${270}^{\circ }$). It is worth noting that when $z/h_s$ is less than 0.4, $|p_s|/p_0$ in the upstream and downstream pipelines decreases with depth as the period increases. When $z/h_s$ is greater than 0.4 until the bottom of the seabed, $|p_s|/p_0$ in the upstream and downstream pipelines increases in depth as the period increases.

For soil characteristics, as reported in the literature for the case with a single pipeline [4], $k_s$ and $S_r$ will have notable influences on the responses caused by waves in the seabed. Therefore, this study will discuss the influences of different $k_s$ and $S_r$ on the pore pressure around pipelines.

Figure 22 illustrates the significant influence of $k_s$ on the maximal pore pressure ($|p_s|/p_0$) around the twin pipelines. The $|p_s|/p_0$ around the pipeline would decrease with the $k_s$. For $k_s$=1 × 10${}^{-5}$ m/s, $|p_s|/p_0$ continuously decreases around the upstream pipeline circumference (${180}^{\circ }$ < $\theta$ < ${360}^{\circ }$). Figure 23 shows that the increase in the degree of saturation $\left(S_r\right)$ causes the decrease of $|p_s|/p_0$ around the twin pipelines. The above two figures illustrate that $k_s$ and $S_r$ greatly affect the distribution of pore pressure around pipelines.

Figure 20

Effects of different $d_w$ on the pore pressure around the pipelines. (a) Upstream pipeline; (b) downstream pipeline; (c) upstream pipeline; (d) downstream pipeline.

[Figure omitted. See PDF]

Figure 21

Effects of different $T_w$ on the pore pressure around the pipelines. (a) Upstream pipeline; (b) downstream pipeline; (c) upstream pipeline; (d) downstream pipeline.

[Figure omitted. See PDF]

Figure 22

Effects of different $k_s$ on the pore pressure around the pipelines. (a) Upstream pipeline; (b) downstream pipeline; (c) upstream pipeline; (d) downstream pipeline.

[Figure omitted. See PDF]

Figure 23

Effects of different $S_r$ on the pore pressure around the pipelines. (a) Upstream pipeline; (b) downstream pipeline; (c) upstream pipeline; (d) downstream pipeline.

[Figure omitted. See PDF]

4.4. Liquefaction Assessment around Fully Buried Pipelines

4.4.1. Seabed Consolidation

The seabed undergoes consolidation due to its own weight; the initial consolidation will have an impact on the seabed response around the marine structure. Therefore, it is necessary to determine the initial effective stress inside the seabed when analyzing the soil liquefaction. Figure 24 illustrates the $|{\sigma }_0^{\prime }|$ due to differences in the pipeline configuration. It can be indicated that the own weight of the pipeline would significantly impact $|{\sigma }_0^{\prime }|$ and further impact the stability of the pipeline.

The modified liquefaction standard based on the average normal effective stress is represented as [32]:

(19)$p_s-p_b\ge {\sigma }_0^{\prime }=\frac{{\sigma }_{x0}^{\prime }+{\sigma }_{z0}^{\prime }}{2},$

4.4.2. Distribution of Liquefaction Zones under Various Conditions

Figure 25 shows the impacts of different wave parameters on the liquefaction distribution around pipelines. The figure shows that as $H_w$, $T_w$, and $d_w$ increase, the liquefaction range of the seabed also increases, although the trend is not very obvious. Therefore, when $H_w$ and $T_w$ are large, the seabed soil is more prone to liquefaction.

Figure 26 represents the influence of different soil parameters on the liquefaction range of soil. It can be seen from the figure that smaller $S_r$ and $k_s$ result in a larger liquefaction range. Therefore, liquefaction is more likely to occur around pipelines with lower $S_r$ and $k_s$.

Figure 25

Distribution of the liquefied zone near submarine pipelines under different wave characteristics. (a) $H_w$; (b) $d_w$; (c) $T_w$.

[Figure omitted. See PDF]

4.4.3. Effect of Pipeline Configuration

This section mainly investigates the effects of G and $R_1/R_2$ on the response of soil around pipelines and seabed liquefaction. The impacts of different G values on the seabed response around the pipeline are shown in Figure 27. Wherein, Figure 27a,b show the distribution of $|p_s|/p_0$ along the lower surface of the pipelines (${180}^{\circ }$ < $\theta$ < ${360}^{\circ }$), respectively. The $|p_s|/p_0$ around the pipeline increases with G. Meanwhile, Figure 27c,d show the variation of $|p_s|/p_0$ below the pipeline with seabed depth, respectively. Similarly, $|p_s|/p_0$ increases with G.

To demonstrate the influence of G on the seabed liquefaction range, Figure 28 shows the seabed liquefaction near the pipelines, with four different gaps between the pipelines. It can be clearly seen that the larger the spacing between the two pipelines, the larger the liquefaction range around the pipeline. However, the liquefaction range of the upstream and downstream pipelines would be independent and similar to that of a single pipeline when G > 7.5 m.

Figure 29 represents the influences of the pipeline radius ratio ($R_1/R_2$) on the seabed response around the pipeline. It is noted that the upstream pipeline diameter is fixed at 2 m, while the downstream pipeline diameter ranges from 1 to 3m, respectively. As shown in the figure, $|p_s|/p_0$ around the upstream pipeline would decrease with the increasing $R_1/R_2$, while $|p_s|/p_0$ from the downstream pipeline to the surrounding area is the opposite. As the diameter of the downstream pipeline decreases, the influence of the upstream pipeline on $|p_s|/p_0$ around the downstream pipeline increases.

Figure 27

Effects of different G values on the pore pressure around pipelines. (a) Upstream pipeline; (b) downstream pipeline; (c) upstream pipeline; (d) downstream pipeline.

[Figure omitted. See PDF]

Figure 28

Distribution of the liquefied zone near submarine pipelines under different gap values. (a) G = 2.5 m; (b) G = 5 m; (c) G = 7.5 m; (d) G = 10 m.

[Figure omitted. See PDF]

[Figure omitted. See PDF]

Figure 29

Effects of different $R_1/R_2$ values on the pore pressure around pipelines. (a) Upstream pipeline; (b) downstream pipeline; (c) upstream pipeline; (d) downstream pipeline.

[Figure omitted. See PDF]

To demonstrate the impact of the pipeline radius ratio $(R_1/R_2)$ between twin pipelines on the liquefaction range near the pipeline, Figure 30 shows the seabed liquefaction near the pipeline; there are three different pipeline radius ratios between the pipelines. It should be noted that the pipeline diameter of the upstream pipeline is fixed, and only the downstream pipeline changes. From the figure, it can be seen that as $R_1/R_2$ increases, the liquefaction range near the upstream pipeline remains almost unchanged, while the liquefaction range near the downstream pipeline gradually decreases. When $R_1/R_2$ = 2, the liquefaction depth above the downstream pipeline is the smallest, but there is no liquefaction phenomenon under the pipeline.

5. Conclusions

This article adopts a two-way coupling algorithm to investigate the liquefaction potential around twin pipelines under the combined action of waves and currents. To verify the effectiveness of this model, a series of comparisons was also conducted with previous experimental data. At the same time, the influences of different wave–current parameters, soil parameters, and pipeline configurations to the response of seabed soil were analyzed. Based on the numerical examples presented above, it can be concluded that:

• This model adopts a two-way coupling algorithm to study the seabed response around pipelines. Experimental verification shows that the model is effective.

• The liquefaction range is greatly influenced by $H_w$ and $T_w$. The results prove that when $H_w$ and $T_w$ are small, liquefaction hardly occurs.

• The soil characteristics greatly affect the seabed response and liquefaction zone distribution around the pipeline. This study mainly explores the effects of $S_r$ and $k_s$. The liquefaction hardly occurs when $S_r$ and $k_s$ are large enough.

• The pipeline configurations, including the gap $\left(G\right)$ and pipeline radius ratio $(R_1/R_2)$, also have obvious effects on the liquefaction zone around the twin pipelines. The liquefaction range around the pipeline gradually increases with the gap $\left(G\right)$. However, the liquefaction range of the upstream and downstream pipelines is independent and similar to that of a single pipeline when G > 7.5 m.

• As the pipeline radius ratio $(R_1/R_2)$ gradually increases, the liquefaction range above the downstream pipeline gradually decreases, and no liquefaction occurs below the pipeline, so the pipeline is relatively stable.

There are a few limitations of the present model that could be further improved in future work. Only the oscillatory soil response and the corresponding transient liquefaction were considered in this paper; a more advanced soil constitutive model that could capture the elasto-plastic behavior of the seabed is desired and the residual mechanism and residual liquefaction could be simulated.

Author Contributions

J.Z.: methodology; validation; data curation; visualization; and writing—original draft; L.C.: methodology; validation, writing—review and editing; supervision; H.Z.: methodology; validation; data curation; writing—review and editing; D.-S.J.: conceptualization; funding; supervision; writing—review and editing. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available upon request from the corresponding author.

Acknowledgments

The authors thank Jisheng Zhang from Hohai University for providing the experimental data used to validate the present model.

Conflicts of Interest

The authors declare no conflict of interest.

Notation

$E_s$	Young’s module
G	gap between pipes
$G_s$	shear modules
$H_w$	wave height
$K_w$	true bulk modulus of water
$P_{w0}$	absolute pressure
$R_1$	pipeline radius
$R_2$	pipeline radius
$R_1/R_2$	pipeline radius ratio
$S_r$	degree of saturation
$T_w$	wave period
$\mathit{U}$	velocity vector
$U_c$	current velocity
$d_w$	water depth
$e_1$	burial depth
$e_2$	burial depth
g	gravitational acceleration
$h_s$	seabed thickness
$k_s$	permeability
$n_s$	porosity
$p_s$	fluid pressure
$p_b$	hydrodynamic pressure acting on the seabed surface
t	Time
${\mathit{u}}_{\mathit{r}}$	relative velocity field
${\beta }_s$	compressibility of pore fluid
${\gamma }_w$	unit weight of water
${\mathbf{u}}_{\mathbf{s}}$	soil displacement vector
${\mu }_s$	Poisson’s ration
${\sigma }_0^{\prime }$	initial stresses
$\tau$	viscous stress vector

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Concepts of coupling models. (a) One-way coupling model; (b) two-way coupling model.

Figure 2. The sketch of wave–seabed–pipelines in tandem.

Figure 5. Maximum free surface elevation under different mesh densities. (a) Mesh distribution in the x-direction; (b) mesh distribution in the z-direction.

Figure 6. Maximum pore pressure [Forumla omitted. See PDF.] under different mesh densities. (a) Mesh distribution in the x-direction; (b) mesh distribution in the z-direction.

Figure 7. Time series of wave surface changes.

Figure 9. Comparison of pore pressure around the pipeline with experimental results [29] for test 10.

Figure 10. Comparison of wave-induced pore pressure with experimental results [29] for Test 10 at different depths: (a) [Forumla omitted. See PDF.] and (b) [Forumla omitted. See PDF.].

Figure 11. Schematic diagram of the experimental setup [18] for validation no. 2.

Figure 12. Comparison of excess pore pressure [Forumla omitted. See PDF.] at P3 with experimental results [18].

Figure 13. Comparison of the excess pore pressure amplitude [Forumla omitted. See PDF.] around the upstream pipeline in tandem with the experimental data [18].

Figure 15. Comparison of pore pressure around the pipeline with experimental results [30].

Figure 16. Comparison of seabed surface velocity between the single and twin pipelines. (a) [Forumla omitted. See PDF.]; (b) [Forumla omitted. See PDF.].

Figure 17. Distributions of the maximum pore pressure [Forumla omitted. See PDF.]. (a) around the pipeline; (b) along the vertical direction below the pipeline.

Figure 18. Effects of different [Forumla omitted. See PDF.] on the pore pressure below the pipeline. (a) Upstream pipeline; (b) downstream pipeline.

Figure 19. Effects of different [Forumla omitted. See PDF.] on the pore pressure around the pipelines. (a) Upstream pipeline; (b) downstream pipeline; (c) upstream pipeline; (d) downstream pipeline.

Figure 24. The distribution of the initial average effective stress [Forumla omitted. See PDF.] in the consolidated state for various pipeline configurations. (a) Single pipeline; (b) twin pipelines.

Figure 26. Distribution of the liquefied zone near submarine pipelines under different soil characteristics. (a) [Forumla omitted. See PDF.]; (b) [Forumla omitted. See PDF.].

Figure 30. Distribution of the liquefaction zone near submarine pipelines under different pipeline radius ratios. (a) [Forumla omitted. See PDF.] = 2/3; (b) [Forumla omitted. See PDF.] = 1; (c) [Forumla omitted. See PDF.] = 2.

Table 2

Parameters used in the numerical simulation.

Characteristics	Value	Unit
Wave characteristics
$U_c$	0 or various	[m/s]
$H_w$	4.0 or various	[m]
$d_w$	10.0 or various	[m]
$T_w$	8.0 or various	[s]
Seabed characteristics
$k_s$	1 × 10${}^{-4}$ or various	[m/s]
${\mu }_s$	0.3	-
$n_s$	0.425	-
$S_r$	0.98 or various	-
$G_s$	1 × 10${}^7$	[N/m${}^2$]
$E_s$	2.6 × 10${}^7$	[N/m${}^2$]
$h_s$	15	[m]
Pipeline characteristics
$R_1$	1	[m]
$R_2$	0.5–1.5	[m]
$e_1$	2	[m]
$e_2$	2.5 or various	[m]
G	5 or various	[m]

References

1. Sumer, B.M. Liquefaction around Marine Structures; World Scientific: Hackensack, NJ, USA, 2014; [DOI: https://dx.doi.org/10.1142/7986]

2. Lin, Z.; Guo, Y.; Jeng, D.S.; Liao, C.; Rey, N. An integrated numerical model for wave–soil–pipeline interactions. Coast. Eng.; 2016; 108, pp. 25-35. [DOI: https://dx.doi.org/10.1016/j.coastaleng.2015.11.003]

3. Duan, L.; Liao, C.; Jeng, D.S.; Chen, L. 2D numerical study of wave and current-induced oscillatory non-cohesive soil liquefaction around a partially buried pipeline in a trench. Ocean Eng.; 2017; 135, pp. 39-51. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2017.02.036]

4. Liang, Z.; Jeng, D.S.; Liu, J. Combined wave–current induced seabed liquefaction around buried pipelines: Design of a trench layer. Ocean Eng.; 2020; 212, 107764. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2020.107764]

5. Afzal, M.S.; Holmedal, L.E.; Myrhaug, D. Sediment transport in combined wave–current seabed boundary layers due to streaming. J. Hydraul. Eng.; 2021; 147, 04021007. [DOI: https://dx.doi.org/10.1061/(ASCE)HY.1943-7900.0001862]

6. Afzal, M.S.; Holmedal, L.E.; Myrhaug, D. Effect of Wave Skewness and Sediment Particle Size on Sediment Transport Due to Combined Wave–Current Seabed Boundary Layer Streaming. J. Hydraul. Eng.; 2022; 148, 06022011. [DOI: https://dx.doi.org/10.1061/(ASCE)HY.1943-7900.0002007]

7. Dutta, D.; Bihs, H.; Afzal, M.S. Computational Fluid Dynamics modelling of hydrodynamic characteristics of oscillatory flow past a square cylinder using the level set method. Ocean Eng.; 2022; 253, 111211. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2022.111211]

8. Dutta, D.; Afzal, M.S.; Alhaddad, S. 3D CFD study of scour in combined wave–current flows around rectangular piles with varying aspect ratios. Water; 2023; 15, 1541. [DOI: https://dx.doi.org/10.3390/w15081541]

9. Jeng, D.S. Porous Models for Wave-Seabed Interaction; Springer: Berlin/Heidelberg, Germany, 2012; [DOI: https://dx.doi.org/10.1007/978-3-642-33593-8]

10. Yamamoto, T.; Koning, H.; Sellmeijer, H.; Hijum, E.V. On the response of a poro-elastic bed to water waves. J. Fluid Mech.; 1978; 87, pp. 193-206. [DOI: https://dx.doi.org/10.1017/S0022112078003006]

11. Hsu, J.R.C.; Jeng, D.S. Wave-induced soil response in an unsaturated anisotropic seabed of finite thickness. Int. J. Numer. Anal. Methods Geomech.; 1994; 18, pp. 785-807. [DOI: https://dx.doi.org/10.1002/nag.1610181104]

12. Seed, H.B.; Rahman, M.S. Wave-induced pore pressure in relation to ocean floor stability of cohesionless soils. Mar. Geotechnol.; 1978; 3, pp. 123-150. [DOI: https://dx.doi.org/10.1080/10641197809379798]

13. Sumer, B.M.; Fredsøe, J. The Mechanics of Scour in the Marine Environment; World Scientific: Hackensack, NJ, USA, 2002; [DOI: https://dx.doi.org/10.1142/4942]

14. Rados, K.; Pitt, D.; Macfarlane, D.; Dimla, D. On the interaction of pipelines on the seabed. Proceedings of the Tenth International Offshore and Polar Engineering Conference; Seattle, WA, USA, 28 May–2 June 2000; pp. 155-161.

15. Zhao, M.; Cheng, L. Numerical modeling of local scour below a piggyback pipeline in currents. J. Hydraul. Eng. ASCE; 2008; 134, pp. 1452-1463. [DOI: https://dx.doi.org/10.1061/(ASCE)0733-9429(2008)134:10(1452)]

16. Jo, C.H.; Shin, Y.S.; Min, K.H. Numerical and Experimental Studies of Dual Subsea Pipelines in Trench. J. Ship Ocean Technol.; 2002; 6, pp. 12-22.

17. Zhang, J.; Zheng, J.; Zhang, C.; Jeng, D.S.; Guo, Y. Numerical study on the interaction between waves and twin pipelines in sandy seabed. J. Coast. Res.; 2013; 65, pp. 428-433. [DOI: https://dx.doi.org/10.2112/SI65-073.1]

18. Zhai, Y.; Zhang, J.; Guo, Y.; Tang, Z.; Zhang, T. Study of wave-induced seabed response around twin pipelines in sandy seabed through laboratory experiments and numerical simulations. Ocean Eng.; 2022; 244, 110344. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2021.110344]

19. Chen, X.; Jeng, D.S.; Liang, Z.D. Numerical analysis of seabed liquefaction in the vicinity of two tandem pipelines in a trench. Ocean Eng.; 2022; 266, 112656. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2022.112656]

20. Gao, Y.; Lin, J.; Zhang, J.; Guo, Y. Numerical modeling of combined wave and current-induced residual liquefaction around twin pipelines. Ocean Eng.; 2022; 261, 112134. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2022.112134]

21. Zhai, H.; Jeng, D.S. Two-way coupling model for wave-induced oscillatory soil response in a porous seabed. Ocean Eng.; 2022; 249, 110791. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2022.110791]

22. Zhang, J.S.; Zhang, Y.; Jeng, D.S.; Liu, P.F.; Zhang, C. Numerical simulation of wave–current interaction using a RANS solver. Ocean Eng.; 2014; 75, pp. 157-164. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2013.10.014]

23. Higuera, P.; Lara, J.L.; Losada, I.J. Realistic wave generation and active wave absorption for Navier-Stokes models: Application to OpenFOAM. Coast. Eng.; 2013; 71, pp. 102-118. [DOI: https://dx.doi.org/10.1016/j.coastaleng.2012.07.002]

24. Liang, Z.D.; Jeng, D.S. PORO–FSSI–FOAM model for seafloor liquefaction around a pipeline under combined random wave and current loading. Appl. Ocean Res.; 2021; 107, 102497. [DOI: https://dx.doi.org/10.1016/j.apor.2020.102497]

25. Jeng, D.S.; Cha, D. Effects of dynamic soil behavior and wave non-linearity on the wave-induced pore pressure and effective stresses in porous seabed. Ocean Eng.; 2003; 30, pp. 2065-2089. [DOI: https://dx.doi.org/10.1016/S0029-8018(03)00070-2]

26. Ulker, M.B.C.; Rahman, M.S.; Jeng, D.S. Wave-induced response of seabed: Various formulations and their applicability. Appl. Ocean Res.; 2009; 31, pp. 12-24. [DOI: https://dx.doi.org/10.1016/j.apor.2009.03.003]

27. Ye, J.H.; Jeng, D.S. Response of Porous Seabed to Nature Loadings: Waves and Currents. J. Eng. Mech.; 2011; 138, pp. 601-613. [DOI: https://dx.doi.org/10.1061/(ASCE)EM.1943-7889.0000356]

28. Karunarthna, K.; Lin, P. Numerical simulation of wave damping over porous seabeds. Coast. Eng.; 2006; 53, pp. 845-855. [DOI: https://dx.doi.org/10.1016/j.coastaleng.2006.05.003]

29. Sun, K.; Zhang, J.; Gao, Y.; Jeng, D.S.; Guo, Y.; Liang, Z. Laboratory experimental study of ocean waves propagating over a partially buried pipeline in a trench layer. Ocean Eng.; 2019; 173, pp. 617-627. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2019.01.022]

30. Chen, H.; Zhang, J.; Tong, L.; Sun, K.; Guo, Y.; Wei, C. Experimental study of soil responses around a pipeline in a sandy seabed under wave-current load. Appl. Ocean. Res.; 2023; 130, 103409. [DOI: https://dx.doi.org/10.1016/j.apor.2022.103409]

31. Zhang, Y.; Jeng, D.S.; Gao, F.; Zhang, J.S. An analytical solution for response of a porous seabed to combined wave and current loading. Ocean Eng.; 2013; 57, pp. 240-247. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2012.09.001]

32. Zhao, H.Y.; Jeng, D.S.; Guo, Z.; Zhang, J.S. Two-dimensional model for pore pressure accumulations in the vicinity of a buried pipeline. J. Offshore Mech. Arct. Eng. ASME; 2014; 136, 042001. [DOI: https://dx.doi.org/10.1115/1.4027955]