E ﬀ ect of Gain in Soil Friction on the Walking Rate of Subsea Pipelines

: Subsea pipelines are commonly employed in the o ﬀ shore oil and gas industry to transport high-pressure and high-temperature (HPHT) hydrocarbons. The phenomenon of pipeline walking is a topic that has drawn a great deal of attention, and is related to the on-bottom stability of the pipeline, such as directional accumulation with respect to axial movement, which can threaten the security of the entire pipeline system. An accurate assessment of pipeline walking is therefore necessary for o ﬀ shore pipeline design. This paper reports a comprehensive suite of numerical analyses investigating the performance of pipeline walking, with a focus on the e ﬀ ect of increasing axial soil resistance on walking rates. Three walking-driven modes (steel catenary riser (SCR) tension, downslope, and thermal transient) are considered, covering a wide range of inﬂuential parameters. The variation in walking rate with respect to the e ﬀ ect of increased soil friction is well reﬂected in the development of the e ﬀ ective axial force (EAF) proﬁle. A method based on the previous analytical solution is proposed for predicting the accumulated walking rates throughout the entire service life, where the concept of equivalent soil friction is adopted. tension scenario; (2) downslope condition—vertical component G 1 of pipeline gravity G ( G 1 = G × cos φ , where φ is the angle of seabed slope) and tangential component G 2 of gravity G ( G 2 = G × sin φ ) are applied on the pipeline elements to simulate the inﬂuence of seabed slope; and (3) thermal transient condition—the user-deﬁned subroutine UTEMP is used to simulate the thermal transient process with temperature as a function of time and space coordinates. In this study, linear thermal transient proﬁles are adopted, considering the constant thermal gradient. For the continuous thermal transient, a su ﬃ ciently small temperature increment is applied in the static analysis, with total increments of 100 adopted in this study. Figure 3 shows an example of a thermal transient proﬁle for a 1000 m length of pipeline during the heating-up process. After the early 100 temperature increments in the heating-up analysis step, the point at x / L = 0 (where x is the position along the pipeline and L is the pipeline length) reaches the target temperature of T = 100 ◦ C with the point at x / L = 1 starting to be heated. The corresponding four separate lines extrapolated by every other 25 computational increments are plotted in Figure 3a for the proﬁles of thermal transient along the pipeline. After that, the other 100 temperature increments follow until the temperature of the whole pipeline arrives at the target. Figure 3b demonstrates the changes in temperature at positions x / L = 0, 0.5, and 1 with the interval of 25 computational increments throughout the entire heating-up analysis. SCR tension scenario; (2) downslope condition — vertical component G 1 of pipeline gravity G ( G 1 = G × cos ϕ , where ϕ is the angle of seabed slope) and tangential component G 2 of gravity G ( G 2 = G × sin ϕ ) are applied on the pipeline elements to simulate the influence of seabed slope; and (3) thermal transient condition — the user-defined subroutine UTEMP is used to simulate the thermal transient


Introduction
Submarine pipelines are important facilities for subsea production systems to transport the products of high-pressure and high-temperature (HPHT) hydrocarbon sources in the offshore oil and gas industry. On-bottom stability analysis is required in subsea pipeline design by a number of industry design codes (e.g., DNV GL-OS-F101, DNV GL-RP-F114), and this has been an attractive topic among the research community. The evaluation of global axial displacement of a pipeline, i.e., "pipeline walking" [1], throughout its entire service life has attracted increasing attention, since significant accumulation of axial displacement could threaten the safety of not only the pipeline itself, but also the end connections and other auxiliary flowline devices. For instance, an accumulated axial displacement of over 7 m was unexpectedly observed for a 2000 m long pipeline in the North Sea [2], leading to the collapse of the end connections.
The early development of hydrocarbon sources was mainly localized nearshore, with pipelines fully buried in the seabed, and sometimes paved with stones, so that the pipeline could be restricted, resulting in negligible movement. This may not, however, be the cause of the phenomenon of pipeline walking. There have only been a few investigations related to potential axial movements in pipelines, focusing on aspects such as the relative movement between the pipeline and the concrete coating [3], and the soil friction along the pipeline [4,5].

