Coupled Flow-Seepage-Elastoplastic Modeling for Competition Mechanism between Lateral Instability and Tunnel Erosion of a Submarine Pipeline

: The instability of a partially embedded pipeline under ocean currents involves complex ﬂuid–pipe–soil interactions, which may induce two typical instability modes; i.e., the lateral instability of the pipe and the tunnel erosion of the underlying soil. In previous studies, such two instability modes were widely investigated, but separately. To reveal the competition mechanism between the lateral instability and the tunnel erosion, a coupled ﬂow-seepage-elastoplastic modeling approach was proposed that could realize the synchronous simulation of the pipe hydrodynamics, the seepage ﬂow, and elastoplastic behavior of the seabed soil beneath the pipe. The coupling algorithm was provided for ﬂow-seepage-elastoplastic simulations. The proposed model was veriﬁed through experimental and numerical results. Based on the instability criteria for the lateral instability and tunnel erosion, the two instability modes and their corresponding critical ﬂow velocities could be determined. The instability envelope for the ﬂow–pipe–soil interaction was established eventually, and could be described by three parameters; i.e., the critical ﬂow velocity ( U cr ), the embedment-to-diameter ratio ( e/D ), and the non-dimensional submerged weight of the pipe ( G ). There existed a transition line on the envelope when switching from one instability mode to the other. If the ﬂow velocity of ocean currents gets beyond the instability envelope, either tunnel erosion or lateral instability could be triggered. With increasing e/D or concurrently decreasing G , the lateral instability was more prone to being triggered than the tunnel erosion. The present analyses may provide a physical insight into the dual-mode competition mechanism for the current-induced instability of submarine pipelines.


