The Role of 2D Seepage on Sediment Incipient Motion around a Pipeline

Pipelines have been used as one of the main transportation methods for the offshore industry, with increasing activities in marine resources recently. Prediction of seabed instability is one of key factors that must be taken into consideration for an offshore pipeline project. As the first step of the scour process, sediment incipient motion has been intensively studied in the past. Most previous investigations didn’t consider the wave-induced seepage in the elevation of sediment motion. In this paper, two-dimensional seepage was considered to modify the conventional Shields number and its associated impact on sediment incipient motion around the trenched pipeline was investigated. Both flat and sloping seabeds are considered. The numerical results indicated that a peak or valley of the modified Shields number was formed below the pipeline and horizontal seepage flow tremendously impact the sediment motion in the vicinity of the pipeline. Parametric analysis concludes: the influence of the seepage around the pipeline becomes more significant in a large wave, shallow water in a seabed with large shear modulus and permeability, and larger pipeline diameter and smaller flow gap ratio. This will make soil particles be more easily dragged away from the seabed.


Introduction
Offshore Pipelines have been one of key offshore installations for the offshore oil and gas industry, due to the growth of development of marine resources. In the marine environments, the existence of a pipeline changes the flow pattern and promotes the ability of sediment transport, which further affects the scour process [1][2][3]. The development of scour beneath pipelines can cause a free spanning pipeline that will cause further structural damage.
In general, sediment incipient motion occurs at the beginning stage of local scour around marine installations [4]. Sediment threshold is not only important for the general problems of sediment transport and scour, but also impacts the scour protection and the design of mitigation measures in engineering. Hence, it has been extensively studied by scientists and engineers [5]. Shields number originally proposed by Shields [6] has been widely used as an index for the sediment to move from the static condition [4].
Since Shields [6], numerous modifications of conventional Shields number have been proposed, based on a series of laboratory tests and field measurements. Among these, Madsen and Grant [7] determined the maximum bottom shear stress under an oscillatory fluid motion according to the wave friction factor concept [8] and proposed the Shields criterion for oscillatory unsteady flow. Whitehouse [9] and Chiew and Parker [10] investigated the threshold condition for the onset of cohesionless sediment transport on a sloping seabed. Later, Juez et al. [11] proposed a new Shields number considering gravity currents over steep slopes and investigated the bed load transport at the action of gravity currents through both experimental and numerical models, based on their previous studies for bed load transport over the sloping river [12][13][14]. However, only the basic forces (i.e., included submerged weight force, left and drag force) were included, and the seabed surface was assumed as an impermeable and rigid boundary in their studies, in which the seepage force was ignored.
In general, the seepage velocity in a seabed is rather small, compared with the freestream velocity [15]. Therefore, it has been considered as a minor effect on the flow fields. However, its associated change with the hydrodynamic force on the seabed surface could be visible for some cases with a structure [16]. Hence, it is necessary to consider the impact of seepage to the sediment motion. For example, as pointed out by Cheng [17] and Sumer et al. [18,19], seepage has significant effects on sediment incipient motion based on their experiments. Qi and Gao [20] reported that the upward seepage brings the sand-bed more susceptible to scouring based on their wave flume experiments. In addition to the experimental studies, Niven [21] theoretically examined the onset of sediment motion with upward seepage. Recently, Yang [22] further proposed a formulas for sediment transport due to vertical flows.
To date, only a few studies have considered seepage on incipient sediment motion quantitatively. Cheng and Chiew [23] considered the upward seepage in the force analysis for the threshold condition of sediment incipient motion and verified their theoretical model by a series of laboratory experiments. Recently, Guo et al. [24] further proposed a new Shields number with the upward seepage, and integrated it into their numerical model (COMSOL model) to investigate the impact of upward seepage under an oscillatory flow. Later, Li et al. [16] further included the upward seepage in their CFD scour model (OpenFOAM model). However, these studies only considered the upward seepage, i.e., onedimensional. In fact, the real scientific problem should be either 2D or 3D.
The combination of scour and soil response has always been a major difficulty in the field of marine engineering [25]. The purpose of this study is to make a preliminary connection of scour and soil response by considering seepage on sediment incipient motion, which would be further drawn into the process of scour in future research. By the way, the process of scour around marine structures would be more accurately predicted, which has a great effect on the scour protection and the design of mitigation measures.
In this paper, we proposed a new modified Shields number by considering 2D wave (current)-induced seepage. The paper is organized as follows: in Section 2, the previous numerical model (PORO-FSSI-FOAM, [26,27]) for the wave/current induced soil response in a porous seabed in the vicinity of a submarine pipeline is outlined. In Section 3, the new Shields number will be derived. Then, the effects of 2D seepage on sediment motion will be investigated. In Section 4, comparisons between the new Shields number with the conventional Shields number will be presented through the parametric study. In Section 5, the proposed model will be applied in the cases with a submarine pipeline in a trenched layer. Finally, the key finding will be summarized in Section 6.