Scope of This Study
Subsea pipelines require inspection and repair during their service life, therefore intermittent start-up and shut-down are generally applied. When the pipeline is started up, HPHT contents flow through the whole pipeline, which leads to obvious increases in temperature and pressure of the flowline, thus axial expansion in a short pipeline due to insufficient constraint force from surrounding soil friction. In contrast, the temperature and pressure of the flowline decrease during the shut-down operation, which can cause axial shrinkage. If the pipeline expands and shrinks symmetrically about the midpoint throughout an operational procedure of start-up following shut-down, the pipeline exhibits no accumulation of axial movement. However, complex loading conditions of pipelines, such as SCR tension force loaded at one endpoint, seabed slope along the flowline, thermal transient along the flowline route, etc. [1], could lead to asymmetry in the EAF profile. As a result, deformations of expansion and shrinkage are symmetric about different points. In addition, the direction of expansion in some sections of the pipeline during the start-up stage could be the same as that during the shut-down stage, which can lead to globally accumulated axial movement after the operational loading cycle.
Subsea pipelines are subjected to a significant amount of start-up and shut-down cycles [16][17][18] and may suffer large walking rates during their whole service life. An operational loading cycle is generally characterized by four typical phases: Heating up, operation, cooling down, and maintenance. Figure 1a demonstrates possible behaviors for axial pipe-soil interactions during an operational loading cycle, where a segment of pipeline is considered to be already laid on the clayey seabed (common in offshore areas). More details are provided below. flowline, thermal transient along the flowline route, etc. [1], could lead to asymmetry in the EAF profile. As a result, deformations of expansion and shrinkage are symmetric about different points. In addition, the direction of expansion in some sections of the pipeline during the start-up stage could be the same as that during the shut-down stage, which can lead to globally accumulated axial movement after the operational loading cycle. Subsea pipelines are subjected to a significant amount of start-up and shut-down cycles [16][17][18] and may suffer large walking rates during their whole service life. An operational loading cycle is generally characterized by four typical phases: Heating up, operation, cooling down, and maintenance. Figure 1a demonstrates possible behaviors for axial pipe-soil interactions during an operational loading cycle, where a segment of pipeline is considered to be already laid on the clayey seabed (common in offshore areas). More details are provided below.   Heating-up phase: The axial soil resistance is first mobilized, and excess pore pressure (i.e., additional generation of pore water pressure due to undrained shearing that is defined relative to the in situ geostatic water pressure [14,19,20]) around the pipe accumulates synchronously.
Operation phase: The excess pore pressure is dissipated, with the temperature and pressure of the whole pipe being maintained without any axial movement or axial shearing. This causes a gain in axial soil resistance due to consolidation of soil around the pipe.
Cooling-down phase: After operation, further axial movement of the pipe occurs due to shrinkage during the cooling-down phase, leading to the generation of additional excess pore pressure.
Maintenance phase: In this phase, the dissipation of excess pore pressure leads to a regain in axial soil resistance, which will alter the walking rate during the subsequent heating-up phase in the next cycle.
As is evident, throughout a pipeline's operational loading cycle, the surrounding soil could be enhanced during the intermittent periods, such as in the operation and maintenance phases, which is mainly attributed to the occurrence of consolidation in these stages [15,[21][22][23][24]. Hence, reasonable adoption of soil friction in the design of subsea pipelines is necessary in order to capture a reliable estimation of walking rate. The concept of equivalent friction coefficient is employed, which is commonly used in the design code of subsea pipelines, such as DNV GL-RP-F114, and other research related to pipeline stability, such as [1,10,[25][26][27][28][29]. Following their descriptions, the equivalent friction coefficient is the production of axial soil resistance F divided by pipeline submerged weight W per unit length of pipeline (µ = F/W ). However, µ is generally taken as constant through an operational loading cycle, neglecting the reconsolidation effect during the intermittent operation and maintenance phases. An updated concept of equivalent friction coefficient is shown in Figure 1b, which is correlated with the episodic nature of realistic pipeline scenarios. Figure 1b shows a concept of the equivalent soil friction correlated with the episodic nature of realistic pipeline operation. For instance, if soil friction coefficient µ 1 is adopted in pipeline design in the heating-up phase, a greater soil friction coefficient µ 2 can be used on the pipeline in the following cooling-down phase due to the enhanced soil resistance.
Considering such gains in soil friction, the pipeline's walking behavior can be changed. In the following sections, the development of walking rate during a typical operational loading cycle is investigated and interpreted. A wide range of gains in axial friction are considered, with other key parameters that affect pipeline walking behavior also included. Then, the corresponding evaluation of walking rate is implemented and extended to the scenario of several operational cycles.

Finite Element Modelling
The entire analysis was undertaken using commercially available Abaqus software with standard analyses adopted [30].

