Reliability Modelling of Pipeline Failure under the Impact of Submarine Slides-Copula Method

: The instability of seabed slope sediments is the main factor inﬂuencing the safety of marine resource development. Therefore, to ensure the safe operation of submarine pipelines under complex and uncertain seabed rock and soil conditions, a reliability model was developed to elucidate the trend of impact-related pipeline damage due to submarine slides. Then, a risk assessment of the damage process of submarine slides impacting pipelines was conducted, which is of great signiﬁcance for the in-depth safety assessment of pipelines impacted by submarine slides. Based on the copula function, a joint probability distribution model considering the correlation among risk variables was established for rational correlation characterization. A probability analysis method of impact-related pipeline damage attributed to submarine slides based on the copula function was proposed. The Monte Carlo simulation (MCS) method was employed to simulate the random uncertainty in limited observation values and accurately determine the reliability of safe pipeline operation under the action of submarine slides. The conclusions were as follows: (1) Based on the copula function, a joint probability distribution model of risk variables with any marginal distribution function and related structure could be developed. (2) The copula function could reasonably characterize relevant nonnormal distribution characteristics of risk variables and could simulate samples conforming to the distribution pattern of the risk variables. (3) The failure probability calculated with the traditional independent normal distribution model was very low, which could result in a notable overestimation of the reliability of submarine pipelines.