CFD Model (PORO-FSSI-FOAM)
The problem of wave-seabed interactions in the vicinity of a submarine pipeline has been studied since the 1980s, based on the poro-elastic theory including wave experimental studies [28] and numerical modeling [29,30] for wave loading. Recently, some studies considered both wave and current loading and the pipeline in a trenched layer [26,31].
In this study, the porous model for the fluid-seabed-structure interactions (PORO-FSSI-FOAM) proposed in the previous study [26] is used to determine the wave/current induced seepage within the porous seabed with a pipeline. In the present model, the two-phase fluid flow domain above the seabed is simulated by olaflow later [32], while the porous seabed is simulated by the Biot Quasi-Static (QS) model [26,33].
The problem of fluid-seabed-pipeline interaction considered here is depicted in Figure 1. As shown in the figure, the pipeline is located at the position about 25D from the inlet and h above the seabed surface, which is spanned in the flow field. In the flow model, the two-phase fluids model of air and water is established to simulate the wave propagation and the wave pressure, since the only difference between the present study and the previous work [34] is the inclusion of pipeline and uniform currents. Herein, we do not repeat the governing equations for the olaflow. The details of the governing equations of olaflow can be found in Higuera et al. [32], Higuera [35], Zhai et al. [34].
For the seabed model, the Quasi-Static(QS) model [36] is employed in the soil-pipeline model to describe the soil response. In the model, the soil skeleton is assumed as an isotropic elastic material.
Based on the conservation of mass and momentum, the governing equations for the flat seabed model are listed here: where p s is the pore pressure of the seabed; and β s is the compressibility of pore fluid, which is related to degree of saturation (S r ); u s and w s are the soil displacements in xand zdirections, respectively; n s is the soil porosity; k s is the permeability of the seabed; γ w is the unit weight of pore water; µ s is the Poisson ratio, and G s is the shear modulus of soil. The above governing equations for the seabed model, (1)-(3), are established in global coordinates, which should be redefined in local coordinates for a slop seabed. The sketches of numerical domain are displayed in Figure 1, and both global (x,z) and local (x ,z ) axes are defined. In the model, the projected gravity in global coordinates is adopted [11]: in which β is the angle of slop seabed and cs = cos β is introduced to simplify the expression in the following equations. By considering a control volume with vertical sides in local coordinates, the relationships among local and global variables are given by [11]: With the above transformation, the governing equations for the seabed model in the local axes (x ,z ) should be transferred to the global axes (x, z) as: in which p s is the pore-pressure in the slop seabed and u s and w s are the soil displacements in x -and z -directions, respectively. For the seabed model, the following boundary conditions are specified. At the seabed surface, the pore pressure (p s ) should be equal to the wave pressure (p w ) and the stresses vanish [37], At the seabed bottom and two lateral sides, no vertical flow and displacement occur as: where n is the normal direction of the boundaries. In this study, the pipeline is considered as a rigid impermeable object. Then, the noflow boundary condition is applied to its surface: Based on the seabed model, the pore pressure and soil displacement can be determined. Then, the pore pressure gradient and two-dimensional seepage can be further obtained. In this study, both flow and seabed models are established in the OpenFOAM environment. The seabed model established by the QS model ignores the relative displacement acceleration between the fluid and soil particles. The communication of the two sub-models takes place at the boundaries in one direction. Therefore, this is a one-way integrated model. The flow chart for the one-way integration of two models is illustrated in Figure 2.
Solve Navier-Stokes and Volume of Fluid equations Export velocity field ( ), volume fraction function field (α 1 ), and dynamic wave pressure ( ) Impose dynamic wave pressure to the wave-seabed interface using specific boundary condition Solve Biot's equation using the segregated solution procdure Export pore-water pressure ( ), and soil displacement field ( , , ) Output and save results according to pre-set time interval