Model Description
The pipeline is prescribed with the element type PIPE31. It is only allowed to move axially, with the lateral movement fixed. The seabed is simulated with the element type C3D8R, with all sides and bottom boundaries restricted. The mesh densities of both pipeline and seabed soil are examined by mesh sensitivity analysis, with the results demonstrated in Figure 2b. As is evident, the examined walking rate is independent of the mesh length of the seabed soil, and shows convergence until the mesh length of the pipeline is less than 0.5 m. It is reasonable, since the seabed soil only provides a rigid boundary. The axial interaction force is fully transferred from seabed to pipeline by a setup of pipe-soil interface. More detailed information will be described later. The loading conditions for three driven modes of pipeline walking are specified as follows: (1) SCR tension condition-a concentrated force is applied on one end of the pipeline to simulate the The loading conditions for three driven modes of pipeline walking are specified as follows: (1) SCR tension condition-a concentrated force is applied on one end of the pipeline to simulate the SCR tension scenario; (2) downslope condition-vertical component G 1 of pipeline gravity G (G 1 = G × cosφ, where φ is the angle of seabed slope) and tangential component G 2 of gravity G (G 2 = G × sinφ) are applied on the pipeline elements to simulate the influence of seabed slope; and (3) thermal transient condition-the user-defined subroutine UTEMP is used to simulate the thermal transient process with temperature as a function of time and space coordinates. In this study, linear thermal transient profiles are adopted, considering the constant thermal gradient. For the continuous thermal transient, a sufficiently small temperature increment is applied in the static analysis, with total increments of 100 adopted in this study. Figure 3 shows an example of a thermal transient profile for a 1000 m length of pipeline during the heating-up process. After the early 100 temperature increments in the heating-up analysis step, the point at x/L = 0 (where x is the position along the pipeline and L is the pipeline length) reaches the target temperature of T = 100 • C with the point at x/L = 1 starting to be heated. The corresponding four separate lines extrapolated by every other 25 computational increments are plotted in Figure 3a for the profiles of thermal transient along the pipeline. After that, the other 100 temperature increments follow until the temperature of the whole pipeline arrives at the target. Figure  SCR tension scenario; (2) downslope condition-vertical component G1 of pipeline gravity G (G1 = G × cosϕ, where ϕ is the angle of seabed slope) and tangential component G2 of gravity G (G2 = G × sinϕ) are applied on the pipeline elements to simulate the influence of seabed slope; and (3) thermal transient condition-the user-defined subroutine UTEMP is used to simulate the thermal transient process with temperature as a function of time and space coordinates. In this study, linear thermal transient profiles are adopted, considering the constant thermal gradient. For the continuous thermal transient, a sufficiently small temperature increment is applied in the static analysis, with total increments of 100 adopted in this study. Figure 3 shows an example of a thermal transient profile for a 1000 m length of pipeline during the heating-up process. After the early 100 temperature increments in the heating-up analysis step, the point at x/L = 0 (where x is the position along the pipeline and L is the pipeline length) reaches the target temperature of T = 100 °C with the point at x/L = 1 starting to be heated. The corresponding four separate lines extrapolated by every other 25 computational increments are plotted in Figure 3a for the profiles of thermal transient along the pipeline. After that, the other 100 temperature increments follow until the temperature of the whole pipeline arrives at the target. Figure 3b demonstrates the changes in temperature at positions x/L = 0, 0.5, and 1 with the interval of 25 computational increments throughout the entire heating-up analysis.  For modelling of the pipe-soil interaction, the general contact algorithm obeying the Coulomb friction law is used to specify the pipe-soil interface [26][27][28][31][32][33][34][35][36]. To represent a gain in axial soil friction due to the consolidation hardening effect during the intermittent operation and maintenance phases, the equivalent frictional coefficient µ described in Figure 1b is used.
Regarding a typical operational loading cycle, an initial frictional coefficient µ 1 is specified in the heating-up phase, then a larger µ 2 is applied in the following cooling-down phase. These yield a gain in axial friction of η = µ 2 /µ 1 . A wide range of η is considered in the modelling to investigate the effect of η on walking rate during a typical loading cycle. Also, the conventional configuration of µ is presented for comparison, where no change in µ occurs through the typical loading cycle.
Regarding the simulation of accumulated walking rates through the whole service life of the pipeline, a hardening rule representing the gain in axial friction with the pipeline operational cycles is necessary. Yan et al. [15] proposed a relationship between the gain in undrained shear strength and axial shearing cycles (note that a pipeline operational cycle comprises two axial shearing cycles, associated with a heating-up followed by a cooling-down phase), which can be used to represent the equivalent gain in axial friction throughout the shearing cycles, given by where µN is the equivalent axial frictional coefficient, indicating the values of frictional coefficient in the heating-up and cooling-down phases in a specified operational cycle; µ 1 is the initial equivalent axial frictional coefficient, associated with the initial intact undrained shear strength of the seabed; N is the shearing cycles, with one operational cycle comprising two shearing cycles, corresponding to the heating-up and cooling-down phases; R is a state parameter defined in the critical state soil model; and κ/λ is the elasto-plastic volumetric stiffness ratio. Here, R and κ/λ can be taken as 3 and 0.2, respectively, for typical marine clay. Although their variation may influence the response of the equivalent friction coefficient, that is beyond the scope of this study. Following Equation (1), a user FRIC subroutine is adopted in the simulation of accumulated walking rates to replicate the nonlinear mobilization of equivalent soil friction throughout the whole life of the pipeline. The key process of setting up the FRIC subroutine is to capture a reasonable equivalent friction coefficient against operational loading cycles. Figure 4 shows the transformation process from gain in soil friction against shearing cycles to that against operational loading cycles. Following the description of Yan et al. [15], each operational loading cycle comprises two shearing cycles, associated with a heating-up and a cooling-down phase. The relationship between them is indicated in Figure 4, where the dashed black line represents the response of gain in friction against shearing cycles and the dashed blue line its corresponding gain in friction against operational loading cycles. The traditional adoption of equivalent friction coefficient is also included, which is a constant response through cycles.