Introduction
With the continuous progress in development and exploration technology, the exploitation of oil and natural gas has been increasingly promoted from land to sea. In recent years, the number of submarine pipelines has significantly increased, and the exploitation of offshore oil and gas resources has become a new field of oil and gas exploration worldwide. However, offshore oil and gas development faces more risks and challenges than those associated with onshore oil and gas development. Among these issues, the most important problem affecting the safety of marine resource exploitation is the instability of submarine slope sediments. Under the action of earthquakes and faults, gas hydrate disassociation, waves and currents, rock and soil masses, and sediments are susceptible to sliding, thereby forming submarine slides that could impact submarine pipelines [1,2]. Hance [3] noted that a submarine slide can be characterized by a large volume, large distance, and high speed. For example, the maximum value could reach approximately 20,331 km 3 , the maximum sliding distance could reach 850 km, and the maximum speed could reach 10,100 km/h. Submarine slides could easily impact submarine pipelines, destroy oil and gas transportation pipelines and exploitation facilities, and threaten the safety of offshore oil and gas development [4,5]. Therefore, to ensure the safe operation of submarine pipelines located on the ocean floor under highly complex and uncertain rock-soil mass conditions, a reliability model was constructed to elucidate the pattern of pipeline impact damage due to submarine slides, after which a risk assessment of the damage process of submarine slides impacting pipelines was conducted, which can provide an important theoretical basis for an in-depth safety evaluation of the effect of submarine slides on pipelines.
Since the 1970s, many scholars in China and abroad have studied the impact force of submarine slides exerted on pipelines, such as Demars [6], Randolph and Houlsby [7], and Zakeri and Hawlader [8], and various bearing capacity calculation equations have been proposed based on geotechnical mechanics theory. The impact force of a landslide was effectively combined with the rock and soil mass strength. Pazwash and Robertson [9], Chehata et al. [10], and Liu et al. [11] regarded the landslide mass as a fluid and analysed the impact force of submarine slides based on classical cylindrical fluid flow theory. Randolph and White [12], Dong [13], Dutta et al. [14], and Fan et al. [15] systematically analysed the impact force of submarine slides on pipelines by combining the theory of rock and soil mechanics with fluid mechanics theory. Although achievements have been made in recent years, the attained progress has promoted the development of submarine landslide-pipeline interaction research to a certain extent. However, the marine environment of submarine slides is complex and changeable. Compared to landsides, many uncertain factors exist, and the formation location is difficult to determine. In addition, field monitoring and sampling are extremely difficult and costly operations, and the obtained test data are limited, which further increases the uncertainty in seabed rock and soil mass parameters. Under the condition of limited data and a large number of uncertain factors, it remains difficult to apply the traditional analysis method to analyse the damage impact of submarine slides exerted on pipelines. Reliability theory can quantitatively consider multiple uncertain factors in a scientific and reasonable manner and can effectively overcome the limitation of a single index for structural safety evaluation. This approach has received increasing attention in the field of civil engineering and has been widely employed in structural safety design and analysis in bridge, structural, and other engineering fields [16][17][18][19]. However, in the marine engineering field, reliability research involving the impact of submarine slides on pipelines remains lacking. Therefore, the combination of reliability analysis theory and a method to evaluate the safety of pipelines under the damage impact of submarine slides could provide an important theoretical and scientific basis for disaster prevention and a mitigation design of marine energy exploitation systems.
The vertical force of submarine slides exerted on pipelines constitutes one of the important indexes used to evaluate the safety of submarine oil and gas pipelines. The vertical force fluctuates, which could pose a potential resonance risk to a given pipeline. When there exists a narrow gap between the pipeline and seabed, the load fluctuation magnitude is large, which matches that of the horizontal load. Moreover, the vertical force component could affect the pipeline vertical position, causing horizontal force fluctuations [20]. Fan et al. [21,22] found, through physical model tests and a large number of numerical calculations, that the impact force of submarine slides exerted on pipelines can be divided into two mechanical stages, including the instantaneous stage (the peak impact force is considered to represent the destructive effect of submarine slides) and the stable stage (the steady impact force is considered to represent the continuous effect of submarine slides on pipelines), and there exists a certain correlation between these two types of destructive forces. Through numerical analysis, a prediction model of the impact-related damage of pipelines due to submarine slides could be established, which could provide a reference for the safety and disaster prevention design of submarine pipelines.
However, the correlation between two or more risk factors has seldom been considered in numerical analysis methods, and risk variables are not simple and isolated quantities in complex geological environments. Moreover, if the correlation among risk variables is ignored and variables are assumed to follow independent normal distributions to simplify the research problem, the obtained reliability analysis results cannot truly reflect the impact of the correlation among the considered risk variables. Therefore, it is necessary to establish a reasonable joint probability distribution model of the risk variables for reliability analysis improvement. However, due to the high uncertainty in the submarine geological environment, there exists no unified opinion on the expression of the relationship between the peak impact force of submarine slides on pipelines and the steady impact force. At present, the primary task involves the urgent introduction of an accurate and reliable method to characterize the relationship between risk variables and construct the optimal joint probability distribution function in a scientific, reasonable, and comprehensive manner. In recent years, with the deepening of mathematical theory, the development of copula theory has provided a new method to establish joint probability distribution models of relevant nonnormal variables. The core idea of copula theory entails the separate construction of the marginal distribution function and copula function. There exists an unlimited variety of marginal distribution functions and corresponding structural types, and varied joint distribution models under arbitrary combinations can be established within the framework of copula theory. Copula theory has been widely applied in finance [23,24], hydrology [25,26], ecological sciences [27,28], reliability analysis [29,30], geotechnical engineering, and other fields [31][32][33][34][35][36] due to its incomparable flexibility and applicability in the establishment of a joint distribution model of variables. Currently, no study has fully considered the correlation among risk variables or has proposed a framework to analyse the reliability of safe pipeline operation under the action of submarine slides.
In summary, modelling was performed while considering the correlation among risk variables under the condition of incomplete probability information, the reliability of observations was evaluated through visualization methods, and the accuracy of the obtained reliability evaluation results was ensured. This paper fully considered two risk variables (peak and stable vertical force values) in evaluating the importance of submarine slides to pipeline safety. Based on the copula function, a joint probability distribution model considering the correlation among risk variables was established for reasonable correlation characterization. A probability analysis method for pipeline impact-related damage due to submarine slides based on the copula function was proposed. The Monte Carlo simulation (MCS) method was used to simulate the random uncertainty in limited observation values, and combined with big data analysis and visualization technology, the reliability of safe pipeline operation under the action of submarine slides was accurately analysed.