Introduction
Offshore exploitations of gas and oil have been turning from shallow waters to deep waters. By the end of 2018, about 14 large and medium-sized deep-water oil and gas fields had been discovered in the northern South China Sea; e.g., the Liwan 3-1 and the Lingshui 17-2 gas fields [1]. The submarine pipeline is an effective tool for the transport of oil and gas in the deep sea, and thus is termed as the lifeline of ocean resources. According to statistics, approximately 51 accidents occurred involving 315 submarine pipelines of CNOOC from 1986 to 2016 [2]. The structural failure of submarine pipelines usually involves complex flow-pipe-soil interactions. The instability of a submarine pipeline incorporates the lateral instability of the pipe, the local scour and seepage failure of soil, etc. (see Fredsøe [3]; Gao [4]; Drumond et al. [5]). As is well known, with the increase of the water depth, the surface-wave-induced water particle oscillations near the seabed are gradually weakened and finally vanish in the water approximately deeper than half of one wavelength, while ocean currents always exist even in the deep waters (see Shanmugam [6]). As deep-water pipelines are often laid directly on the seabed, an appropriate assessment of the instability of the shallowly embedded pipeline is particularly vital for the engineering design and safe operation of deep-water pipelines.
Under the action of ocean currents, the flow over a partially embedded pipeline and the seepage flow within the underlying soil can be generated concurrently, as illustrated in Figure 1. Thus, the on-bottom stability of a submarine pipeline involves complex flow-pipesoil interactions [7,8], but not only limits to the pipe-soil interaction [9]. When the lateral soil resistance provided by the seabed soil is insufficient to balance the hydrodynamic drag force on the pipeline, the pipeline will breakout laterally from its original location; i.e., lateral instability takes place (see DNV [10]). Meanwhile, if the seepage failure is triggered within the underlying soil due to the pressure drop from the upstream to the downstream of the pipe, tunnel erosion beneath the pipe can also be induced, and further result in the pipeline spanning (see DNVGL [11]). J. Mar. Sci. Eng. 2021, 9, x FOR PEER REVIEW 2 of 25 known, with the increase of the water depth, the surface-wave-induced water particle oscillations near the seabed are gradually weakened and finally vanish in the water approximately deeper than half of one wavelength, while ocean currents always exist even in the deep waters (see Shanmugam [6]). As deep-water pipelines are often laid directly on the seabed, an appropriate assessment of the instability of the shallowly embedded pipeline is particularly vital for the engineering design and safe operation of deep-water pipelines.
Under the action of ocean currents, the flow over a partially embedded pipeline and the seepage flow within the underlying soil can be generated concurrently, as illustrated in Figure 1. Thus, the on-bottom stability of a submarine pipeline involves complex flow-pipe-soil interactions [7,8], but not only limits to the pipe-soil interaction [9]. When the lateral soil resistance provided by the seabed soil is insufficient to balance the hydrodynamic drag force on the pipeline, the pipeline will breakout laterally from its original location; i.e., lateral instability takes place (see DNV [10]). Meanwhile, if the seepage failure is triggered within the underlying soil due to the pressure drop from the upstream to the downstream of the pipe, tunnel erosion beneath the pipe can also be induced, and further result in the pipeline spanning (see DNVGL [11]). Since the 1970s, the lateral instability of pipe and the tunnel erosion of soil have been investigated separately from either geotechnical or hydrodynamical perspectives. An accurate estimation of soil lateral resistance is crucial for the design of the on-bottom stability of submarine pipelines [10]. Early studies regarding the lateral instability of pipe were mainly carried out through mechanical-loading model tests, and several empirical pipe-soil interaction models were proposed based on the results of mechanical-actuator tests (e.g., Wagner et al. [9]; Verley and Sotberg [12]; Verley and Lund [13]). Some models are currently encapsulated in recent codes; e.g., DNV-RP-F109 and DNV GL-RP-F114. It should be noticed that only the hydrodynamic forces on the pipe (including the drag and the lift forces) were simulated with a mechanical actuator mounted onto the test pipe, but the influence of the waves and (or) the current on the seabed was ignored. Since the 1970s, the lateral instability of pipe and the tunnel erosion of soil have been investigated separately from either geotechnical or hydrodynamical perspectives. An accurate estimation of soil lateral resistance is crucial for the design of the on-bottom stability of submarine pipelines [10]. Early studies regarding the lateral instability of pipe were mainly carried out through mechanical-loading model tests, and several empirical pipe-soil interaction models were proposed based on the results of mechanical-actuator tests (e.g., Wagner et al. [9]; Verley and Sotberg [12]; Verley and Lund [13]). Some models are currently encapsulated in recent codes; e.g., DNV-RP-F109 and DNV GL-RP-F114. It should be noticed that only the hydrodynamic forces on the pipe (including the drag and the lift forces) were simulated with a mechanical actuator mounted onto the test pipe, but the influence of the waves and (or) the current on the seabed was ignored.
It can be recognized that the lateral instability of a shallowly embedded pipeline involves not only the pipe-soil interaction, but also flow-pipe-soil interactions. Based on similarity analyses and flume experiments under waves or steady flow conditions, three characteristic stages were observed during the process of a pipeline losing stability: (I) onset of local scour; (II) pipe rocking; and (III) pipe breakout [7,14]. Moreover, it was found that the Froude number is the non-dimensional controlling parameter for the pipe's lateral instability. Then, an empirical relationship between the Froude number and the dimensionless submerged weight of the pipe was established for predicting the lateral instability of very shallowly embedded pipelines (e.g., e/D < 0.05, in which e/D is the embedment-to-diameter ratio) [8,14]. Large-scale O-tube tests were recently conducted to further investigate the local-scour effects on the on-bottom stability of submarine pipelines [15,16]. The lateral-stability mechanism neglecting local scour may not well reflect the actual behavior of flow-pipe-soil interactions [3,16]. Nevertheless, a quantitative characterization of the tripartite flow-pipe-soil interactions is still lacking in the existing design codes for the on-bottom stability design of submarine pipelines [10].
In addition to the aforementioned lateral instability, tunnel erosion beneath a partially embedded pipeline can also be triggered under the action of hydrodynamic loads. The triggering mechanism for the tunnel erosion was attributed to the development of local scour of sand particles at the seabed surface around the pipe. Nevertheless, flume experiments showed that when an impermeable plate was laid at the upstream of the test pipe in a steady flow, the tunnel erosion was effectively suppressed [17], indicating that the hydraulic gradient played a key role in the occurrence of tunnel erosion [18]. A series of flume observations [19] indicated that during the process of the pipe being suspended, there usually existed three characteristic stages; i.e., stage I: local scour around pipe; stage II: tunnel erosion triggered by seepage failure; and stage III: complete suspension of the pipe. While the local scour always emerged during the pipe suspension process, its development was observed much slower than the occurrence of tunnel erosion [20]. It has been well recognized that tunnel erosion is essentially due to the seepage failure of the underlying soil resulting from a pressure drop, rather than the progressive development of local scour, even for an uneven seabed [21]. Note that the model pipe was fixed in the previous physical modeling of tunnel erosion, so that the lateral instability of the pipe could not be involved.
In actual submarine geological and hydrodynamic environments, multi-mechanical processes could be involved in the flow-pipe-soil interactions, as illustrated in Figure 1, including the shear flow passing over the pipe, the seepage flow within the seabed, and the elastoplastic soil deformation around the pipe, etc. Previous studies mainly concentrated on the individual aspects of the flow-pipe-soil interactions, such as the pipe-soil interaction, the flow-pipe interaction, and the flow-soil interaction [3,4,22]. It was recently found that the lateral instability of pipe and the tunnel erosion of soil are closely correlated with each other [23]. That is, with increasing flow velocity, both of these two instability modes could take place, but only the weaker one (i.e., the mode requiring lower critical flow velocity) would emerge eventually. As an uncoupled approach, such correlation analysis [23] was made on the basis of the existing critical relationship for the lateral instability of pipe [24] and that for the tunnel erosion of soil [20]. Coupled modeling is needed for simulating the multi-mechanical processes in such flow-pipe-soil interactions. A poro-elastoplastic model was recently developed to simulate the buried pipe-soil interaction under wave loading, and indicated that the upward seepage flow under wave troughs may significantly weaken the effective stress of the soil and further reduce the uplift resistance to the buried pipe [25]. Similarly, in the present study, the ocean current may induce seepage flow in the soil beneath the pipe, which would have effects on both the lateral instability and the tunnel erosion of the partially embedded pipe. For the case of the buried pipeline, the wave pressure on the seabed mudline was generally analytically treated as a preset load boundary condition in the existing numerical models [25,26]; nevertheless, for the case of a partially embedded pipeline in the present study, the flow-induced pressure drop on the mudline near the pipe and the corresponding hydrodynamic loads on the pipe could not be predicted analytically. This motivated us to design a coupled flow-seepageelastoplastic simulation to reveal the competition mechanism for the aforementioned two instability modes.
In the present study, a coupled flow-seepage-elastoplastic modeling approach was proposed that is capable of implementing the sequential simulation of the flow field over a pipe, the seepage-flow field, and the elastoplastic stress-strain behavior of the underlying soil. The numerical model was firstly validated with previous experimental results. The competition mechanism between the lateral instability and tunnel erosion was then investigated for a partially embedded pipeline under the action of an ocean current. Eventually, the instability envelope for the flow-pipe-soil interaction system was established. The flow-pipe-soil interaction for a partially embedded pipeline under the action of ocean currents can be assumed as a plane strain problem. To simulate the incompressible viscous flow over the pipe laid on a porous seabed, the Reynolds-averaged Navier-Stokes (RANS) equations and the continuity equation were employed; these can be expressed in the two-dimensional (2D) Cartesian coordinate system as follows:

A Coupled
where u i , u j are the averaged velocities of fluid; u i and u j are the fluctuating velocities; x i (or x j ) are the coordinates in the horizontal and vertical directions, respectively (see Figure 2); t is the time; ρ f is the mass density of fluid; p is the flow pressure; and ν f is the kinematic viscosity of fluid. For the steady flow, the term ∂u i /∂t in Equation (1) vanishes. The term of turbulent fluxes can be approximated by the Boussinesq assumption as: where δ ij is the Kronecker delta with 1 for i = j and 0 for i = j, k is the turbulent kinetic energy (i.e., k = u i u i /2), and ν t is the turbulent viscosity. A turbulence model is necessary to provide a value for the turbulent viscosity (ν t ) in Equation (3) (see Shih [27]; Durbin [28]). For the simulation of the turbulent flow around a free-spanning pipeline, several turbulence models were compared in [29], including the standard k-ε (high) and the low Reynolds number k-ω. It was indicated that the standard k-ε model could well predict the mean velocity, and the k-ω model with a no-slip boundary on the cylinder surface may provide a better prediction for the vortex shedding. For the instability of the integrated flow-pipe-soil system, hydrodynamics exerted on the pipeline and the flow pressure drop are key parameters. The relative error calculated with the k-ε and the k-ω model was found to be generally less than 15% through a series of trial tests. As such, in the present study, the standard k-ε turbulence model [30] was employed to calculate the turbulent viscosity coefficient of the flow-field for its advantageous of convergence rate and low memory requirements; i.e.,: , ν t is defined as ν t = C µ k 2 /ε, with ε denoting the turbulent energy dissipation rate; G k is defined as G k = −u i u j ∂u i /∂x j ; and the constants in the model are calibrated by comprehensive data fitting for a wide range of turbulent flows [31]: In this study, the seabed soil was assumed to be a saturated, homogeneous isotropic porous medium. It was assumed that the seepage flow in the seabed ob Darcy's law. According to Biot's consolidation theory [32], the continuity equation f isotropic porous soil is:

Seepage Flow within the Soil
In this study, the seabed soil was assumed to be a saturated, homogeneous, and isotropic porous medium. It was assumed that the seepage flow in the seabed obeyed Darcy's law. According to Biot's consolidation theory [32], the continuity equation for an isotropic porous soil is: where k s is the permeability coefficient of the soil; n is the porosity of the soil; K is the apparent bulk modulus of pore water; and ε ii is the volumetric strain of the soil; i.e., ε ii = u s i,i (i = 1, 2). In the simulation of the steady-current-induced pore pressure around the partially embedded pipe, the transient effects of volumetric strain (∂ε ii /∂t) and excess pore pressure (∂u w /∂t) can be further ignored, so Equation (6) can be simplified as the Laplace equation: ∇ 2 u w = 0.

Elastoplastic Behavior of the Soil
In the present study, to simulate the multi-mechanical processes of a partially embedded pipeline, the seepage flow and elastoplastic behaviors of the soil were considered concurrently.
The stress balance equation in the soil can be expressed as: where σ ij is the total stress tensor (positive in compression) and the subscripts i and j (i, j = 1, 2) indicate the horizontal and vertical directions; ρb i is the body force, b i is the body acceleration per unit mass (b 1 = 0, b 2 = g is the gravitational acceleration), and ρ is the mass densities of the soil; i.e., ρ = nρ f + (1 − n)ρ s , in which, ρ s and ρ f are the mass density of the soil particles and the pore water, respectively. Based on Terzaghi's effective stress principle, the total stress can be divided into two components: where σ ij and u w are the effective stress tensor and the excess pore pressure in the soil, respectively. Based on the deformation continuity condition, the total strain tensor can be written in terms of displacement gradients: where ε ij represents the strain tensor, and u s i denotes the soil-displacement component. The elastoplastic constitutive relation was adopted, and can be expressed in terms of infinitesimal increments (see Potts and Zdravkovic [33]): where D ep ijkl stands for the elements of the elastoplastic constitutive matrix; and dε kl is the incremental strains, which can be split into an elastic dε e kl and a plastic dε p kl component: Alternatively, the incremental stress dσ ij can also be computed by the incremental elastic strain dε e kl with the elastic constitutive matrix D e ijkl in the following form: (12) in which the Lame parameters are respectively. E s is the elastic modulus of soil, and ν s is the Poisson's ratio of soil. In Equation (12), the plastic strain can be calculated on the basis of yield function and the flow rule of the elastoplastic model: where Q p is the plastic potential. For a perfectly plastic material, when substituting Equations (12) and (13) into the consistency condition, the plastic multiplier dλ can be determined: where: In Equation (15), F y is the yield function. The yield surface encloses the elastic region defined by F y < 0, while the plastic flow occurs when F y ≥ 0 for an elastic trial stress. Numerous advanced constitutive models have been proposed during last decades [34][35][36][37]. For simplicity in the coupled flow-seepage-elastoplastic modeling, an elastic-perfectly plastic model was utilized to simulate the elastoplastic behavior of such a seabed soil under drained conditions. In the present model, the Drucker-Prager (D-P) yield criterion [38] was chosen, which can be expressed as: where J 2 is the second deviatoric stress invariant (the von Mises equivalent stress q = √ 3J 2 ), I 1 is the first invariant of Cauchy's stress tensor (the equivalent mean stress σ m = −I 1 /3), and α and K are the material parameters of the Drucker-Prager criterion, which can be matched to the coefficients in the Mohr-Coulomb criterion. For the case of 2D plane-strain applications, they can be expressed as: in which ϕ is the internal friction angle of soil, and c is the cohesion of soil. In such an elastic-perfectly plastic model, no hardening/softening law is assumed. The sand at shallow depth is under a low stress level, which normally shows dilation during shearing, with a high dilation angle possibly close to internal friction angle [39]. For simplicity, the associate flow rule implying the friction angle is equal to the dilation angle was assumed (Q p = F y ).

Geometry and Computational Meshes
As illustrated in Figure 2a, the geometry of the coupled finite element (FE) model mainly consisted of three domains; i.e., the water, the pipeline, and the soil. In the present simulation, a pipe was located on the seabed surface with an embedment (e), the water depth was set as 8D, and the soil depth was chosen as 10D, in which D is the outer diameter of pipe. The left and the right boundaries were 10D and 15D from the center of the pipe, respectively. Such sizes of the water and soil domains were proved to be sufficient to eliminate blockage and boundary effects (see Gao and Luo [20]).
The free triangular meshes were used in the water and the pipe domains, and the structured quadrilateral meshes were generated in the soil domain through a mapped meshing technique. Lagrange elements with order 2 were chosen for the entire domains, and more refined grids were used in proximity to the pipe to ensure the computation accuracy, as shown in Figure 2b. The mesh resolution and mesh element quality were key aspects for the coupled simulation. The mesh quality measures included skewness, maximum angle, volume versus circumradius, etc. The mesh quality could be described with a dimensionless parameter between 0 and 1, where 1 represents a perfectly regular element, and 0 represents a degenerated element [40]. In the present model, the number of all elements was approximately 3 × 10 4 , and the highest resolution was about 0.001 m. A high mesh quality was measured through element statistics (the minimum value of a series of average mesh quality values for different measured aspects was about 0.9), which proved sufficient for guaranteeing the convergence and efficiency of the results.