Key Parameters Considered in Modelling
Throughout the modelling of pipeline walking, a relative pipeline walking rate β is introduced, given by where, SR and ST are pipeline walking rates per operational loading cycle with and without considering the gain in axial soil friction, respectively. They are directly captured from FEM (finite element modelling) analysis in this study. Following the description by Carr et al. [1] of pipeline walking, the axial displacements at the middle point of the pipeline in the analysis are used to indicate the walking rates. If the value of β is negative, it means the gain in axial soil resistance during pipeline operation can reduce the walking rate, and vice versa.
The ranges of 11 parameters adopted throughout the entire finite element analysis are listed in Table 1. The parameters involved in the examination of walking rate include six related to pipeline properties, three relevant to walking-driven mode, and two soil property parameters. The six pipeline parameters describe the basic information covering cross-section area A, Young's modulus E, thermal expansion coefficient α, submerged weight of pipeline per unit length W′, temperature difference ΔT, and total length of the whole pipeline L. The three key factors triggering walking are SCR tension force FSCR, seabed slope ϕ, and temperature gradient q. The two soil property parameters are initial soil friction coefficient μ1 and its gain η. Note that the ranges for each parameter tabulated in Table 1 represent typical conditions in offshore pipeline design, and refer to [1,2,9,10,18,25,26]. For instance, the length of the pipeline is taken as 1 to 4 km, since the pipeline cannot be too long to reach a condition of full constraint, which is traditionally less than 5 km [25].

Key Parameters Considered in Modelling
Throughout the modelling of pipeline walking, a relative pipeline walking rate β is introduced, given by where, S R and S T are pipeline walking rates per operational loading cycle with and without considering the gain in axial soil friction, respectively. They are directly captured from FEM (finite element modelling) analysis in this study. Following the description by Carr et al. [1] of pipeline walking, the axial displacements at the middle point of the pipeline in the analysis are used to indicate the walking rates. If the value of β is negative, it means the gain in axial soil resistance during pipeline operation can reduce the walking rate, and vice versa. The ranges of 11 parameters adopted throughout the entire finite element analysis are listed in Table 1. The parameters involved in the examination of walking rate include six related to pipeline properties, three relevant to walking-driven mode, and two soil property parameters. The six pipeline parameters describe the basic information covering cross-section area A, Young's modulus E, thermal expansion coefficient α, submerged weight of pipeline per unit length W , temperature difference ∆T, and total length of the whole pipeline L. The three key factors triggering walking are SCR tension force F SCR , seabed slope φ, and temperature gradient q. The two soil property parameters are initial soil friction coefficient µ 1 and its gain η. Note that the ranges for each parameter tabulated in Table 1 represent typical conditions in offshore pipeline design, and refer to [1,2,9,10,18,25,26]. For instance, the length of the pipeline is taken as 1 to 4 km, since the pipeline cannot be too long to reach a condition of full constraint, which is traditionally less than 5 km [25]. Table 1. Ranges of adopted parameters for parametric analysis.

Procedure of Numerical Analysis
To capture the relative walking rate β under conditions with combinations of the above described parameters, the simulation procedure is as follows: (1) First, a reference case is implemented with all parameters choosing the middle values in their corresponding ranges tabulated in Table 1; (2) the value of η is maintained at its minimum value, η = 1, and other parameters start changing gradually within their ranges, to investigate their individual influence and to capture the values of S T ; and (3) the value of η starts to change in its given range, and different values of other parameters are adopted for each change in η to examine the responses of S R . There are 700 simulation cases in total in this study.

Model Validation
To validate the proposed finite element model, the walking rates caused by SCR tension, downslope, and thermal transient captured from numerical modelling are compared with those calculated using an analytical method [1]. The pipeline walking rates driven by SCR tension, seabed slope, and thermal transient in the analytical solutions are presented in Equations (3)-(5), respectively, expressed as where ∆S k is axial displacement at the center of the pipeline, given by with where k is the number of heating increments and x A is the location of the anchor point indicating the position of the pipeline with zero movement through the heating-up stage, and the initial value of x θk is taken with x θ0 = x A (more details presented in Carr et al. [1]). Note that the gain in axial soil friction cannot be accounted for in the analytical solutions [1], therefore a constant soil friction coefficient is adopted in the validation analysis (such that the initial frictional coefficient µ 1 is taken through the analysis). The relevant input parameters of numerical modelling are listed in Table 2.

Variation of β Considering Gain in Axial Friction of Pipeline
The following three subsections mainly focus on the performance of relative pipeline walking rate β during one typical operational loading cycle, examining the individual effects of multiple influential factors on β. For each case demonstrated in these subsections, only one corresponding parameter changes in the analysis, with other parameters taken as the middle value given in Table 1, except for a value of η = 1.1 adopted throughout all simulation cases, denoted as "Ref. case" in this study. Figure 6 shows the variation of β with different influential factors under the SCR tension condition. As is evident, β is negative in each case, implying that the gain in axial soil friction can lead to a reduction in pipeline walking rate. Compared to the cases with different E, the variations in parameters of α, μ1, ΔT, W′, L, and A cause more significant changes in β. As α, ΔT, and A increase, the reductions in pipeline walking rates decrease (see the gradual increase of negative β in Figure 6), while with increasing μ1, W′, and L, the reductions in walking rates increase. The changing SCR tension (see the curve labeled "FSCR") has a negligible influence on β. In addition, Figure 6b shows two curves of β with μ1 (one is for the typical case with A taken as the middle value of 0.02 m 2 , the other for an additional case with A = 0.03 m 2 ). They show different responses of β, and this phenomenon can be also found in other examined groups. It means that the influences on β of all parameters adopted here are not independent. The coupling effect of the parameters evidently exists. The results of pipeline walking rate captured from finite element analysis reported in [10,28,37] are also demonstrated in Figure 5. The parameters adopted in each case are listed in Table 2. Compared to the corresponding results obtained by the proposed numerical model in this paper, it is clear that very close responses can be found in general.