Grid Independence Tests
In general, the density of grid tremendously affects the accuracy of numerical results. In order to ensure the accuracy of numerical simulations, different grid densities of wave and seabed models were adopted to perform the validation of grid independence. Figure 3a,b illustrate the variations of the wave free surface versus time with different grid sizes in the xand z-directions, respectively. The input data are tabulated in Table 1. It is clear that the maximal grid size is ∆x = L/80 in the x-direction and ∆z = H/40 in the z-direction, when the numerical results are not affected by the grid size. Therefore, the grid size is designed as ∆x = L/80 and ∆z = H/40 for the flow model.  Similarly, for the seabed model, the gradient of pore pressure in the z-direction (∂p s /∂z) varying time with different grid size is indicated in Figure 4. As shown in the figure, ∆z = h s /200 can satisfy the requirements of accuracy. Therefore, it is used in the present numerical seabed model.

Model Validations
Since the study on the impact of wave(current)-induced seepage to the sediment motion is limited, no experimental data have been available in the literature. In the section, we can only validate the wave-induced pore pressure in a seabed for the proposed model. Both flow and seabed models will be validated by comparing the numerical results with previous data from other numerical solutions and experimental data. The following validations will be presented in this section.

•
Validation of the flow model: Comparison with the experimental data for the wave models with a pipeline laid on a seabed [38]. • Validation of seabed model #1: Comparison with the laboratory experimental data with a fully buried pipeline in a seabed [28] and the boundary element model [29]. • Validation of seabed model #2: Comparison with the wave flume tests for the seabed model with a trenched layer [31].
In addition to the above validations, numerous validation of the present numerical model, including flow and seabed sub-models, for the cases without a structure are available in the authors' previous work [34].

Validation of the Flow Model
The validation of the flow model is to compare the flow velocity above the pipeline calculated by the numerical model with the experimental results of Mattioli et al. [38]. In this comparison, the present model is simplified to the case without a seabed, i.e., wavepipe interactions. This is to ensure that the sketch of the flow model is the same as that of Mattioli et al. [38]. The input data of the validation are listed in Table 2. Table 2. Input data used in the comparison with Mattioli et al. [38].

Parameters Value
Wave

Validations of the Seabed Model
In the section, the seabed model will be validated by comparing it with previous data from the literature. Among these, in the first validation for the seabed model, the seabed model with a fully buried pipeline is compared with the experimental data of Turcotte et al. [28] and the previous numerical model of Cheng and Liu [29]. In the experiments, the trench layer is surrounded by impermeable walls. The input data of this comparison are given in Table 3, and the distributions of the pore pressure (|p s |/p 0 ) along the pipeline (θ p ) under the wave loading are represented in Figure 6. As shown in the figure, the present model can capture the essential trends of the previous numerical model [29] and the experimental data [28]. Table 3. Input data used in the validation with Turcotte et al. [28] and Cheng and Liu [29].

Parameters Value Parameters Value
The second validation of the seabed model is for the case of a half buried pipeline in a trenched layer [31]. In their experiments, a sediment basin with the size of 2 × 1 × 0.58 m (L × W × D) was used to explore the seabed response around the pipeline. The results of case 49 were used to compare with the present model. The detailed input data were tabulated in Table 4.

Parameters Value Parameters Value
Wave height (H) 0.12 (m) Wave theory (T) 1.6 (s) The comparisons of the normalized excess pore pressure (| u e | /σ 0 ) around a pipeline are presented in Figure 7a, while the distribution of the pore-water pressure (| u e |) versus seabed depth (z/h) is plotted in Figure 7b. Here, σ 0 = (γ s − γ w )z(1 + 2K 0 )/3, where K 0 is the coefficient of lateral earth pressure at rest. As seen from Figure 7, the trend of the numerical solution is consistent with the experimental data [31]. Comparisons of the excess pore pressures (|u e |/σ 0 ) between the the experimental data [31] and the present model.
In the above validations, both flow and seabed models in different conditions were validated. This demonstrates the capacity of the present model for flow and seabed fields, which would further calculate the seepage force correctly.