Boundary Conditions
The boundary conditions were set as follows (see Figure 2a): W1 (the inlet of the water domain): a constant undisturbed flow velocity u 1 | x=0 = U was specified. The inlet value for the turbulent kinetic energy and dissipation could be evaluated with k = 1.5(U I t ) 2 , ε = C µ 3/4 k 3/2 /L t , in which turbulent intensity I t = 0.05 and turbulence length scale W2 (the top of the water domain): a no-flow symmetry boundary was set at y = 18D. 3.
W3 (the outflow boundary): the pressure was given a reference value p = 0, whereas the other flow variables (velocity, turbulent kinetic energy, dissipation) were allowed to adjust freely with zero x-gradient conditions. The suppressed backflow was selected to prevent fluid from entering the domain through the boundary. 4.
WS4, WS5 (the water-seabed interface at the left and right of the pipe): the logarithmic wall function [42] was implemented; The effective normal stress and the shear stress vanished: σ ij i(j)=1,2 = 0; and the excess pore pressure was equal to the water pressure of the flow field at the surface of the seabed: u w y=−10D = p.

5.
WP6 (the pipe-water interface): no flow through the impermeable wall of the pipe (i.e., ∂p/∂n = 0); similar to WS4 and WS5, the logarithmic wall function was implemented at the pipe-water interface. 6.
S7, S9 (the left and the right lateral boundaries of the soil domain): no seepage was induced at the lateral boundaries: ∂u w /∂x = 0; and the normal displacement was set as zero: S8 (the bottom of the soil domain): an impermeable fixed bottom was set at y = 0; i.e., zero displacement: u s i | i=1,2 = 0; and no vertical seepage: ∂u w /∂y = 0. 8.
PS10 (the pipe-soil interface): the wall of the pipe was assumed to be impermeable, and there was no pore pressure gradient at the pipe surface: ∂u w /∂n = 0. Both the rolling and the sliding motions of the pipe were permitted for the freely laid pipe. 9.
A contact-pair algorithm was adopted to describe the pipe-soil interfacial behavior [40]. The pipe and seabed surface were assigned as the source and destination boundaries in the contact pair, respectively. Hard contact was chosen for the normal contact, and normal tensile stress was not allowed along the pipe-soil interface. When the pipe-soil contact surface was closed, it conveyed tangential stress. The Coulomb friction theory was used for the interfacial tangent contact. The pipe and the soil in the contact pair adhered to each other if the frictional shear stress (τ) was less than the critical one (τ crit ). Once τ crit was exceeded, slippage along the interface occurred. In the Coulomb friction model, the friction coefficient (µ F ) is defined as where ϕ µ is the friction angle of the pipe-soil interface. The value of ϕ µ generally ranges from 0.5 to 1.0 of the soil friction angle ϕ, depending on the interface characteristics and relative movement between the pipe and the soil [43]. For a smooth pipe on a sandy seabed, a value of ϕ µ = 0.58ϕ was herein considered, thus µ F = 0.37 was adopted in the present numerical simulations.

Properties of Materials
For the numerical simulations, the properties of the steel pipe, the saturated sandy seabed, and the steady flow are summarized in Table 1.

Coupling Algorithm
The aforementioned governing equations were solved using the software COMSOL Multiphysics with a coupling algorithm to obtain all the variables sequentially, as illustrated in Figure 3.

1.
Firstly, the responsibility of the flow field over the pipe (Equations (1)-(5)) in one computational step was to determine the pressure around the pipe, and to provide pressure/force acting on the seabed and structures as a surface boundary condition for the geotechnical simulations. As such, the continuity of the water-soil interface pressure (p-u w ) in every computational step was ensured, and the one-way coupling of flow and seepage fields could be implemented.

2.
Then, by solving the continuity equation (Equation (6)), the pore pressure could be obtained. Substituting the Terzaghi's effective stress principle (Equation (8)), the relationship between total strain and displacement (Equation (9)), and the elastoplastic constitutive model (Equation (10)) into the stress balance equation (Equation (7)), the governing equation for soil deformation with the unknown variables of soil displacement and pore water pressure could be solved through iterative computation.

3.
Concurrently, the soil's total strain, and its elastic and plastic components, could be calculated with the deformation continuity condition (Equation (9)) together with the yield function and the flow rule (Equations (13)- (16)). Based on the constitutive model (Equation (12)), the effective stress within the soil around the partially embedded pipe could be acquired.

4.
Based on the criteria for the lateral instability and tunnel erosion of the pipe (Equations (20) and (21) in Section 4.2), we could then examine whether the two instability modes were triggered.

Verification of Numerical Model
In this section, the proposed numerical model is validated by comparing with the existing experimental or numerical results. Both the pressure distribution at the mudline near the pipe and the critical velocity for the lateral instability of a pipe are examined, respectively.

Pressure Distribution at the Mudline near the Pipe
The proposed model was firstly validated against the existing experimental data regarding pressure distribution at the mudline near the pipe laid above or on the bed surface. The calculated distributions of pressure coefficient (Cp) at the mudline in the proximity of the pipe are plotted in Figure 4a, in which the experimental data of Bearman and Zdravkovich [44] are also provided. The pressure coefficient Cp is defined as follows: where Ps is the steady-flow-induced pressure at the mudline around the pipe, and r P is the reference pressure. In the experiments by Bearman and Zdravkovich [44], two values of e/D were examined; i.e., e/D = −0.4, −0.8; and the examined Reynolds number was ( ) f Re UD ν = = 1.5 × 10 4 . Note that the pipe was located at x/D = 0; the negative value of e/D denoted that the pipe was placed above the bed surface. Similarly, Figure 4b shows the numerical results of the distributions of Cp under the conditions of e/D = 0 (i.e., the pipe was just touching the bed surface) and Re = 4.5 × 10 4 . The experimental measurements by Tsiolakis [45] and the previous numerical results [46][47][48] are also provided in