Under the SCR Tension Condition
These provide sufficient evidence to prove the robustness of the numerical model used in this study.

Variation of β Considering Gain in Axial Friction of Pipeline
The following three subsections mainly focus on the performance of relative pipeline walking rate β during one typical operational loading cycle, examining the individual effects of multiple influential factors on β. For each case demonstrated in these subsections, only one corresponding parameter changes in the analysis, with other parameters taken as the middle value given in Table 1, except for a value of η = 1.1 adopted throughout all simulation cases, denoted as "Ref. case" in this study. Figure 6 shows the variation of β with different influential factors under the SCR tension condition. As is evident, β is negative in each case, implying that the gain in axial soil friction can lead to a reduction in pipeline walking rate. Compared to the cases with different E, the variations in parameters of α, µ 1 , ∆T, W , L, and A cause more significant changes in β. As α, ∆T, and A increase, the reductions in pipeline walking rates decrease (see the gradual increase of negative β in Figure 6), while with increasing µ 1 , W , and L, the reductions in walking rates increase. The changing SCR tension (see the curve labeled "F SCR ") has a negligible influence on β. In addition, Figure 6b shows two curves of β with µ 1 (one is for the typical case with A taken as the middle value of 0.02 m 2 , the other for an additional case with A = 0.03 m 2 ). They show different responses of β, and this phenomenon can be also found in other examined groups. It means that the influences on β of all parameters adopted here are not independent. The coupling effect of the parameters evidently exists.   Figure 7 shows the variation of β against different influential factors under the downslope condition. Comparing the SCR tension condition ( Figure 5) and the downslope condition (Figure 6), the main findings can be summarized as follows: (1) The gain in axial soil friction causes a reduction in pipeline walking rate for both cases; and (2) the increasing α, ∆T, and A reduce the decrement in the pipeline walking rate, and the increasing µ 1 , W , and L give negative increases in walking rates, but the magnitude of change in β for the downslope condition differs from that in the SCR tension condition; and (3) the changes in seabed slope φ have a minimum influence on β.  Figure 7 shows the variation of β against different influential factors under the downslope condition. Comparing the SCR tension condition ( Figure 5) and the downslope condition (Figure 6), the main findings can be summarized as follows: (1) The gain in axial soil friction causes a reduction in pipeline walking rate for both cases; and (2) the increasing α, ΔT, and A reduce the decrement in the pipeline walking rate, and the increasing μ1, W′, and L give negative increases in walking rates, but the magnitude of change in β for the downslope condition differs from that in the SCR tension condition; and (3) the changes in seabed slope ϕ have a minimum influence on β.

Under the Thermal Transient Condition
For the thermal transient condition, the responses of β to different influential factors are shown in Figure 8. The effect of changes in temperature ∆T is not considered here, because thermal transient induced walking is independent of ∆T. Compared with the response of β under SCR tension and downslope conditions, the pipeline presents quite different performance under the thermal transient condition. Both positive and negative β can be observed with the changing influential factors. In particular, for increasing q, a nonmonotonic response of β can be found, with a sharp drop from positive to negative following a relatively moderate decrease in the negative range. A monotonic response of β to the changing influential factors is also obtained. A positive decrease in β is found in the case of increasing µ 1 and E, while a positive increase in β is observed with increasing α. A negative decrease in β is found with increasing W . In addition, the increase of L has a negligible influence on the change in β.

Interpretation of Development of Relative Walking Rate β
As described above, relative pipeline walking rate β shows significantly complex responses when multiple influential factors are involved. Yet, the development of walking rate is essentially related to the build-up of the pipeline EAF. The EAF profile is therefore employed here, with its reference profile taken from the Ref. case. The examination of differences in EAF profile due to changes of influential factors is implemented in this section, in an attempt to interpret the response of β.
For the SCR tension condition, Figure 9a shows a group of EAF profiles, with dark lines for the Ref. case and others for individual cases where only one parameter changes compared to the Ref. case (where η = 1.1). The equivalent thermal loads per unit length of pipeline, which are fully constrained forces as well (see straight lines in Figure 9a), are also included. Considering a typical operational loading cycle, the development of EAF through a cycle can be expressed as a nondimensional quantity, given by ∆s 1 + ∆s 1 where s 1 + s 1 is the sum of maximum EAF under both heating-up and cooling-down phases, and ∆s 2 + ∆s 2 is the sum of thermal loads in these two phases, each of which is denoted as the fully constrained force [1]. It is in fact an axial force applied on the pipeline, ensuring that the axial deformation along the cross-sectional pipeline is fully restricted. Its value is equal to EAα∆T, representing the thermal load applied to the pipeline. Under the SCR tension condition, both ∆s 2 and ∆s 2 are constant throughout the heating-up and cooling-down phases. They are indicated in Figure 9a.