Joint Distribution Model of the Risk Variables Based on the Copula Function
In 1959, Sklar [37] first proposed the copula function and suggested that any multidimensional joint distribution function could be divided into a copula function and a corresponding number of marginal distribution functions. The copula function represents the correlation among variables (including the correlation coefficient and correlation structure). Its essence is the bridge function connecting the marginal distribution function of the variables with the joint distribution function and is often referred to as the combination function, bond function, or connection function. For a detailed theoretical introduction to copulas, please refer to Joe [38], Durante and Sempi [39], and Salvadori and De Michele [40].
The basic concept of the copula function is as follows: under N-dimensional conditions, the copula function can be defined as an N-dimensional joint distribution function with the marginal distribution function in [0,1] N space uniformly distributed in [0,1]. The Sklar theorem [37] can be described as follows: let F(x 1 , x 2 , . . . , x N ) be an N-dimensional joint distribution function with N marginal distribution functions F 1 (x 1 ), F 2 (x 2 ), . . . , F N (x N ). Then, there exists a copula function C(u 1 , u 2 , . . . , u N ) connecting the marginal distribution function F 1 (x 1 ), F 2 (x 2 ), . . . , F N (x N ) and the joint probability distribution function F(x 1 , x 2 , . . . , x N ).
Considering the Sklar theorem and a bivariate distribution [37], the joint distribution function comprises two parts: the distribution function of the variables and the copula function characterizing these variables. If H(x, y) is the joint distribution function with marginal distribution functions G(x) and Q(y), a copula function must exist. For any variables x and y, the copula function satisfies the following: where G(x) and Q(y) are the distribution functions of the risk variables, and θ denotes the parameters of the copula function.
If the marginal distributions G(x) and Q(y) are continuous, the copula function C is unique, and the joint probability density function can then be given as: where g X (x) and q Y (y) are the density functions of G(x) and Q(y), respectively, and D(G(x), Q(y); θ) is the density function of C(G(x), Q(y); θ).
The process of constructing a multivariate joint probability distribution function based on the copula function can be divided into two steps: (1) the marginal distribution function of the variables is determined based on the original data; and (2) the type of copula function properly representing the correlation among the variables is selected. These two steps are independent. Thus, the advantage of establishing a joint probability distribution model based on the copula function is that a normal or nonnormal distribution separates the marginal distribution from related structures, which can overcome the limitations of the traditional model. Moreover, a joint probability distribution model can be constructed with an arbitrary marginal distribution function and the related structure type of the joint probability distribution function, which can reveal the internal regularity of the original data.

Reliability Modelling Method of Pipeline Impact-Related Damage Due to Submarine Slides
To solve the problem whereby the traditional numerical analysis approach based on a single safety evaluation index cannot consider the limitations of various uncertain factors, and to reasonably characterize the correlation among the risk variables of pipeline impact-related damage due to submarine slides, a method based on the copula function was proposed to determine the correlation between submarine slides and pipeline damage risk variables and to evaluate the reliability of pipeline safety. Figure 1 shows a flowchart of the proposed method, which comprises three main steps, as described below.
Step 1: Determination of the optimal marginal distribution function of the risk variables.
The marginal distribution function can accurately describe the probability distribution of the variables. The primary task of the establishment of a joint probability distribution model based on the copula function entails the determination of the optimal marginal distribution function types of the variables. Since the destructive impact force of submarine slides exerted on pipelines is positive, this paper selected five marginal distribution functions commonly considered in engineering, namely, the normal distribution, log-normal distribution, truncated extremum type I distribution, Weibull distribution, and gamma distribution. Table 1 lists the various probability density functions and cumulative distribution functions of the different distribution types, where µ is the mean value and σ is the standard deviation.