New Shields Number with Seepage Flow
In general, the onset of sediment motion will occur when the driving force acting on the particles is greater than the resistance. Based on the force balance of a soil particle, Shields [6] defined the Shields number (θ) as where γ w and γ s are the unit weight of fluid and sediment, respectively; d 50 is the mean grain size and τ b is the bed shear stress, which is determined by: in which u * is the friction velocity at the seabed surface, and ρ w is the density of the fluid.
In this study, a non-cohesive sandy seabed is considered. Considering a coordinate system as shown in Figure 8, the soil particles at the seabed surface are subjected to the following forces: (1) In the z-direction: they are submerged weight (W ), lift force (F L ), and seepage force in the z-direction (F sz ). (2) In the x-direction: they are drag force (F D ) and seepage force in the x-direction (F sx ). The equilibrium equation of the sediment particles at the seabed surface can be expressed as in which f is the static friction coefficient, and F R is the resistant force in the x-direction.
Herein, all soil particles are assumed to be spheres with a uniform diameter of d 50 , then the submerged weight, lift force, and drag force in (14) can be written as: where C L and C D are the lift and drag coefficients, respectively, and u 0 is the flow velocity approaching to the particle at the seabed surface.
Substituting (15)- (17) into (14), the following relation can be obtained: in which τ c is the bed shear stress at the critical condition, R e is the Reynolds number; f 2 (R e )(= u 0 /u * ) is the ratio of velocity of the particles and friction velocity [5].
The left side of (18) is the new Shields number modified of two-dimensional seepage. Compared with (12), the extra part of reflects the impact of two-dimensional seepage. In order to calculate the extra part, seepage force should be calculated firstly, which would be obtained above.
Considering the soil-fluid mixture of the seabed, the seepage within a porous medium can be decomposed into the xand z-directions and expressed as [34,39,40], in which n s is the soil porosity and p s is the pore-pressure in the seabed. It is clear from the right-hand-side of (14) that the resultant force in the xdirection is the same as the horizontal seepage force (F sx ) and drag force (F D ). Therefore, the sign must be consistent when they have the same direction and vice versa. However, the sign of F D is always positive as seen from (17), and its direction is the same as the flow. For that, u 0 |u 0 | is introduced to adjust the sign of F sx to realize that F sx is positive when it has the same direction with F D and F sx is opposite when it has the opposite direction with F D . However, u 0 |u 0 | can be used only when u 0 is not equal to zero and the condition when u 0 = 0 should be discussed separately, which is due to the fact that u 0 |u 0 | is meaningless as the denominator and can not be a parameter to adjust the sign of F sx . For the special condition, drag force is zero, and there is only F sx in the horizontal direction. Hence, only |F sx | should be considered. Now, the horizontal seepage force (F sx ) is rearranged as Here, the scaling factors e i (i = x, z, xz) are introduced to represent the impact of seepage to the modified Shields number: Now, we introduce the modified Shields numbers, θ z and θ xz , to represent the impact of seepage: in which the conventional Shields number(θ) given in (12) and the modified Shields number with vertical seepage given in (25) [24]. Note that d 50 was used to replace the particle size of soil and uniform soil was considered in the above derivation of the modified Shields number. That is, we ignored the influence of wide graded of soil on the modified Shields number in practical situations. This is an important point to further investigate the influence of the widely grade of soil on the modified Shields number in the future. Now, we further consider the influence of seepage on a sloping bed. In a sloping seabed, a coordinate system has a parallel to the slope as the x -direction and perpendicular to the slope as the z -direction, as shown in Figure 9. The force balance on a sloping bed is presented in the figure. Based on the force balance, the equilibrium equation in the x -direction can be written as: Substituting (15)- (17) into (27), the following relation can be obtained: With the similar process, the seepage components in the x and z can be expressed as: Herein, the scaling factors for a sloping seabed (e x and e z ) are introduced below: Similar to the case with a flat seabed, the scaling factor for a sloping seabed, e x z = e x + e z , is introduced to represent the influence of 2D seepage on the Shields number. Meanwhile, similar to (25) and (26), modified Shields numbers θ z and θ x z on a sloping seabed can be written as: Note that, when the angle of slop seabed is zero, i.e., a flat seabed, the results of modified Shields numbers θ z and θ x z are identical to θ z and θ xz in (25) and (26).
It is clear that the influence of the 2D seepage on sediment incipient motion depends on the scaling factor (e xz or θ x z ). The ratio, 1/(1 − e xz ) or 1/(1 − e x z ), is a magnification of the Shields number under the influence of 2D seepage in a flat seabed or sloping seabed. According to this, it is found that 2D seepage would promote the sediment incipient motion when 1/(1 − e x z ) > 1 and restrain the onset of sediment motion when 0 < 1/(1 − e x z ) < 1. In addition, for the special case when 1/(1 − e x z ) < 0, the seepage force is larger than the submerged weight force and the soil particle is in the suspension state, which is beyond the scope of our research.
Equations (25), (26), and (32) used for the evaluation of the influence of 2D seepage on sediment incipient motion could further affect the evaluation of the process of scour and erosion of sediment by gravity currents under non negligible slopes, which are extremely significant in the design of marine structure. In practical offshore engineering, the scour protections are an extremely important part to ensure the stability of the superstructure [41]. Hence, more accurate predictions of sediment initiation motion and scour can make the scour protections design more effective. It is obvious that the seepage promotes the onset of sediment when 1/(1 − e x z ) > 1, which means that the protective layer of scour design should be large enough. In addition, seepage restrains the onset of sediment, and the protective layer of scour design can be smaller when 0 < 1/(1 − e x z ) < 1.
Recently, two studies considered the seepage in the sediment motion or local scour. First, Guo et al. [24] re-derived the Shields number with the vertical seepage. Second, Li et al. [16] further considered the vertical seepage in their CFD scour model for a sloping seabed. However, their studies ignored the multidimensional of seepage. For that, the influence of horizontal is added in the paper. It is noted that the definition of F sz used in Li et al. [16] was only considered the magnitude of the upward seepage. Since they used a 1D approach, they only focused on the case under wave troughs, where the seepage (F sz ) is always upward. In this study, the whole cycle is considered, and the results are displayed below.