Interpretation of Development of Relative Walking Rate β
As described above, relative pipeline walking rate β shows significantly complex responses when multiple influential factors are involved. Yet, the development of walking rate is essentially related to the build-up of the pipeline EAF. The EAF profile is therefore employed here, with its reference profile taken from the Ref. case. The examination of differences in EAF profile due to changes of influential factors is implemented in this section, in an attempt to interpret the response of β.
For the SCR tension condition, Figure 9a shows a group of EAF profiles, with dark lines for the Ref. case and others for individual cases where only one parameter changes compared to the Ref. case (where η = 1.1). The equivalent thermal loads per unit length of pipeline, which are fully constrained forces as well (see straight lines in Figure 9a), are also included. Considering a typical operational loading cycle, the development of EAF through a cycle can be expressed as a nondimensional quantity, given by where s1 + s′1 is the sum of maximum EAF under both heating-up and cooling-down phases, and Δs2 + Δs′2 is the sum of thermal loads in these two phases, each of which is denoted as the fully constrained force [1]. It is in fact an axial force applied on the pipeline, ensuring that the axial deformation along the cross-sectional pipeline is fully restricted. Its value is equal to EAαΔT, representing the thermal load applied to the pipeline. Under the SCR tension condition, both Δs2 and Δs′2 are constant throughout the heating-up and cooling-down phases. They are indicated in Figure  9a.   As evident from both Figure 9a and Equation (8), compared to the Ref. case, increased E, A, α, and ∆T result in a reduction of (s 1 + s 1 )/(s 2 + s 2 ), which implies that the applied thermal load increases relative to soil friction (here it is equal to EAF). This results in increased deformation induced by thermal load. Such an increase in thermal deformation could contribute to a rising walking rate. Considering an identical EAF profile (as gain in axial friction η = 1.1), the relative walking rate β is a negative decrease. For cases of increased L, µ 1 , and W', the increase of (s 1 + s 1 )/(s 2 + s 2 ) reflects a relative rise in soil friction compared to the applied thermal load, generating lower walking rates and a negative increase of β.
Similarly, for the downslope condition, the nondimensional quantity for the EAF development can be expressed as Using Equation (9), the responses of walking rates under the downslope condition can be well explained similar to those under the SCR tension condition.
For the thermal transient condition, the EAF profiles after 20 temperature increments in the heating-up phase and their corresponding thermal loads are plotted in Figure 9c. The EAF profile at the completion of the last cooling-down phase (see dashed dark line in Figure 9c) is also included as a backbone for analysis. It is considered that no walking rate is generated in the first operational loading cycle under the thermal transient condition. The nondimensional quantity of EAF build-up can be given by Here, s 1 + s 1 and s 2 + s 2 are selected from the sections of activated pipeline due to thermal transient, which are both indicated in Figure 9c. As evident from Equation (10), L and ∆T have no influence on walking rate, and the development of walking rate is also well interpreted similar to the above described conditions.

Evaluation of Relative Walking Rate β
With respect to Equations (8)- (10), relative loading ratio f r with a strong linkage to β can be found, and can be expressed as where f 1 is the initial axial soil resistance acting on the unit length of the pipeline, given by In addition, f * is the thermal load acting on the unit length of the pipeline, given by for SCR tension or downslope condition Based on the described relative loading ratio f r , all results obtained from the parametric analysis are reprocessed, and the responses of β are replotted against f r . Figure 10 reports the variation of β in response to f r under different gains in axial soil friction η for the SCR tension condition. A series of clear linearly decreasing curves can be found, with the maximum β of~−5% at η = 1.1 and f r = 0.06, and the minimum β of~−26% at η = 1.5 and f r = 0.62. This implies that the higher η (due to intermittent consolidation) and the relative loading ratio are, the lower the pipeline walking rate will be. These tendencies are expressed by the curve-fitting process: where R 1 and R 2 are the fitting factors, given by Good agreement between the results from numerical analysis and calculated by Equation (14) is obtained, with the maximum discrepancy being less than 3%.
Based on the described relative loading ratio fr, all results obtained from the parametric analysis are reprocessed, and the responses of β are replotted against fr. Figure 10 reports the variation of β in response to fr under different gains in axial soil friction η for the SCR tension condition. A series of clear linearly decreasing curves can be found, with the maximum β of ∼−5% at η = 1.1 and fr = 0.06, and the minimum β of ∼−26% at η = 1.5 and fr = 0.62. This implies that the higher η (due to intermittent consolidation) and the relative loading ratio are, the lower the pipeline walking rate will be. These tendencies are expressed by the curve-fitting process: 1 r 2 R f R Δ = + (14) where R1 and R2 are the fitting factors, given by Good agreement between the results from numerical analysis and calculated by Equation (14) is obtained, with the maximum discrepancy being less than 3%.  Figure 11 demonstrates the responses of β against fr at different η for the downslope condition. A similar tendency of β can be observed compared to that for the SCR tension condition. The maximum β for the downslope condition achieved at η = 1.5 and fr = 0.62 is slightly higher than that Eq.14 Figure 10. Relative walking rate β against relative loading ratio f r under SCR tension condition. Figure 11 demonstrates the responses of β against f r at different η for the downslope condition. A similar tendency of β can be observed compared to that for the SCR tension condition. The maximum β for the downslope condition achieved at η = 1.5 and f r = 0.62 is slightly higher than that for the SCR tension condition, but the maximum β is almost identical. Considering a better estimation of β, these responses can be expressed as where R 3 and R 4 are the fitting factors, given by: 3 r 4 β R f R = + (17) where R3 and R4 are the fitting factors, given by:  Figure 11. Relative walking rate β against relative loading ratio fr under downslope condition. Figure 12 illustrates the responses of β in correlation with fr at different η for the thermal transient condition. It is clear that all responses of β at different η intersect at the axis of β = 0 when fr is equal to around 0.15. With increasing η, positive Δ shows an increasing tendency when fr is lower than 0.15. The maximum β can be found at η = 1.5 and fr = 0.01. With increasing η, negative β shows a decreasing tendency when fr is lower than 0.15. The minimum β can be found at η = 1.5 and fr = 0.44. This implies that pipelines with larger gain in axial friction have a larger walking rate when fr < 0.15, and a lower walking rate when fr > 0. 15. By fitting all data in Figure 12, the expression of β under the thermal transient condition can be given by 5 r 6 β R f R = + (20) where R5 and R6 are the fitting factors, given by Eq.17 Figure 11. Relative walking rate β against relative loading ratio f r under downslope condition. Figure 12 illustrates the responses of β in correlation with f r at different η for the thermal transient condition. It is clear that all responses of β at different η intersect at the axis of β = 0 when f r is equal to around 0.15. With increasing η, positive ∆ shows an increasing tendency when f r is lower than 0.15. The maximum β can be found at η = 1.5 and f r = 0.01. With increasing η, negative β shows a decreasing tendency when f r is lower than 0.15. The minimum β can be found at η = 1.5 and f r = 0.44. This implies that pipelines with larger gain in axial friction have a larger walking rate when f r < 0.15, and a lower walking rate when f r > 0. 15. By fitting all data in Figure 12, the expression of β under the thermal transient condition can be given by where R 5 and R 6 are the fitting factors, given by