Distribution Type Probability Distribution Function Probability Density Function Note
Notes: Φ() denotes the standard cumulative distribution function, φ() is the probability density function of the normal distribution, and Γ() denotes the factorial.
In engineering, the Akaike information criterion (AIC) [41] and Bayesian information criterion (BIC) [42] are commonly adopted to determine the optimal marginal distribution function. These criteria require that the marginal distribution function with the lowest calculated AIC or BIC value is the optimal marginal distribution function of the fitting variable. The above two criteria are simple in principle, provide a suitable stability, can be easily implemented in calculations, are widely applied in engineering, and can facilitate accurate and reliable data fitting. Therefore, the optimal marginal distribution function of the risk variables was identified and determined with the above method, and the specific expressions are as follows: where p and q are distribution parameters associated with µ and σ, respectively, n is the total number of data, and k 1 denotes the number of parameters of the alternative distribution types. Moreover, f (x i ; p, q) is the probability density function of the alternative distribution types. Table 1 indicates that all five distribution types contain two distribution parameters. Therefore, k 1 = 2, which is the minimum value considered in the calculation results, determines the optimal distribution type. Based on the above principles, the steps to determine the optimal marginal distribution function of the risk variables are as follows: (1) the mean and variance in the original data are calculated; (2) the parameters of the different types of marginal distribution functions are determined; (3) according to Equations (3) and (4), the optimal marginal distribution function of each risk variable is obtained based on the minimum calculated AIC or BIC value. When the alternative distribution types consider the same parameter samples, the identification results based on the AIC and BIC are the same. In this section, the AIC is used to identify the optimal edge distribution types.
Step 2: Determination of the optimal copula function fitting the correlation among the risk variables.
The correlation between parameters can be captured with the correlation coefficient and correlation structure type. In terms of the correlation coefficient, the Pearson linear correlation coefficient and Kendall rank correlation coefficient are mainly adopted. The Pearson linear correlation coefficient is an index used to measure the degree of linear correlation between the considered parameters. The Kendall rank correlation coefficient is based on the rank of the original parameter data and can describe the correlation between the parameters. The related structure types can be described according to the different copula functions. The θ parameter is the key to copula function determination and can be obtained based on the Pearson correlation coefficient and Kendall rank correlation coefficient θ [43]. According to the definition of the correlation coefficient, the relationship between parameter θ of the copula function and the Pearson correlation coefficient ρ is: where µ 1 and µ 2 are the average values of the two variables x 1 and x 2 , respectively, and σ 1 and σ 2 , respectively, are the standard deviations of these two variables.
With the above equation, parameter θ of the copula function can be determined. However, except for the Gaussian copula function, most copula functions are difficult to solve via integration. According to a previously reported method in the literature [44,45], parameter θ of the Gaussian copula function can be computed through the Pearson correlation coefficient, as follows: where ζ 1 = Φ −1 (u 1 ) and ζ 2 = Φ −1 (u 2 ) are variables of the standard normal distribution. Φ() is the standard normal distribution function, and Φ −1 () is the inverse of the standard normal distribution function. After parameter θ of the Gaussian copula function has been obtained, the Kendall rank correlation coefficient τ can be calculated with the following equation: Finally, the following equation can be employed to determine the parameter θ values of the different copula functions: In copula theory, many copula functions are available [43] that can be employed to describe variable correlation structure types, but most types can only facilitate simulations within a limited range of correlation coefficient values [46], e.g., elliptic copula functions such as the Gaussian copula function and t copula function, the Plackett copula function, and Archimedean copula functions, such as the Frank, Clayton, CClayton, Gumble, No. 16 and No. 17 copula functions.
The copula function couples the marginal distribution function and joint distribution function. After determination of the optimal marginal distribution function, the optimal copula function can be obtained based on the AIC and BIC. The corresponding calculation methods are expressed as Equations (9) and (10), respectively. Considering the above principles, the process of optimal copula function fitting is as follows: (1) the Pearson correlation coefficient and Kendall rank correlation coefficient are calculated; (2) parameter θ of the copula function is determined according to Equation (8); (3) AIC and BIC values can be computed with Equations (9) and (10), respectively, to determine the optimal copula function.
where (u i , v i ), i = 1, 2, . . . , n denotes the test data of the parameters, n is the total number of data, D(u i , v i ; θ) is the probability density function of the alternative copula function, k 2 is the number of parameters of the alternative copula function, and the minimum value in the calculation results determines the optimal copula function. The evaluation criterion is the same as that in step 2, and the AIC can be used to determine the optimal copula function.
Step 3: Reliability analysis based on simulations.
In practical engineering, due to the limitations of engineering technology and economic conditions, the available test and measurement data are very limited, and the joint probability distribution function of variables, which requires complete probability information, cannot be obtained. Only the marginal distribution function and correlation coefficient of the considered variables can be determined under limited data conditions, i.e., incomplete probability information. Especially in ocean engineering, it is very difficult to collect a large amount of high-quality test data. The MCS method can randomly generate sufficient samples based on the characteristics of the original data and has become a robust statistical tool [47]. This method is widely applied in probability analysis and provides an important technical means for the reliability analysis of submarine slide-impacted pipelines. Reliability analysis based on the MCS method directly calculates the failure probability by combining random simulation and statistical tests, and the calculation equation is as follows: where n f is the number of samples in the failure domain and n is the total number of samples. The method exhibits a clear concept, simple application, and few limitations.
With an increasing sample number, the calculation results become increasingly accurate and reliable. Under the condition of extremely limited data, the joint probability model of the considered risk variables was established based on the copula function, and relevant variables were effectively simulated with the MCS method. The corresponding probability was determined according to a large amount of simulation data. For example, Wang and Kulhawy [48] analysed the reliability of the normal service state and limit state of a given structure based on the MCS method. Aladejare and Wang [49] adopted the MCS method to generate a large amount of rock strength data and combined this method with correlation analysis to evaluate the reliability of rock slope stability. Pan et al. [50] examined the reliability of the tunnel driving face with the MCS method. Therefore, this study aimed to establish a joint probability model of the relevant risk variables based on the copula function and employed the MCS method to simulate a large amount of simulation data to approximate the actual variable distribution characteristics and process the limited available observation data.