Verification of Numerical Model
In this section, the proposed numerical model is validated by comparing with the existing experimental or numerical results. Both the pressure distribution at the mudline near the pipe and the critical velocity for the lateral instability of a pipe are examined, respectively.

Pressure Distribution at the Mudline near the Pipe
The proposed model was firstly validated against the existing experimental data regarding pressure distribution at the mudline near the pipe laid above or on the bed surface. The calculated distributions of pressure coefficient (C p ) at the mudline in the proximity of the pipe are plotted in Figure 4a, in which the experimental data of Bearman and Zdravkovich [44] are also provided. The pressure coefficient C p is defined as follows: where P s is the steady-flow-induced pressure at the mudline around the pipe, and P r is the reference pressure. In the experiments by Bearman and Zdravkovich [44], two values of e/D were examined; i.e., e/D = −0.4, −0.8; and the examined Reynolds number was Re(= UD/ν f ) = 1.5 × 10 4 . Note that the pipe was located at x/D = 0; the negative value of e/D denoted that the pipe was placed above the bed surface. Similarly, Figure 4b shows the numerical results of the distributions of C p under the conditions of e/D = 0 (i.e., the pipe was just touching the bed surface) and Re = 4.5 × 10 4 . The experimental measurements by Tsiolakis [45] and the previous numerical results [46][47][48] are also provided in this figure for comparison. As shown in Figure 4, the results of C p calculated with the present model overall agreed well with the existing experimental and numerical results, indicating that the present model was capable of predicting the flow-induced pressure on the seabed surface in the vicinity of the pipeline. Experimental data [45] Numerical results [46] Numerical results [47] Numerical results [48] Present numerical results

Critical Velocity for the Lateral Instability of a Pipe
To further verify the proposed numerical model, the steady-current-induced lateral instability of a freely laid pipe was investigated. Recently, a series of tests was conducted  Figure 5 shows the numerically simulated and the measured critical velocity for the lateral instability (Ucr) with various submerged weights of the pipe (Ws) on the silt bed for two values of e/D; i.e., e/D = 0.01, 0.05. As shown in Figure 5, the present numerical results were in reasonable agreement with the data for flume measurements [49]. The  Figure 5 shows the numerically simulated and the measured critical velocity for the lateral instability (U cr ) with various submerged weights of the pipe (W s ) on the silt bed for two values of e/D; i.e., e/D = 0.01, 0.05. As shown in Figure 5, the present numerical results were in reasonable agreement with the data for flume measurements [49]. The critical flow velocity for pipe instability (U cr ) increased approximately linearly with the submerged weight of the pipe (W s ) for a given value of e/D. For a fixed value of W s , the critical flow velocity U cr rose for a larger pipe embedment (e/D). It should be noted that, for the examined fine-grained silts with certain cohesion, tunnel erosion beneath the pipe was not observed in the flume experiments, even for the small embedment; e.g., e/D = 0.01. The proposed coupled FE model was capable of predicting the lateral-instability process by employing the aforementioned contact-pair algorithm for describing the pipe-soil interfacial behavior (see Section 2.

Flow over and Seepage Flow Beneath a Partially Embedded Pipe
The proposed coupled FE model was employed to explore the competition mechanism between the lateral instability and the tunnel erosion of a freely laid pipeline on a non-cohesive seabed. The input data of all parameters for the sandy seabed, the steady flow, and the pipe are listed in Table 1.
As shown in Figure 6a, for a certain inflow velocity (U = 1.0 m/s), the flow over and the seepage flow beneath the partially embedded pipe could be obtained synchronously, so that the pressure distribution at the water-soil interface was kept continuous. Similar to the observations by Mao [50], due to the presence of the pipe, a large lee-wake vortex was generated at the downstream of pipe, and two small vortices were formed in front and at the rear of the pipe, respectively (see Figure 6b). The flow pressure in front of the pipe was higher than that at the rear of it, which further induced the seepage flow within the soil beneath the pipe. Figure 6b,c give the seepage-flow direction and the contour of hydraulic gradients (i) around the pipe, respectively, indicating the peak values of the hydraulic gradient that emerged at the upstream and downstream corners of the pipe-soil interface (points A and B in Figure 6c). The seepage vectors were generally upward at the downstream zone (see Figure 6b), which could lead to an obliquely upward seepage flow at the seepage exit (i.e., point B in Figure 6c).

Flow over and Seepage Flow Beneath a Partially Embedded Pipe
The proposed coupled FE model was employed to explore the competition mechanism between the lateral instability and the tunnel erosion of a freely laid pipeline on a noncohesive seabed. The input data of all parameters for the sandy seabed, the steady flow, and the pipe are listed in Table 1.
As shown in Figure 6a, for a certain inflow velocity (U = 1.0 m/s), the flow over and the seepage flow beneath the partially embedded pipe could be obtained synchronously, so that the pressure distribution at the water-soil interface was kept continuous. Similar to the observations by Mao [50], due to the presence of the pipe, a large lee-wake vortex was generated at the downstream of pipe, and two small vortices were formed in front and at the rear of the pipe, respectively (see Figure 6b). The flow pressure in front of the pipe was higher than that at the rear of it, which further induced the seepage flow within the soil beneath the pipe. Figure 6b,c give the seepage-flow direction and the contour of hydraulic gradients (i) around the pipe, respectively, indicating the peak values of the hydraulic gradient that emerged at the upstream and downstream corners of the pipe-soil interface (points A and B in Figure 6c). The seepage vectors were generally upward at the downstream zone (see Figure 6b), which could lead to an obliquely upward seepage flow at the seepage exit (i.e., point B in Figure 6c).   Figure 7a,b give the pressure distributions along the mudline (P S ) and around the pipe surface (P P ) for various inflow velocities (U), respectively. It is indicated that the flow-induced pressure was generally positive at the upstream of the partially embedded pipe, while negative at the downstream. As the flow velocity increased, the pressure drop increased gradually. The hydrodynamic forces on the partially embedded smooth pipe could be calculated by integrating the pressure distribution around the pipe surface in the water domain (see Figure 2a), which could be decomposed into the in-line (drag force F D ) and cross-flow (lift force F L ) components. Figure 8 gives the variations of drag force on the pipe (F D ) and the hydraulic gradient at the seepage exit (i ex ) within the soil. As shown in Figure 8, the values of both F D and i ex increased nonlinearly with the increase of flow velocity (U), which played decisive roles for triggering either the lateral instability or the tunnel erosion of the pipe.  Figure 7a,b give the pressure distributions along the mudline (PS) and around the pipe surface (PP) for various inflow velocities (U), respectively. It is indicated that the flow-induced pressure was generally positive at the upstream of the partially embedded pipe, while negative at the downstream. As the flow velocity increased, the pressure drop increased gradually. The hydrodynamic forces on the partially embedded smooth pipe could be calculated by integrating the pressure distribution around the pipe surface in the water domain (see Figure 2a), which could be decomposed into the in-line (drag force FD) and cross-flow (lift force FL) components. Figure 8 gives the variations of drag force on the pipe (FD) and the hydraulic gradient at the seepage exit (iex) within the soil. As shown in Figure 8, the values of both FD and iex increased nonlinearly with the increase of flow velocity (U), which played decisive roles for triggering either the lateral instability or the tunnel erosion of the pipe.

Competition between Tunnel Erosion and Lateral Instability
As stated above, the seepage failure beneath the pipe resulting from a pressure drop was recognized as the dominant cause of the onset of tunnel erosion. To quantitatively predict the onset of tunnel erosion under the action of a steady current, the critical hydraulic gradient for the vertical seepage failure, − , was adopted in previous analyses (e.g., Sumer et al. [18]; Zang et al. [47]; Zhang et al. [21]). Note that Gs is the specific gravity of soil grains. As indicated in Figure 6b, the direction of seepage flow at the exit (point B) was tangent to the pipe surface instead of upward vertically. As such, in the present study, the criterion for the oblique seepage failure proposed by Gao and Luo [20] was employed; i.e., the hydraulic gradient at the seepage exit (iex) was beyond the critical value (icr): where θ is the embedment angle of the pipe (see Figure 8).
With increasing flow velocity of ocean currents, if the lateral soil resistance (FR) was insufficient to balance the drag force (FD), i.e., FD > FR, the pipe would break out from its original location. As observed in the flume experiments (also see Gao et al. [14]), the pipe breakout (or the lateral instability) generally occurred suddenly with a relatively large lateral displacement. In the present numerical simulations, the criterion for the lateral instability of the pipe was adopted as follows: where UcL is the critical flow velocity for the lateral instability of the pipe, s is the lateral displacement of the pipe at an increasing flow velocity (U) with a velocity increment ). Figure 9a shows the lateral displacement of the pipe (s/D) and the maximum upward hydraulic gradient (iex) with the increase of flow velocity (U), respectively. It is indicated that the values of s/D increased gradually with the increase of U until its critical value was reached (e.g., UcL = 1.20 m/s). The embedded pipe moved obliquely upward progressively (see Figure 9b). As shown in Figure 9b, at a relatively lower flow velocity, the plastic-strain zone was generated dominantly beneath the pipe, especially at the pipe-mudline corners; with further increasing U, the plastic strains extended successively and the soil resistance (FR) was increased correspondingly, resulting in a squeezed plastic zone at the rear of the pipe. According to the criterion (Equation (21)), the critical flow velocity (UcL) for lateral instability could be obtained. It should be noted that, in the