Approach to Assessing Pipeline Walking Rate Considering Gain in Axial Friction
Based on the discussion above, this section proposes an approach by developing the traditional solution [1] to estimate the pipeline walking rate considering the gain in axial soil friction. The main calculation procedure during one arbitrary operational cycle can be outlined as shown in Figure 13.

Approach to Assessing Pipeline Walking Rate Considering Gain in Axial Friction
Based on the discussion above, this section proposes an approach by developing the traditional solution [1] to estimate the pipeline walking rate considering the gain in axial soil friction. The main calculation procedure during one arbitrary operational cycle can be outlined as shown in Figure 13.
Calculating the walking rate without considering the gain in axial friction can be referred to the analytical solution [1], where the detailed computation procedure is outlined for SCR tension, downslope, and thermal transient conditions.
(1) Calculate relative loading ratio f r in terms of individual walking-driven condition based on Equations (11)-(13). (2) Calculate relative pipeline walking rate β according to the corresponding walking-driven modes based on Equations (14)- (22). (3) Calculate the pipeline walking rate with the gain in axial friction taken into account S T by using Equation (2).
For the accumulated walking rates, a reasonable response of equivalent soil friction throughout the whole service life of the pipeline is required, and Equation (1) [15] is adopted. The total walking rates can be obtained by a simple summary of the walking rates in each operational loading cycle. The details will be described in the following section. Calculating the walking rate without considering the gain in axial friction can be referred to the analytical solution [1], where the detailed computation procedure is outlined for SCR tension, downslope, and thermal transient conditions.

1)
Calculate relative loading ratio fr in terms of individual walking-driven condition based on Equations (11)-(13).

3)
Calculate the pipeline walking rate with the gain in axial friction taken into account ST by using Equation (2).
For the accumulated walking rates, a reasonable response of equivalent soil friction throughout the whole service life of the pipeline is required, and Equation (1) [15] is adopted. The total walking rates can be obtained by a simple summary of the walking rates in each operational loading cycle. The details will be described in the following section. Relative loading factor f r = f 1 /f* Relative pipeline walking rate β Pipeline walking rate S R (considering gain in axial friction) According to walking-driven mechanisms Figure 13. Procedure of proposed method to assess pipeline walking rate.