Failure Probability Estimation under Submarine Slide-Induced Pipeline Impact Damage
In the structural system of civil engineering, Equations (12) and (13) can be applied to calculate the failure probability P f . It can be considered that Y(Q) denotes the minimum limit value, Y − Y(Q) < 0 is the failure condition, and the failure probability can be calculated with Equation (12). Alternatively, Y(Q) denotes the maximum limit value, and Y − Y(Q) > 0 is the failure condition. The failure probability can be calculated according to Equation (13).
where Y(Q) denotes the minimum limit value and Y is the actual observed value.
where Y(Q) is the maximum limit value and Y denotes the actual observed value. Compared to one-way flow through a cylinder, the process of submarine slide impacting pipelines is more complicated, and there are many risk factors. Determination of the risk level posed by submarine slides to pipelines enables the safety control of submarine pipelines. In this study, the computational fluid dynamics (CFD) method was adopted to simulate a series of multiphase flows, and peak and stable values of the vertical sliding pipe forces were determined under different landslide velocity conditions. Considering the peak and stable values of vertical forces, the failure probability of submarine slides impacting pipelines can be expressed as: where G(X 1 ) = X 1 − X 1 , G(X 2 ) = X 2 − X 2 , X 1 , and X 2 are the actual peak and stable values, respectively, and X 1 and X 2 are the maximum limits of the peak and stable values, respectively. For G(X 1 ) > 0 or G(X 2 ) > 0, the submarine pipeline is considered to be invalid. By establishing the joint probability distribution model of the peak and stable values considering different related structure types, the varied failure probabilities of pipelines under the impact of submarine slides can be determined.

Numerical Model
Induced by earthquakes, sedimentation, hydrate decomposition, and waves, a given submarine slope first becomes unstable, and the landslide mass begins to slide upon detachment from the unstable area. Under the influence of the water environment, the landslide mass gradually changes into a flow and continues to slide across a certain distance before stopping. In the whole landslide flow process, the flow velocity of the landslide mass gradually increases, reaching a maximum of up to 30 m/s [51]. Compared to the trigger start-up stage, the strength of the landslide mass at the flow slide stage is lower and the function rate is higher, which can significantly impact submarine pipelines. Therefore, based on the flow slip stage of submarine slides, this study employed the commercial CFD platform ANSYS-CFX to conduct numerical simulations of the impact of submarine slides on pipelines. The pre-and postprocessing tools required for CFD simulations include a 3D model, meshing scheme (ICEM-CFD), and a CFX solver based on the finite volume (FV) method. The CFD numerical method is very useful in fluid-structure interaction analysis. This study employed CFD software ANSYS 14.5 (CFX 2010a: CFX solver models; CFX program (version 13.0) physical modelling documentation, Canonsburg, PA, USA; ANSYS Inc. 2010b: CFX solver theory, CFX program (version 13.0) theory documentation, Canonsburg, PA, USA), which is a general-purpose CFD program including a solver based on the FV method for unstructured grids. A Euler-Euler multiphase flow model with nonuniform two-phase separation was applied to simulate the interaction between the submarine slide mass and seawater.
The established submarine slide-pipeline interaction numerical model is shown in Figure 2. The numerical calculation domain exhibits dimensions of 15.5 m × 8 m × 0.5 m (height × width × thickness), the pipeline diameter D pipe is 0.5 m, and the distance from the horizontal entrance is 2.5 m (5D pipe ). Moreover, the gap between the pipeline and seabed is H ps , and the landslide mass enters from a height of 10.5 m. The pipeline was considered fixed (or pipeline position variation could be ignored). The centre of the pipe is 2.75 m away from the inlet. The grid adopts tetrahedral elements, and the maximum mesh size depends on the pipe diameter. The maximum mesh size is 0.5D pipe , and the number of elements in numerical analysis exceeds 270,000. The grid in the area within a radius of 0.5-0.75 m (1.5D pipe ) was refined, and five layers of refined grids (with a total thickness of 0.05 m) were set up near the pipeline. The inlet was set as a velocity boundary, and the outlet was defined as an open boundary. The top and upper boundaries of the entrance were set as free-slip surfaces, and the surfaces of the pipeline and seabed were rough, each with an equivalent roughness k s of 0.0015 mm. The submarine slide flow was assumed to involve continuous free surface flow considering buoyancy and was simulated as an incompressible two-phase flow. All high-speed water and slide flow motions were determined based on the extended standard k-ε turbulence model. The landslide entrance was defined as the velocity boundary, the exit was set as an open boundary, the top of the computational domain was established as a free-slip boundary, the bottom and pipe surface were defined as rough no-slip boundaries, and the surface equivalent roughness k s values were 0.5 and 0.0015 mm, respectively. When the sliding distance reached 48 m (96D pipe ), the simulation calculation was terminated, and the calculation process adopted the second-order, high-precision upwind difference format. By varying the flow velocity and Reynolds number, the peak and stable vertical force values of 67 groups of submarine slides impacting pipelines were obtained.
Due to the complex marine environment and numerous factors influencing the stability of submarine slopes, submarine slides have become high-frequency geological disasters with a wide impact and potential threats. The highly notable impact produced seriously threatens the stability and safety of submarine pipelines. Therefore, it is of great significance to effectively simulate the impact of submarine slides on pipelines, reasonably determine the impact force of submarine slides exerted on pipelines, especially the vertical impact force must be improved, and study the correlation between the peak and stable vertical force values of submarine pipelines to accomplish more convincing safety and reliability evaluations of submarine pipeline projects.