Impact of 2D Seepage on Sediment Incipient Motion
In offshore environments, pipelines are particularly important for the transportation of oil and gas. Numerous studies for the wave-seabed-pipeline interaction have been reported in the literature in the past [42,43]. This includes the numerical [26,[44][45][46] and experimental investigations [31,47,48] for the wave/current-induced seabed liquefaction around a pipeline. Here, the proposed numerical model of wave-seabed interaction around a pipeline (see Section 2) is employed to examine the impact of 2D seepage on the sediment incipient motion.
The scaling factors (e z and e xz for a flat seabed; e z and e xz for a sloping seabed) reflect the effects of 2D seepage on the sediment motion. The impact of seepage on the sediment incipient motion around a pipeline will be analyzed first by the present numerical model. Then, a parametric analysis will be adopted to explore the influence of horizontal seepage. The input data of the numerical examples are listed in Table 5. The applicable ranges of various wave theories are marked in Figure 10.

Effect of Seepage on Shields Number
In a flat seabed, the scaling factors (e x , e z and e xz ) play extremely important roles in evaluating and researching the influence of seepage on the sediment incipient motion. Variations of the wave-induced scaling factors without a structure are plotted in Figure 11. It is clear from the figure that the positive e z would occur around the wave trough and promote the onset of sediment motion. On the contrary, negative e z would occur around the wave crest and restrain the onset of sediment motion. However, e x would have a discontinuity that occurs at the position where pore-pressure is zero, in which the flow velocity is also equal to zero. The special condition corresponds to the theoretical results when u 0 = 0 in (22), and the value of the position is marked by * in the figure. Furthermore, the discontinuity is extended in e xz . The pipeline would impact the hydraulic environment and further impact the soil response and seepage below the pipeline. To explore the influence of pipeline on seepage, three different scaling factors (e x , e z and e xz ) for the cases with and without pipeline in different times are illustrated in Figure 12. The pipeline is located in the position where x/L = 0. As seen from the figure, the value of e x changes obviously under the pipeline. The discontinuity occurs in e x , and further occurs in e xz , which are caused by the special case with zero flow velocity, and the values of the special positions are represented by a special symbol ( * ). When the discontinuity appears below the pipeline, the discontinuous position is different compared with the case without a pipeline referring to Figure 12b,d, which is caused by the different flow velocity with and without the resistance of the pipeline. On the contrary, when the discontinuity appears far from the pipeline, the discontinuous positions are consistent with the cases without a pipeline due to the fact that the influence of pipeline on flow velocity is negligible in the position. The trend of scaling factor e z is opposite before and behind the pipeline compared with the case without a pipeline and the extent of the change is smaller than e x . Hence, the variety of e xz below the pipeline mainly depends on the e x . Therefore, the horizontal seepage is essential and could not be ignored, especially around the pipeline. The previous studies with vertical seepage [16,24] need to be re-visited.
Considering that the scaling factors (e x , e z and e xz ) vary at different times, as illustrated in Figure 12, the corresponding envelop curves for scaling factors in a wave cycle are represented in Figure 13. It is obvious that the scaling factors e x , e z , and e xz change in the range from the minimum to the maximum. The scaling factors e x and e xz have an obvious increase below the pipeline, but the change of e z is smaller, which is consistent with Figure 12.   The e x has a slight increase behind the pipeline because of the lee-wake vortex. Meanwhile, e x has an obvious increase below the pipeline, which is caused by the resistant force of the pipeline on flow water. The results of θ xz /θ with and without pipeline are plotted in Figure 14b. As seen from the figure, the trends of θ xz /θ are similar to the scaling factors. Furthermore, the maximum of θ xz /θ can reach 3.2, which is large enough to impact the sediment incipient motion, and further processing of scour around the pipeline.

Effects of Pipeline Characteristics
The diameter of pipeline is one of the key parameters to affect the onset of sediment motion around a pipeline, which impacts the degree of hindrance on the flow water. Therefore, four different pipeline diameters (D = 1, 2, 3, 4 m) were adopted in the numerical simulation. The θ xz /θ obtained from the corresponding envelop curves of e xz are plotted in Figure 15a. In the cases, the flow gap ratio (h/D) is prescribed as 0.1. Hence, the gap from the bottom of pipeline to the seabed surface h increases with the pipeline diameter, which causes the value of θ xz /θ below the pipeline to decrease. However, θ xz /θ reaches a peak before the pipeline and a decrease behind the pipeline, which increases with the pipeline diameter.   Figure 15b. It is observed that θ xz /θ increases with the reduced flow gap ratio. That is, the effect of seepage on sediment incipient motion below the pipeline would increase with the reduced flow gap ratio. This is due to the enhancement of disturbed fluid below the pipeline when the pipeline is near the seabed surface. However, the position where θ xz /θ starts increasing before the pipeline is not affected by the flow gap ratio (as seen from Figure 15a) and pipeline diameter (as seen from Figure 15b). For the case of a partial buried pipeline, i.e., e = −0.5, a small peak occurs behind the pipeline that is caused by the lee-wake vortex.

Effects of Wave Characteristics
It has been reported that wave parameters will affect the seepage. In this section, two wave parameters, dimensionless wave height (H/gT 2 ), and water depth (d/gT 2 ) would be considered in the parametric study. As shown in Figure 16a, the ratio of modified Shields number and conventional Shields number (θ xz /θ) obtained from the corresponding envelop curves of e xz increase with the H/gT 2 . However, θ xz /θ increases below the pipeline, and the amplification is larger with the increasing H/gT 2 . When the H/gT 2 is small enough, for example, the case of H/gT 2 = 0.00319, the existence of pipeline hardly affects sediment incipient motion below the pipeline. In addition, the increase of θ xz /θ before the pipeline is larger with the increasing H/gT 2 , which is due to the larger turbulence before the pipeline.  Water depth is another key wave parameter to affect the pressure at the bottom of flow field. To investigate the relationship between water depth and the modified Shields number, different dimensionless water depths (d/gT 2 = 0.05079, 0.06371, 0.07645, 0.08919, 0.10194) are adopted, which are in the range of regular waves as shown in Figure 10. The results are plotted in Figure 16b. It is found that θ xz /θ obtained from the corresponding envelop curves of e xz decreases with the increasing of d/gT 2 . However, θ xz /θ is closer to 1 with the larger dimensionless water depth (d/gT 2 ). That is, the effect of seepage on sediment incipient motion is quite small in deep water. The influence of pipeline on θ xz /θ increases with the reduced d/gT 2 , which is characterized as the increasing amplification below the pipeline.

Effects of Uniform Current (u)
As discussed previously, the authors' previous work [34] could be the first attempt to consider 2D seepage in the definition of Shield number. However, neither currents nor the pipeline were considered in their work. In this section, we further investigate the effects of current on the new Shield number. The change of current velocity will tremendously affect the dynamic pressure at the seabed surface and the soil response within the seabed, which would further affect the seepage field. The distribution of flow velocity near the pipeline is different for different combinations of wave and current. Figure 17 illustrates the distribution of velocity around the pipeline with different current velocities (u = −0.4 m/s, 0 m/s and 0.4 m/s). Note that the negative sign denotes that current and wave have opposite directions. As seen from Figure 17b, the flow velocity would decrease firstly and then increase above the pipeline by the action of wave only. However, a similar trend occurs at the upper left of the pipeline when wave and current have an opposite direction, and at the upper right of the pipeline when wave and current have the same direction, as seen from Figure 17a,c. Meanwhile, flow velocity forms a ring before and behind the pipeline when the direction of wave and current are opposite and identical, which are caused by the lee-wake vortex and lead to the change of seepage before and behind the pipeline. However, the lee-wake vortex occurring in the case when wave and current have opposite directions is smaller than the condition when wave and current have the same direction.
Herein, different current velocities u = −0.6, −0.4, −0.2, 0.0, 0.2 and 0.4 m/s are adopted to find the influence of wave and current co-induced seepage on sediment incipient motion. The compared results are plotted in Figure 18a. It is clear that θ xz /θ obtained from the corresponding envelop curves of e xz increase with the current velocity, integrally. However, a different trend is observed below the pipeline. A peak is formed, and the amplitude of θ xz /θ reduces with the increased current velocity below the pipeline when the direction of wave and current are opposite. By contrast, a valley is formed below the pipeline and the θ xz /θ increases with the current velocity, when wave and current have the same direction.
To further clarify the reasons for the formation of the peak and valley, the comparison of scaling factors (e z and e xz ) in different combinations of wave and current is illustrated in Figure 18b. It is known that the peak and valley mentioned above are due to the change of vertical seepage below the pipeline and the shake of e xz and θ xz /θ are caused by the horizontal seepage.

Effects of Seabed Characteristics
Shear modulus (G) is one of the key seabed properties affecting the seepage in the seabed [50]. In this section, the shear modulus is varied from G = 5 × 10 5 N/m 2 to G = 5 × 10 7 N/m 2 , which practically include all kinds of sand and silt. The numerical results are represented in Figure 19a, which shows the influence of shear modulus on θ xz /θ. It is obvious that the increase of shear modulus leads to an increased θ xz /θ, which means that the sediment incipient motion would occur more easily with the larger shear modulus. However, the increase of θ xz /θ is not obvious when the shear model is larger than 5 × 10 6 N/m 2 . The θ xz /θ increases from approximately 2.2 (at x = −2.5D) to 3.2 at the peak near x = 0, which represents the dramatic effect of seepage on sediment incipient motion, especially near the pipeline (x/D = −0.5 to 0.5). To clarify the influence of soil permeability in the z-direction (k z ), different permeability, k z = 10 −2 , 10 −3 , 10 −4 , 10 −5 , 10 −6 , and 10 −7 m/s are selected in numerical simulation. The influence of k z on θ xz /θ is illustrated in Figure 19b. As shown in the figure, the decrease of k z leads to the increase of θ xz /θ, and the peak value below the pipeline increases significantly. At the peak, close to x/D = 0, the maximum of θ xz /θ increases from about 1.15 to 3.7, which means that the conventional Shields number increases from 15% to 270% by considering two-dimensional seepage. This has demonstrated the significant effect of seepage on the Shields number.
Different degrees of saturation (S r = 0.9, 0.92, 0.94, 0.96, 0.98 and 1) are adopted to explore the influence of saturation on θ xz /θ, as illustrated in Figure 19c. Similar to the effect of k z , the θ xz /θ increases with the decreased degree of saturation. The increase before and behind the pipeline becomes more obvious with the decreased saturation, which are caused by the back-water and lee-wake vortex, respectively. Meanwhile, the point where θ xz /θ begins to increase before the pipeline gets closer to the pipeline with the increasing of saturation. At the peak near x/D = 0, the maximum of θ xz /θ increases from about 1.7 to 4.2, which indicates that the seepage tremendously impacts the sediment incipient motion in an unsaturated seabed.

Effects of Slope Angle (β)
In the offshore environment, the seabed surface may not always be flat. Therefore, the influence of slope angle (β) on sediment incipient motion needs to be investigated. The modified Shield number, θ x z , has been defined in (28). In this section, different slope angles (β = 0, π 72 , π 36 , π 24 and π 18 ) are used to investigate the influence of slope angle. The results of θ x z /θ are displayed in Figure 20. It is found that θ x z /θ would increase with the slope angle of seabed. However, a slope of θ x z /θ would occur before the pipeline and increase with the slope angle of seabed. In addition, the slope angle of θ x z /θ is closer to the slope angle of seabed. In addition, a valley would occur before the pile and the amplitude of the valley would increase with the slope angle.

Engineering Application: A Buried Pipeline in a Trenched Pipeline
Nowadays, the trenched layer has been commonly used in marine oil and gas projects in order to prevent the damage of the pipeline. In the study, the case of a partially buried pipeline in a trenched layer is further considered to simulate the practical engineering application. Although this problem has been studied recently [26,46,51], they mainly focused on the wave (current)-induced soil response and liquefaction near the pipeline but ignores the effect of seepage on the sediment motion. Hence, determining the influence of seepage on sediment motion around the half buried trenched pipeline has great engineering significance.
Herein, the influence of seepage on sediment motion around the trenched pipeline will be examined by the proposed model. The sketch of the numerical model with a trenched pipeline is displayed in Figure 21, and detailed input data for the trench are listed in Table 6. The other soil parameters are the same as the cases with the spanning pipeline above.  (N/m 2 ) Degree of Saturation (S r ) 98 (%) Permeability (k s ) 2.382 × 10 −6 (m/s) The distributions of corresponding enveloped scaling factors (e x , e z and e xz ) around the trenched pipeline are plotted in Figure 22. It is found that the e z decrease is approaching the pipeline. However, e x has a tremendous increase around the pipeline, which is due to the complex fluid environment and pore pressure. In addition, e x and e xz are larger before the pipeline than the position behind the pipeline. Meanwhile, e xz also increases around the pipeline, which is mainly dependent on the horizontal seepage flow.  To explore the different effects of seepage on sediment incipient motion between a trenched layer and a flat layer, Figure 23 compares θ z /θ and θ xz /θ around a trenched pipeline and a semi buried pipeline in a flat layer. It is observed from the figure that both θ xz /θ and θ z /θ in the trenched pipeline are larger than the pipeline in a flat layer, which are due to the faster flow velocity and larger pressure in the trench. However, the horizontal seepage has a greater contribution to the trenched pipeline than the pipeline in a flat layer by comparing the division of θ xz /θ and θ z /θ. The embedding methods of the trenched pipeline are multifarious. The full, half, and not buried trenched pipeline are widely used in the engineering environment, which corresponds to the case of e = −1, e = −0.5, and e = 0 in Figure 24. Note that the different ranges of the x-axis for the three cases are due to the different buried depths and the different lengths of the trench bottom. The full buried trenched pipeline hardly changes around the pipeline, which is due to the unaffected flow environment caused by the pipeline. However, the influence of seepage on the other cases becomes more significant around the pipeline.

Conclusions
In this paper, we derived a new Shields number with 2D seepage and discussed its effects on sediment incipient motion around the pipeline. The key findings of numerical examples are summarized here.
(1) The seepage presents periodic change along the seabed surface and tremendously changes around the pipeline. Meanwhile, the horizontal seepage has a significant contribution below the pipeline and should not always be ignored. (2) A discontinuity of the horizontal seepage is observed at the position where the flow velocity is equal to zero, and the absolute value is selected in the position. The discontinuity also extended in the scaling factors e x and e xz . (3) A peak of the modified Shields number is formed below the pipeline in most cases.
However, a valley formed below the pipeline only when the wave and current have the same direction, which is due to the changed fluid pattern around the pipeline.
(4) The effects of seepage become significant with large wave, shallow water, and large current in a seabed with lower saturation but larger shear modulus and permeability. However, the pipeline parameters only impact the seepage force near the pipeline. (5) The impact of seepage on sediment incipient motion increases clearly around the trenched pipeline, and the modified effect of Shields number is obviously better than the spanning pipeline.
In this paper, we attempt to connect the soil response with the Shield number and the onset of sediment motion. The future study will be further linked with a local scour model to simulate the scour process around the offshore structure. Furthermore, the gravity currents have a significant impact on the entrainment, transport, and deposition of sediments in a sloping seabed, as reported in Zordan et al. [52]. The proposed model requires further modification to include the effects of gravity currents to examine its effects with seepage and associated local scour, sinking, and floatation of a pipeline in the future.