Competition between Tunnel Erosion and Lateral Instability
As stated above, the seepage failure beneath the pipe resulting from a pressure drop was recognized as the dominant cause of the onset of tunnel erosion. To quantitatively predict the onset of tunnel erosion under the action of a steady current, the critical hydraulic gradient for the vertical seepage failure, i cr0 = (1 − n)(G s − 1), was adopted in previous analyses (e.g., Sumer et al. [18]; Zang et al. [47]; Zhang et al. [21]). Note that G s is the specific gravity of soil grains. As indicated in Figure 6b, the direction of seepage flow at the exit (point B) was tangent to the pipe surface instead of upward vertically. As such, in the present study, the criterion for the oblique seepage failure proposed by Gao and Luo [20] was employed; i.e., the hydraulic gradient at the seepage exit (i ex ) was beyond the critical value (i cr ): i ex > i cr = (sin θ + cos θ tan ϕ)i cr0 (20) where θ is the embedment angle of the pipe (see Figure 8).
With increasing flow velocity of ocean currents, if the lateral soil resistance (F R ) was insufficient to balance the drag force (F D ), i.e., F D > F R , the pipe would break out from its original location. As observed in the flume experiments (also see Gao et al. [14]), the pipe breakout (or the lateral instability) generally occurred suddenly with a relatively large lateral displacement. In the present numerical simulations, the criterion for the lateral instability of the pipe was adopted as follows: where U cL is the critical flow velocity for the lateral instability of the pipe, s is the lateral displacement of the pipe at an increasing flow velocity (U) with a velocity increment ∆U (=0.05 m/s), and δ is a judging constant (choosing δ = 2.0). Figure 9a shows the lateral displacement of the pipe (s/D) and the maximum upward hydraulic gradient (i ex ) with the increase of flow velocity (U), respectively. It is indicated that the values of s/D increased gradually with the increase of U until its critical value was reached (e.g., U cL = 1.20 m/s). The embedded pipe moved obliquely upward progressively (see Figure 9b). As shown in Figure 9b, at a relatively lower flow velocity, the plastic-strain zone was generated dominantly beneath the pipe, especially at the pipe-mudline corners; with further increasing U, the plastic strains extended successively and the soil resistance (F R ) was increased correspondingly, resulting in a squeezed plastic zone at the rear of the pipe. According to the criterion (Equation (21)), the critical flow velocity (U cL ) for lateral instability could be obtained. It should be noted that, in the process of such lateral instability, the seepage flow was generated concurrently within the soil (see Figure 6), although i ex was lower than its critical value for tunnel erosion (i.e.,i ex < i cr (=1.10); see Figure 9a) in this case study. That is, under the action of ocean currents, the competition between the tunnel erosion and the lateral instability always existed for a pipeline partially embedded in a non-cohesive seabed. process of such lateral instability, the seepage flow was generated concurrently within the soil (see Figure 6), although iex was lower than its critical value for tunnel erosion (i.e., ex cr i i < (= 1.10); see Figure 9a) in this case study. That is, under the action of ocean currents, the competition between the tunnel erosion and the lateral instability always existed for a pipeline partially embedded in a non-cohesive seabed.  Figure 9. Competition between tunnel erosion and lateral instability: Developments of (a) the pipe lateral displacement (s/D) and the maximum upward hydraulic gradient (i ex ) and (b) the lateral soil resistance (F R ) and the corresponding pipe embedment (e y /D) with increasing flow velocity (U).