Application Example Using Proposed Method
A case study of pipeline walking induced by the thermal transient is demonstrated in this section for further interpretation of calculation details and for comparison with the traditional method [1].
Taking an example of 1 km long pipeline, the relevant parameters necessary for pipeline walking design are illustrated in Table 3. Table 3. Parameters for the application study. First, pipeline walking rate S T during the individual operational loading cycle is calculated by using the analytical solution [1]. Table 4 illustrates the corresponding results of S T over 10 cycles. Then, the hardening model of axial soil friction [15] is adopted to determine the gain in axial friction for individual cycles. The detailed process is demonstrated in Figure 14 with 10 typical operational loading cycles. Note that two separate stepping increases in each cycle can be observed in the response of axial frictional coefficient (according to Equation (1)), because the pipeline would shear the ambient soil twice during the heating-up and cooling-down stages. Based on the response of µ N , the response of η for each operational cycle can be captured (see Figure 14). First, pipeline walking rate ST during the individual operational loading cycle is calculated by using the analytical solution [1]. Table 4 illustrates the corresponding results of ST over 10 cycles. Then, the hardening model of axial soil friction [15] is adopted to determine the gain in axial friction for individual cycles. The detailed process is demonstrated in Figure 14 with 10 typical operational loading cycles. Note that two separate stepping increases in each cycle can be observed in the response of axial frictional coefficient (according to Equation (1)), because the pipeline would shear the ambient soil twice during the heating-up and cooling-down stages. Based on the response of μN, the response of η for each operational cycle can be captured (see Figure 14).  Afterwards, using Equations (20)- (22), relative walking rate β can be calculated. The corresponding key parameters (η, R 5 , and R 6 ) and the results of β for 10 operational cycles are illustrated in Table 5. Finally, the walking rate considering gain in axial friction S R can be obtained using Equation (2). Table 6 presents the values of S R for 10 individual operational cycles. In addition, the accumulated walking rate is the sum of all S R . Figure 15 shows a comparison of the accumulated walking rates using the conventional analytical method (without considering the gain in axial soil friction), numerical modelling, and the proposed method, with the accumulative discrepancies in estimation included. As is evident, the numerical method and the proposed method in this study give almost identical estimations of walking rates. Compared to the conventional analytical method, with increased cycles, the accumulated discrepancies get larger, with the maximum discrepancy of around 6 mm at the end of 10 cycles. This is consistent with the changes of β illustrated in Table 5, where all values of ∆ are positive, from the largest value of 9.3% in the first cycle to the lowest value of 2.9% in the tenth cycle. After 10 operational cycles, the method proposed in this study gives a larger estimation of accumulated walking rate, with the cumulative difference between the two methods reaching 6 mm.  Figure 15. Comparison of accumulated walking rates with and without considering the gain in axial friction.

Conclusion
This paper reports an investigation of the pipeline walking rate considering the gain in axial friction by means of finite element modelling. Through the analysis, the main findings can be concluded as follows: 1) The pipeline walking rate is significantly influenced by gain in axial friction η and relative loading ratio fr, which are examined through comprehensive parametric study with 11 key influential factors taken into account.

Conclusions
This paper reports an investigation of the pipeline walking rate considering the gain in axial friction by means of finite element modelling. Through the analysis, the main findings can be concluded as follows: (1) The pipeline walking rate is significantly influenced by gain in axial friction η and relative loading ratio f r , which are examined through comprehensive parametric study with 11 key influential factors taken into account. (2) The performance of pipeline walking differs with walking-driven conditions. Under both SCR tension and downslope conditions, relative walking rate β always has a negative performance, implying that the total pipeline walking would be lower than that estimated by the conventional method [1]. Therefore, techniques to mitigate pipeline walking (e.g., suction caisson along the flowline) can be designed with smaller geometry in pursuit of lower cost with a sufficient margin of safety. However, under the thermal transient condition, relative walking rate β is governed by both η and f r , with a positive increase in β as η increases at lower f r , and a negative decrease in β as η increases at larger f r . As is evident, the conventional method may underestimate the pipeline walking rate and result in potential risk in routine design. Therefore, all cases with low relative loading ratios leading to positive β should be considered carefully. (3) For the accumulated pipeline walking rate through the whole service life of the pipeline, a reasonable and realistic hardening model representing the gain in axial friction due to the effect of consolidation is necessary, through the exponential hardening model [15] used in this study; other soil hardening models available to characterize the real field condition are also compatible with the framework proposed in this study.
Overall, this paper provides some fundamental understanding and a practical framework to extend the current solution to address the pipeline walking rate, but for the best estimation, more work (i.e., elasto-plastic mobilization in friction) will be undertaken in future research.

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

A
Cross-sectional area of pipeline E Young's modulus F SCR Tension force introduced by SCR f 1 Axial soil resistance acting on unit length of pipeline in heating-up process f r Relative loading ratio f * Thermal load acting on the unit length of pipeline G Gravity of the pipeline G 1 Component of pipeline gravity perpendicular to slope seabed G 2 Component of pipeline gravity tangent to slope seabed L Total length of the whole pipeline N Shearing cycles R Parameter in critical state soil model R 1 Fitting factor under SCR tension condition R 2 Fitting factor under SCR tension condition R 3 Fitting factor under downslope condition R 4 Fitting factor under downslope condition R 5 Fitting factor under thermal transient condition R 6 Fitting factor under thermal transient condition S T Pipeline walking rate without considering gain in axial soil resistance S R Walking rate considering friction gain q Temperature gradient in heating process W Submerged weight of pipeline per unit length x A Distance the virtual anchor point moves forward in a heating step x θk Heating length of pipeline α Thermal expansion coefficient β Walking rate index ∆T Temperature difference ∆s 1 Selected EAF during heating-up phase ∆s 1 Selected EAF during cooling-down phase ∆s 2 Selected thermal load during heating-up phase ∆s 2 Selected thermal load during cooling-down phase η Growth rate of equivalent frictional coefficient κ/λ Elasto-plastic volumetric stiffness ratio µ 1 Equivalent initial soil friction coefficient µ N Equivalent axial frictional coefficient in the nth shearing cycle φ Seabed slope