Joint Distribution Model of the Risk Variables
Through calculation, 67 sets of data samples of peak and stable vertical force values were obtained (linear correlation coefficient R = 0.8004). Calculated peak and stable vertical force values over time, considering a flow rate of V = 1.0 m/s and Reynolds number of Re non-Newtonian = 8.7, are shown in Figure 3. Based on the original data scatter plot depicted in Figure 4a, the peak and stable vertical force values exhibited a lower tail correlation and were linearly positively correlated, but the linear relationship was not sufficiently obvious. The original sample data could be converted into uniformly distributed data with the semiparametric method based on maximum likelihood estimation, and the calculation process is expressed in Equation (14).
where (x i , y i ) is the original data sample value, Rank is a sorting function, which can be used to arrange the original sample data in ascending order, and (u i , v i ) is a standard uniformly distributed random variable after transformation. Figure 4b shows the standard uniformly distributed random variable after transformation. Compared to Figure 4a, Figure 4b reveals a more obvious linear positive correlation. Therefore, the candidate copula function selected in this study should be symmetric and must provide a good ability to describe the positive correlation structure of the random variables.
Many types of two-dimensional copula functions exist. The common two-dimensional copula functions can be divided into three types: (1) Gaussian copula functions; (2) twodimensional Plackett copula function; (3) two-dimensional Archimedean copula functions (for example, the Frank, Clayton, CClayton, No. 16. and No. 17 copula functions). To select the optimal copula function capturing the correlation among the risk variables, a copula function with a similar correlation structure to that of the measured data is usually selected in advance as the alternative copula function. Therefore, the copula function type selected in this paper can not only capture all copula function types, but can also capture top-and bottom-tail correlations to comprehensively analyse the original data and determine all possible correlation distribution types of the original data. The five considered alternative copula functions, i.e., the Gaussian, Plackett, Frank, CClayton, and No. 16 copula functions, are classic copula function families, and these functions can suitably describe both positive and negative correlations among variables. The absolute values of the correlation coefficients of these five copula functions all approached 1, which can meet the requirements depicted in Figure 4. Details of these five copula functions are summarized in Table 2.

Copula Function Type
Copula Distribution Function C(u 1 ,u 2 ;θ) Copula Density Function D(u 1 ,u 2 ;θ) ϕ θ (t,θ) In the third section, a method was introduced to determine the optimal marginal distribution function and copula function. Data statistics constituted the basis for the determination of the optimal marginal distribution function, and an important theoretical basis was provided to analyse data distribution characteristics. Table 3 lists statistical information on the peak and stable vertical force values of the 67 groups of submarine slides impacting pipelines. As indicated in Table 3, the mean value of the peak forces was 2.65 times that of the stable forces, and the standard deviation was 2.25 times that of the stable forces. The fluctuation range was larger than that of the stable forces, and the variation coefficient value of the peak forces was much lower than that of the stable forces. Therefore, the change in the stable forces was significantly greater than that in the peak forces.