Instability Envelope: Effects of the Initial Embedment and the Submerged Weight of the Pipe
As aforementioned, if the basic properties of the seabed and the pipeline are given, the critical flow velocity for the pipeline instability (U cr ) will be mainly related to the two dimensionless parameters; i.e., the embedment-to-diameter ratio (e/D) and the nondimensional submerged weight of the pipeline (G = W s / γ D 2 (see Shi and Gao [23]). The proposed coupled FE model was employed to examine the effects of e/D and G on U cr . Figure 10a,b show the distributions of the pressure along the pipe surface (P P ) and along the mudline near the pipe (P s ), respectively, for various values of e/D and the flow velocity U. As shown in Figure 10, the flow-induced pressure distributed non-uniformly along the pipe surface and along the mudline. For a fixed value of U (=0.5, or 1.0 m/s), the values of P P increased gradually with increasing e/D from 0.1 to 0.3, especially at the top of the pipe (0 < β < 180 • ; see Figure 10a); accordingly, the values of P s increased slightly (see Figure 10b). As aforementioned, if the basic properties of the seabed and the pipeline are given, the critical flow velocity for the pipeline instability (Ucr) will be mainly related to the two dimensionless parameters; i.e., the embedment-to-diameter ratio (e/D) and the non-dimensional submerged weight of the pipeline (G = ( ) 2 s ' W D γ (see Shi and Gao [23]). The proposed coupled FE model was employed to examine the effects of e/D and G on Ucr. Figure 10a,b show the distributions of the pressure along the pipe surface (PP) and along the mudline near the pipe (Ps), respectively, for various values of e/D and the flow velocity U. As shown in Figure 10 Figure 10). Such a reduction in the drag force on the pipe due to increasing embedment was suggested in the existing Det Norske Veritas recommended practice [10]. It seems that the pipe embedment (e/D) had a slight effect on the pressure drop (P s ) along the mudline (see Figure 10b), but its effect on the hydraulic gradient at the exit (i ex ) became much more significant due to the enlargement of the seepage path with a slight increase of e/D (e.g., from 0.1 to 0.3; see Figure 11b). J. Mar. Sci. Eng. 2021, 9, x FOR PEER REVIEW 18 of 25 Figure 11a,b show the variations of the drag force on the pipe (FD) and those of the hydraulic gradient at the seepage exit (iex) with the increase of flow velocity (U), respectively, for various values of e/D. It is indicated that both FD and iex increased nonlinearly with increasing U. For a fixed value of U, the values of both FD and iex decreased with increasing e/D, which could be attributed to the variations of pressure distributions of PP and Ps with e/D (see Figure 10). Such a reduction in the drag force on the pipe due to increasing embedment was suggested in the existing Det Norske Veritas recommended practice [10]. It seems that the pipe embedment (e/D) had a slight effect on the pressure drop (Ps) along the mudline (see Figure 10b), but its effect on the hydraulic gradient at the exit (iex) became much more significant due to the enlargement of the seepage path with a slight increase of e/D (e.g., from 0.1 to 0.3; see Figure 11b).  The submerged weight of the pipe (W s ) may induce plastic strains within the soil beneath the pipe (see Figure 9), and further, has influence on the lateral resistance (F R ) to the partially embedded pipe under the action of ocean currents. Figure 12 shows the variations of F R with the flow velocity (U) for various values of the non-dimensional submerged weight of the pipe (G). As indicated in Figure 12, for a fixed value of G, the higher values of F R could be mobilized to balance the drag force (F D ) in a flow with increasing U. A heavier pipe with larger G would become more stable; i.e., a larger value of F R was mobilized to resist a higher inflow velocity U.
J. Mar. Sci. Eng. 2021, 9, x FOR PEER REVIEW The submerged weight of the pipe (Ws) may induce plastic strains within beneath the pipe (see Figure 9), and further, has influence on the lateral resistance the partially embedded pipe under the action of ocean currents. Figure 12 sho variations of FR with the flow velocity (U) for various values of the non-dime submerged weight of the pipe (G). As indicated in Figure 12, for a fixed value o higher values of FR could be mobilized to balance the drag force (FD) in a flow w creasing U. A heavier pipe with larger G would become more stable; i.e., a large of FR was mobilized to resist a higher inflow velocity U. To quantitatively characterize the competition mechanism between lateral i ity and tunnel erosion, a series of simulations was conducted with the proposed In the following parametric study, the values of e/D and G were varied, and the va other parameters are given in Table 1. As aforementioned, following the calc procedure for the coupled flow-seepage-elastoplastic model (see Figure 3), the flow velocity (Ucr) for the lateral instability or the tunnel erosion could be obtain effects of the embedment-to-diameter ratio (e/D) and non-dimensional subm weight of the pipe (G) on the critical flow velocity (Ucr) were further examined n cally. Figure 13a,b show the variations of Ucr with e/D for various values of G, a variations of Ucr with G for various values of e/D, respectively. It is indicated in 13a that, for a given value of G, there existed a transition point from the "tunnel-e (T-mode) to the "lateral-instability" (L-mode) on the critical instability line; the erosion was more prone to being triggered than the lateral instability for smaller of e/D. With the increase of G, the values of Ucr and the corresponding e/D for the tion points were both increased concurrently. For a heavy pipe; e.g., G = 1.29, o T-mode was triggered (that is, the L-mode was suppressed) for the examined r e/D (0.01 < e/D < 0.30; see the blue line in Figure 13a). Similarly, as shown in Figu for a pipe with a relatively larger embedment; e.g., e/D = 0.25, the L-mode wa easier to trigger in the examined range of G (0.22 < G < 1.29).
The instability envelope for the flow-pipe-soil interaction could be further To quantitatively characterize the competition mechanism between lateral instability and tunnel erosion, a series of simulations was conducted with the proposed model. In the following parametric study, the values of e/D and G were varied, and the values of other parameters are given in Table 1. As aforementioned, following the calculation procedure for the coupled flow-seepage-elastoplastic model (see Figure 3), the critical flow velocity (U cr ) for the lateral instability or the tunnel erosion could be obtained. The effects of the embedment-to-diameter ratio (e/D) and non-dimensional submerged weight of the pipe (G) on the critical flow velocity (U cr ) were further examined numerically. Figure 13a,b show the variations of U cr with e/D for various values of G, and the variations of U cr with G for various values of e/D, respectively. It is indicated in Figure 13a that, for a given value of G, there existed a transition point from the "tunnel-erosion" (Tmode) to the "lateral-instability" (L-mode) on the critical instability line; the tunnel erosion was more prone to being triggered than the lateral instability for smaller values of e/D. With the increase of G, the values of U cr and the corresponding e/D for the transition points were both increased concurrently. For a heavy pipe; e.g., G = 1.29, only the T-mode was triggered (that is, the L-mode was suppressed) for the examined range of e/D (0.01 < e/D < 0.30; see the blue line in Figure 13a). Similarly, as shown in Figure 13b, for a pipe with a relatively larger embedment; e.g., e/D = 0.25, the L-mode was much easier to trigger in the examined range of G (0.22 < G < 1.29). It was implied that, for a lighter (smaller G) and more deeply embedded (larger e/D) pipe, lateral instability (L-mode) would be more prone to occur; otherwise, tunnel erosion (T-mode) would be more prone to occur.  The instability envelope for the flow-pipe-soil interaction could be further established, as shown in Figure 14, and could be described by three controlling parameters; i.e., U cr , e/D, and G. As illustrated in this figure, the entire instability envelope was a smoothly curved surface and had two components; i.e., Model-I: Lateral-instability (L-mode), and Model-II: Tunnel-erosion (T-mode); while in between T-mode and L-mode was a transition line. If the flow velocity of ocean currents went beyond the instability envelope, either the T-mode or L-mode of the pipe instability could be induced. It was implied that, for a lighter (smaller G) and more deeply embedded (larger e/D) pipe, lateral instability (L-mode) would be more prone to occur; otherwise, tunnel erosion (T-mode) would be more prone to occur.

Conclusions
Under the action of ocean currents, two typical instability modes could be triggered; i.e., the lateral instability of the pipe and the tunnel erosion of the underlying soil, which involve complex fluid-pipe-soil interactions. Nevertheless, these two instability modes were only investigated separately in the past. In this paper, the competition mechanism between lateral instability and tunnel erosion of a freely laid pipeline on a non-cohesive seabed was investigated numerically. The following conclusions could be drawn: 1. A coupled flow-seepage-elastoplastic modeling approach was proposed that could realize the synchronous simulation of the pipe hydrodynamics, the seepage flow, and the elastoplastic behavior of the soil beneath the pipe. The proposed model was verified through experimental and numerical results. 2. Both the drag force on the pipe and the hydraulic gradient at the seepage exit beneath the pipe increased concurrently with increasing flow velocity, which could induce either the lateral instability of pipe or the tunnel erosion of soil. The competition between such two instability modes always existed for a partially embedded pipe. 3. A parametric study indicated that there generally existed a transition point from tunnel erosion to lateral instability on the critical instability line representing the variation of critical flow velocity (Ucr) with embedment-to-diameter ratio (e/D) for various values of a non-dimensional submerged weight of the pipe (G), or the variation of Ucr with G for various values of e/D. Tunnel erosion was more prone to be-

Conclusions
Under the action of ocean currents, two typical instability modes could be triggered; i.e., the lateral instability of the pipe and the tunnel erosion of the underlying soil, which involve complex fluid-pipe-soil interactions. Nevertheless, these two instability modes were only investigated separately in the past. In this paper, the competition mechanism between lateral instability and tunnel erosion of a freely laid pipeline on a non-cohesive seabed was investigated numerically. The following conclusions could be drawn:

1.
A coupled flow-seepage-elastoplastic modeling approach was proposed that could realize the synchronous simulation of the pipe hydrodynamics, the seepage flow, and the elastoplastic behavior of the soil beneath the pipe. The proposed model was verified through experimental and numerical results.

2.
Both the drag force on the pipe and the hydraulic gradient at the seepage exit beneath the pipe increased concurrently with increasing flow velocity, which could induce either the lateral instability of pipe or the tunnel erosion of soil. The competition between such two instability modes always existed for a partially embedded pipe.

3.
A parametric study indicated that there generally existed a transition point from tunnel erosion to lateral instability on the critical instability line representing the variation of critical flow velocity (U cr ) with embedment-to-diameter ratio (e/D) for various values of a non-dimensional submerged weight of the pipe (G), or the variation of U cr with G for various values of e/D. Tunnel erosion was more prone to being triggered than lateral instability for smaller values of the embedment-to-diameter ratio. 4.
The instability envelope for the flow-pipe-soil interaction system was eventually established, and could be described by three parameters; i.e., U cr , e/D, and G. It was implied that, for a lighter (smaller G) and more deeply embedded (larger e/D) pipe, lateral instability would be more prone to occur; otherwise, tunnel erosion would be more prone to take place.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available upon reasonable request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.

Notations b i Body acceleration per unit mass c
Cohesion of soil C 1ε Constant in Equation (5) C 2ε Constant in Equation (5) C µ Constant in Equations (4) and (5) C p Pressure coefficient at the mudline in the proximity of the pipe d 50 Mean size of soil grains dλ Plastic multiplier in Equations (13)