Optimal Marginal Distribution Function
The optimal marginal distribution function was determined by comparing five marginal distribution functions (truncated normal, log-normal, truncated Gumbel, Weibull, and gamma distribution functions). Equations (3) and (4) were applied to determine the optimal marginal distribution function, and the results are provided in Table 4. Table 4 reveals that the optimal marginal distribution function for both peak and stable variable fitting was the Weibull distribution. To understand the fit between the original data sample and the marginal distribution function more intuitively, Figure 5 shows a histogram of the original data and the five considered marginal distribution functions. It is evident that the shape of the optimal marginal distribution function was better than that of the other marginal distribution functions, which indicates a better agreement with the distribution characteristics of the risk variables. The results in Figure 6 are consistent with those provided in Table 4, confirming the effectiveness of the AIC in determining the optimal marginal distribution function via fitting. To further verify the fitting effect of the marginal distribution function, the Kolmogorov-Smirnov (K-S) method was implemented to assess the fitting degree of the alternative marginal distribution function to the sample data. The K-S test is a probability distribution type test method suitable for small sample data sizes. By measuring the distance D between the known hypothesis probability distribution and the empirical distribution of the measured data, this method evaluates whether the distance occurs within the confidence interval [52]. The specific process of the K-S test method can be summarized as follows: let A 1 (x) denote the theoretical distribution function assumed in advance, while A 2 (x) denotes the actual cumulative distribution function of sample group A. Moreover, D is the maximum value of the gap between A 1 (x) and A 2 (x), i.e., D = max|A 1 (x) − A 2 (x)|. For D ≥ D n,α (D n,α is the rejection threshold), the original hypothesis can be rejected, and, conversely, the original hypothesis can be accepted. Table 5 summarizes the K-S test results for the risk variables (the peak and stable vertical force values). In regard to the risk variable of the peak vertical forces, the D value of the Weibull distribution is 0.0168, and, compared to the other marginal distributions, the D value is the smallest, i.e., the K-S distance is the smallest. According to the basic principle of the K-S test method, when the D value is smaller, it indicates that the two distributions are very similar and that the fitting degree is high, which further verifies the rationality of the Weibull distribution for peak variable fitting. Similarly, the Weibull distribution can be effectively used for stable variable fitting.    Figure 4b shows that the risk variables exhibited a significant positive correlation. If their correlation were ignored, a simple independent distribution could not represent the real distribution characteristics of the original data. Therefore, it is necessary to characterize the correlation among the risk variables based on the copula function and establish a joint probability distribution model. First, Kendall rank correlation coefficient values were calculated with Equation (8) to obtain the parameters of the copula functions with the different structures.
Then, the AIC or BIC was considered to determine the optimal copula function. The results are listed in Table 6. The table demonstrates that, among the five copula functions, the Plackett copula function yielded the lowest AIC values. Hence, the Plackett copula function could be effectively employed to fit the risk variable correlation structure. To further verify that the Plackett copula function is the optimal copula function, Figure 6 shows a frequency histogram of the original data sample, while Figure 7 shows a probability density plot and corresponding contour plot of the Plackett copula function. Figures 6 and 7a clearly exhibit the same shape overall. These two graphs indicate a trend of high values at both ends and low values in the middle, and both graphs exhibit symmetrical tails, suggesting that these graphs are sensitive to changes in the tail correlation between the random variables and reveal a high correlation between the variables. Hence, the symmetric tail correlation between the random variables can be better captured, thus confirming that the Plackett copula function can reasonably represent the structure of the correlation among the risk variables (peak and stable values). Moreover, Figure 7b shows that, in the contour map of the Plackett copula function, the risk variables exhibited a significant symmetry and a positive phase along the diagonal direction. Therefore, both Figures 6 and 7 verify that the Plackett copula function could reasonably represent the correlation characteristics of the original data and could be employed to establish a joint probability distribution model of the risk variables, thereby laying a foundation for subsequent reliability analysis.

Reliability Analysis
After determination of the optimal marginal distribution function and copula function, the joint probability density function of the risk variables can be obtained with Equations (1) and (2). To verify the superiority and importance of establishing a joint probability distribution model of the risk variables based on the copula function, the MCS-copula simulation method was applied to evaluate the reliability of pipelines under impact damage due to submarine slides. This method is helpful to accurately evaluate the safety status of submarine pipelines under the influence of landslide impact processes and provides an important theoretical basis for marine pipeline engineering design.
In this paper, a joint probability distribution model of the risk variables was established based on the copula function, and 10 5 samples were randomly generated with the MCS method, as shown in Figure 8. The figure shows that, without considering the correlation among the variables, the risk variables obey a normal distribution, whereas the simulated data exhibit a discrete uniform distribution, which is significantly different from the distribution characteristics of the original data. Moreover, due to the high variance in the original data sample and wide dispersion range, many negative values occur, which is inconsistent with the actual simulated data. The samples generated based on the above five alternative copula functions can describe the correlation among the variables. The simulated data were roughly distributed along the 45 • diagonal line, which is similar to the distribution characteristics of the original data. The simulated data were matched to the original data. Compared to the other four copula functions, the distribution characteristics of the simulated data obtained with the Plackett copula function were the closest to those of the original data, and this copula function could be adopted to accurately fit the distribution characteristics of the original data. This confirms that the Plackett copula function is the optimal copula function to fit the risk variables.
In this study, two risk variables (peak and stable vertical force values) were adopted as evaluation objects, and a series system [53,54] was employed as a criterion to evaluate the structural failure risk. In other words, when these two conditions were simultaneously satisfied, the structural system was considered to occur in the failure state. Figure 8 shows that, under a safety standard of 2500 N for the vertical peak forces and a stable force safety standard of 1000 N, the area simultaneously meeting these two risk standards is the safe area, indicated as the shaded green area, and the failure point occurs outside the shaded area. The failure probability is the ratio of the number of points in the failure area to the total number of points, as expressed in Equation (14).   The selection of the copula function directly determines the joint probability distribution model of the risk variables and thus greatly influences the structural reliability analysis results. Therefore, based on Section 5.2.2, and with the use of the Plackett copula function as an analysis standard, Table 7 lists the relative errors of the failure probability calculated with the different distribution models. Moreover, to analyse the change pattern of the risk variables and failure probability, this study simplified the problem and only considered the change in one variable. Figure 9 shows the change curve of the submarine pipeline failure probability under the influence of the risk variables. The results revealed that (1) it is unrea-sonable to employ the traditional independent normal distribution model, which ignores the correlation among the variables, and that the established joint probability distribution model could not provide accurate and reliable simulation data for reliability analysis. This could lead to inaccuracy of the calculated failure probability, and the error could even reach 82.54%. (2) The failure probability calculated based on the traditional independent normal distribution model was very low, which could result in a serious overestimation of the structural reliability and could bias the resultant structure design towards danger. The model established based on the copula function could accurately describe the correlation in the original data and could reasonably characterize the distribution characteristics of the original data, which could provide an accurate analysis model for structural reliability analysis and improve the calculation accuracy of the structural failure probability. (3) Due to the different correlation structures of the various copula functions, the calculated failure probability values notably differed. The failure probability error calculated with the No. 16 copula function was the largest, and the failure probability error calculated with the Gaussian copula function was the smallest. As shown in Figure 7, the different types of copula functions notably affected the distribution of the simulated data, which in turn influenced the reliability analysis results. (4) With increasing risk variable value (peak and stable vertical force values), the failure probability of submarine pipelines significantly increased. The five copula functions basically exhibited the same variation trend, but there occurred significant differences in the numerical values, among which the Gaussian copula function yielded the smallest difference.

Conclusions
This paper proposed a probability analysis method of pipeline failure under impact damage due to submarine slides based on the copula function. Under incomplete probability information, reasonable characterization of the correlation among the risk variables, improvement in the calculation precision of the pipeline failure probability under impact damage due to submarine slides, and the accurate assessment of the pipeline reliability under impact damage due to submarine slides are of great importance. The main conclusions were as follows: (1) Based on the copula function, a joint probability distribution model of the risk variables could be established given any marginal distribution function and related structure. The process of reliability analysis through a joint probability analysis model is as follows: a. the optimal marginal distribution function and optimal copula function were determined, and a joint probability distribution model was established and simulated in accordance with the distribution characteristics of the risk variables; b. based on the established joint probability distribution model, the MCS method was applied to generate a large number of random samples to calculate the failure probability of the pipeline impact damage attributed to submarine slides; (2) Under the condition of incomplete probability information, the copula function could reasonably represent the relevant nonnormal distribution characteristics of the risk variables, effectively establish a corresponding joint probability distribution model, simulate data conforming to the distribution pattern of the risk variables, and provide reliable and statistically significant samples for the reliability evaluation of the submarine slide effect on pipeline damage; (3) The traditional independent normal distribution model ignores the nonnormal distribution characteristics of the risk variables, and the calculated failure probability was very low, which could result in the serious overestimation of the reliability of submarine pipelines. Therefore, the correlation and nonnormal distribution characteristics of the risk variables should be comprehensively considered when evaluating the reliability of submarine pipelines.
The method proposed in this paper integrated a variety of uncertain factors, left the uncertain factors unrefined, established a joint probability distribution model of the distribution characteristics of the risk variables, and analysed the reliability of pipelines impacted by submarine slides. Based on existing research results, the first author will conduct a large number of model tests, study the correlation among dual-risk variables, and examine a large number of numerical simulations considering various uncertain factors of submarine slides to supplement the research results obtained in this paper.