Correlation Analysis of Influencing Factors and the Monthly Displacement Increment of a Hydrodynamic Landslide Using a Pseudo-Maximum-Likelihood-Estimation-Mixed-Copula Approach

Heavy rainfall and changes in the water levels of reservoirs directly affect the degree of landslide disasters in major hydropower project reservoir areas. Correlation analyses of rainfall- and water-level fluctuations with landslide displacement changes can provide a scientific basis for the prevention and early warning of landslide disasters in reservoir areas. Because of the shortcomings of the traditional correlation analysis based on linear assumptions, this study proposed the use of a pseudo-maximum-likelihood-estimation-mixed-Copula (MLE-M-Copula) method instead of linear assumptions. We used the Bazimen landslide in the Three Gorges Reservoir Area as a case study to carry out the correlation analysis of the rainfall, water-level fluctuations, and landslide displacement. First, we selected several appropriate influencing factors to construct four candidate Copula models and estimated the parameters using the pseudo-MLE method. After the goodness-of-fit test, we selected the M-Copula model as the optimal model and used this model to study correlations between the monthly displacement increment of the landslide and influencing factors. We then established the joint distribution functions of these correlations. We computed and analyzed the overall and tail correlations between the displacement increment and the influencing factors, and we constructed the conditional probability distribution of the monthly displacement increment for different given conditions. The results showed that the pseudo-MLE-M-Copula method effectively quantified the correlation between the rainfall, reservoir-level fluctuations, and landslide displacement changes, and we obtained the return periods and value at risk of the influencing factors of the Bazimen landslide under different rainfall conditions and reservoir-level changes. Furthermore, the tail correlations between the monthly displacement increment of the landslide and the rainfall- and reservoir-level changes were higher than the overall correlations.


INTRODUCTION
The Three Gorges Hydropower Station is currently the largest hydropower station in the world, and it is the largest project ever constructed in China. Since the start of the experimental water storage in 2003, it has been in continuous operation for 16 years. Because of the complex geological conditions and the frequent human activities in TGRA, this area has been prone to frequent geological disasters over a long period of time (Chen et al., 2008;Yin et al., 2009;Kirschbaum et al., 2010;Miyagi et al., 2011;Ahmed, 2015). The ecological environment of the reservoir has deteriorated as a result of landslides caused by heavy rainfall and reservoir water-level changes, which has become a major potential hazard affecting the long-term operational safety of major hydropower projects and the ecological environment of the reservoir (Varnes, 1996). Important topics in the field of hydrodynamic landslide disaster monitoring and early warning and prevention include: 1) the analysis of the reservoir's hydrodynamic landslide deformation characteristics and the deformation-induced response, 2) the analysis of landslide influencing factors and deformation prediction, and 3) the establishment of a threshold extraction model. To reveal the disaster-causing mechanism of hydrodynamic landslides in reservoir areas, it is particularly important to study the correlation between the influencing factors and the displacement, which holds great significance for the monitoring and early warning of reservoir disasters and the prevention and mitigation of geological disasters.
Currently, correlation analysis of influencing factors of landslide deformation has been based primarily on linear theory. For instance, some traditional correlation analysis methods, such as linear correlation analysis, gray correlation analysis, and hypothesis testing, have been applied to select the influencing factors that affect the displacement of the hydrodynamic landslide (Ohtake, 1986;Kafri and Shapira, 1990;Wang et al., 2004;Pavan et al., 2012;Telesca et al., 2012;Miao et al., 2017a). In addition, some researchers have used autoregressive-moving-average (ARMA) time series analysis and Pearson product-moment correlation coefficients (PPCC) to analyze the lagged correlation between the monthly displacement increments and their influencing factors (Cai et al., 2016;Cao et al., 2016;Zhou et al., 2016;Li et al., 2018). Moreover, based on the daily rainfall data for the TGRA, the finite element method (FEM) and discrete element method (DEM) also have been used to simulate the effect of influencing factors on landslide stability and have verified the correlation between influencing factors and landslide displacement (Kawamoto, 2005;Lollino et al., 2010;Tang et al., 2017). Data-mining technique and clustering method have been applied to interpret landslide monitoring data and to select the influencing factors of the landslide deformation (Shiuan, 2012;Hong et al., 2016;Ma et al., 2016).
These research models and methods based on linear theory, however, can just depict the linear relationship between two variables, and they ignore the thick-tailed variables. Thus, they cannot satisfactorily represent the full relationship between the two variables, and it becomes extremely difficult to accurately measure the correlation structure of the two random variables when the distribution function of the random variables is uncertain or too complex (Iyengar, 1997).
Because the influences associated with landslide hazards are uncertain and dynamic (Bai et al., 2010), a more flexible and robust nonlinear correlation analysis tool (i.e., the Copula theory; (Sklar, 1959) has been used to analyze the correlation structure between variables when it is not certain that linear correlation coefficients can properly measure the correlation. In addition, the Copula method offers many advantages over traditional correlation analysis methods (Sklar, 1959).
Thus, the Copula model often is used to construct complex multidimensional probability distributions. For example, several kinds of Copula models were used to generate the joint probability density distributions of the rainfall, reservoir waterlevel changes, and landslide displacement to study their tail correlation (Motamedi and Liang, 2014;Bezak et al., 2016;Miao et al., 2017;Li et al., 2018;Wang et al., 2020). These methods, however, still have limitations (Hu, 2002). The Copula functions that belong to the elliptic Copula family in the high-dimensional construction case are extremely difficult, whereas construction of the Archimedean Copula function is simple. When the correlation between the random variables is too complex, using a single Archimedean or elliptic Copula function is more likely to produce a one-sided description. To solve this problem, the M-Copula method combing Frank, Clayton, and Gumbel Copula functions has been used for relationships and risk analysis between two random variables (Liu and Zhang, 2016). More important, nearly all of these methods have used the full maximum likelihood estimation (MLE) method to estimate the parameters of the Copula model. The MLE method, however, requires prior assumption of the marginal distribution of the random variables, which directly affects the results of the parameter estimation if the marginal distribution is incorrectly selected.
In summary, we proposed a new pseudo-MLE-M-Copula model to quantitatively investigate the correlation between the displacement of the Bazimen landslide and the influences of the rainfall and reservoir water-level fluctuations in the TGRA. This model not only overcomes the shortcomings of the linear method and more accurately describes the joint distribution relationship between the random variables, but also linearly combines the Gumbel Copula, Frank Copula, and Clayton Copula functions in the Archimedean Copula family to construct the M-copula function model. This model is more comprehensive and flexible, allowing it to describe complex correlations, and also adopts a semi-parametric estimation method. Thus, the pseudo-MLE, for the parameter estimation of the Copula model, eliminates the need to assume the marginal distributions and can effectively reduce the error of the parameter estimation. This research has provided the scientific basis for the study of reservoir hydrodynamic landslide disaster-causing mechanisms, landslide disaster monitoring and early warning, and disaster prevention and mitigation.
function. (Mackay, 1986;Nelsen, 1986). The Archimedean Copula function family has been used widely because of the simplicity of its computation process and the relative clarity of the construction model as well as because it is not limited by the positive or negative correlations between the different variables (Yan, 2006). Given random variables x and y, the twodimensional Archimedean Copula function is defined as follows: where φ(u) and φ −1 (u) are the element generators and their inverse functions for the Archimedean Copula, respectively; u 1 is the marginal distribution of variable x; and u 2 is the marginal distribution of variable y. The Gumbel, Clayton, and Frank Copula functions used in this study are shown below, respectively: In Eqs 2-4, C 1 (u 1 , u 2 ), C 2 (u 1 , u 2 ), and C 3 (u 1 , u 2 ) are the Gumbel, Clayton, and Frank Copula functions, respectively; and α, θ, and c are their corresponding parameters.
In 2002, Hu proved that if these three Archimedean Copula functions were combined linearly by adding coefficients (Hu, 2002), a mixed Copula function could be formed. The mixed function, which has all of the good features of the three Copula functions, can be described from several different angles to depict the more complex correlations between the variables. The expression is as follows: where C 1 , C 2 , and C 3 are the Gumbel, Clayton, and Frank Copula functions; m 1 , m 2 , and m 3 are the weighting coefficients of the three Copula functions, respectively; and m 1 , m 2 , and m 3 are ≥0, and m 1 + m 2 + m 3 1. The tail correlation coefficient (Zhang and Singh, 2007a) includes the lower-tail correlation coefficient and the uppertail correlation coefficient. The lower-tail correlation coefficient represents the effect on another random variable when one random variable takes a smaller value and is given by the following formula: The upper-tail correlation coefficient represents the effect on another random variable when one variable takes a larger value, and it is given by the following formula: where u ∈ [0, 1]; and C 1 − 2u + C(u, u) is the Copula survival function. If λ lo and λ up exist and both are in the interval [0, 1], the larger the value, the stronger the correlation. If λ lo and λ up are equal to 0, the random variables in the upper and lower tails are independent of each other. To further study the influence of each influencing factor on the landslide displacement deformation, we constructed a value-atrisk (VaR) model (Rockafellar and Uryasev, 2002) and tailed value-at-risk (TVaR) based on the Copula parameters and extracted the thresholds. The expression is as follows: where x is the value of each influencing factor; α is the confidence level (0.95 in this study); and F − 1 (x) is the cumulative probability density function of the selected optimal Copula function.

Data Preprocessing
To analyze the correlation between the rainfall, the reservoir water-level fluctuations, and the landslide deformation, we selected the optimal Copula model. From the type of data source, we knew that the monthly displacement increment, the rainfall, the two-month rainfall, the one-month reservoir waterlevel change, and the two-month reservoir water-level change were all discrete data. Thus, a discrete data transformation should be considered when constructing the binary Copula model. In this study, we used the distribution function transformation method proposed by Wu (Wu et al., 2009), which could solve the problem faster. The core idea of this method was that the random variable X (discrete or continuous) was known, and for any real number in x ∈ R, k ∈ I, and the defining functions are as follows: where k obeys a uniform distribution and is independent of X. From this, the generalized distribution function transformation of X is obtained: In this study, we used this method to carry out the generalized transformation of the four impact factor data sets X and the monthly displacement increment data Y: Q y D y, k P Y < y + kP Y y .

Parameter Estimation
The process of estimating the parameters of these four Copula functions using the MLE method (Dempster et al., 1977;Liu and Zhang, 2016) is as follows: Assume that the joint distribution of two-dimensional random variables is H M , and the marginal distribution is continuous; they are F 1 , F 2 respectively. According to Sklar's theorem, where C j is the Archimedes Copula function family with parameter β j , and C j (F 1 , F 2 ; β i ) is one single Archimedean Copula model, 3 j 1 m j C j (F 1 , F 2 ; β i ) is the expression of the M-Copula model. The main difference between MLE and pseudo-MLE is that the marginal distribution will be replaced with the empirical distribution function (Genest, 1987), and the expression of empirical distribution is as follows: So the joint density function with the empirical distribution is: The log likelihood function is: Solving the maximum point of the log-likelihood function, the pseudo-maximum likelihood estimates of m j and β j can be obtained:

Goodness-of-Fit Test
The pattern of the correlations between the multidimensional random variables can be described by the Copula model. The many types of Copula models reflect the different correlation patterns. To select the most appropriate Copula model for this study, we first selected the optimal model based on the Akaike information criterion (AIC), the Bayesian information criterion (BIC), and the root mean squared error (RMSE) values. Then, we used the Kolmogorov-Smirnov test (K-S test) to determine whether the optimal Copula model reflected the correct correlation structure between the variables.
The AIC is a judgment method based on the measurement of information, and it is suitable for testing the Copula models obtained using the MLE method. The BIC is more sensitive to models that are overestimated.
We used the RMSE to describe the error between the predictive Copula model and the empirical Copula model, as follows: where x est k represents the predicted estimated values, x emp k represents the empirical values, N is the capacity of the entire sample, and n is the number of parameters in the model. The K-S test (Christian and Bruno, 2008) can be used to test whether a sample obeys a specified distribution. We constructed the statistics of the K-S test according to the fact that the derivative of the Copula function obeys a uniform distribution within (0, 1). The original assumption H 0 was that a partial derivative of the selected Copula model obeys a uniform distribution within (0, 1).
where F n (x) is one of the partial derivatives of the Copula function; and F(x) is a uniform distribution within (0, 1).
Given a testing level α and a sample capacity n, the fractional table was used to find the critical value D n,α . If D max > D n,α , then H 0 is rejected; otherwise, H 0 is accepted.

Conditional Probability Distribution and Return Period
We further studied the joint distribution of the two random variables under different influencing factors (Zhang and Singh, 2007b). The conditional probability distribution was calculated using Eq. 23, and the periodic variations in the variables were described using the return period concept, which represents the number of time intervals per occurrence averaged over a certain  Frontiers in Earth Science | www.frontiersin.org April 2021 | Volume 9 | Article 637041 5 hydrological variable greater than or equal to a specified value (Eq. 24).

Building the Correlation Model
According to the field monitoring curves for the displacement of the Bazimen landslide and the reservoir's water level and rainfall, we found that the reservoir's water level, rainfall, and monthly landslide displacement increment exhibited different trends with time. First, the annual cycle of the reservoir water-level fluctuations and the rainfall exhibit a certain periodicity because of the existence of a rainy season in this region. Second, the monthly displacement increment of the reservoir landslide's deformation exhibited a weak periodicity under the effects of rainfall and reservoir water-level changes. To reveal the influences of the rainfall and reservoir water-level changes on the landslide deformation, in this study, we used the Copula model to assess the correlations between each influencing factor and the incremental monthly displacement of the landslide. The work flow of the correlation analysis of the landslide deformation and the influencing factors based on the Copula model is shown in Figure 1. It has two parts: model selection and model application and correlation analysis.
The specific implementation process is as follows. First, we selected the appropriate influencing factors to construct the Copula model, selected the optimal model based on the results of the goodness-of-fit test, established the joint probability density distribution functions for the landslide's monthly displacement increment and influencing factors, and analyzed the overall and tail correlations between the influencing factors and the landslide's monthly displacement increment. We constructed the conditional probability distribution model of the M-Copula function, and then calculated the conditional distribution of the displacement increment under different conditions. Accordingly, we obtained the return period of the displacement increment under the corresponding conditions. Finally, we extracted the thresholds of the rainfall and reservoir water-level changes in the reservoir area. This entire process plays an active role in the prevention and mitigation of hydrodynamic landslides in the TGRA.

CASE STUDY: BAZIMEN LANDSLIDE Data Source and Description
The Bazimen landslide is located on the right bank of the Xiangxi River in the TGRA (Figure 2). The slope is oriented north-south, and the landslide is spread in a winnowing fan shape at the foot of the slope. Its elevation is 139-280 m, and it is high in the west and low in the east. It slopes to the east, and the slope of the landslide is 10-30°, which is ladder-like and undulating. The part of the landslide above the water's surface is 380 m long, 100-500 m wide, and 10-35 m thick, and it has a volume of about 2 × 106 m 3 .
Since the Three Gorges Reservoir Area officially began storing water in 2009, the Bazimen landslide has experienced significant FIGURE 4 | Correspondence between the monthly displacement increment of the landslide and the reservoir water-level changes and rainfall at point ZG110. deformation with an increasing trend. The growth is characteristic of a typical step-type hydrodynamic landslide ( Figure 3). The causes of the landslide's formation and development are the topography, geology, geologic structure, other engineering geological conditions, rainfall, earthquakes, and reservoir water-level changes. The Bazimen landslide is strongly influenced by the dispatching of the Three Gorges Reservoir, and its cumulative displacement monitoring data have a strong correlation with reservoir water-level changes and the rainfall conditions ( Figure 3).
On the basis of this understanding, we developed a model to analyze the landslide's deformation triggers and displacementrelated structures using the rainfall and reservoir water-level changes as the basic triggering factors. The Bazimen landslide has multiple monitoring points, among which ZG110 and ZG111 have been monitored for the longest time and are located in the region of the landslide with the largest deformation. The monitoring data for points ZG110 and ZG111 are representative and most accurately represent the development and changes in the landslide displacement. In this study, we

Selection of the Influencing Factors
The curves illustrating the relationship between the landslide's monthly displacement increment, the reservoir's water-level changes, and the rainfall at monitoring point ZG110 on the Bazimen landslide are shown in Figure 4. As can be seen, at the end of each year when the reservoir's water level dropped, the landslide deformation rate was extremely high and underwent periodic changes. In addition, the influence of changes in the reservoir's water level on the landslide deformation was significant.
The effects of the rainfall-and water-level fluctuations on the landslide deformation were as follows. First, the weakening of the landslide's rocks caused by the rainfall infiltration and the reservoir water infiltration led to the deterioration of the mechanical properties of the landslide and resulted in the deformation of the landslide (Yange et al., 2015). Second, the rainfall infiltration and reservoir water-level fluctuations caused changes in the pressure of the seepage water within the landslide body, which led to the deformation of the landslide . Third, the periodic changes and time effects of the rainfall and reservoir water-level fluctuations also had significant impacts on the landslide deformation (Han et al., 2018). In general, because of the low permeability coefficient of the soil slope, the rainfall infiltration is relatively slow. The influence of longterm rainfall should be fully considered when selecting the influencing factor. We selected the rainfall in the reservoir area in the current month and two months as the influencing factor. The impact of reservoir water is usually a slow process. Therefore, based on the selection of the current monthly change of the current reservoir water level, we also used the two-month cumulative change of the reservoir water level.

Model Evaluation and Optimization
In this study, we selected the Gumbel, Clayton, and Frank Copula functions from the Archimedean Copula family and the M-copula function that linearly combines the three Copula functions as candidates.
First, it was necessary to select the Copula function that best reflects the trends and correlations in the raw data. The first step is the parameter estimation of the model. The Gumbel, Clayton, and Frank Copula models contained one unknown parameter each, and the M-Copula had three coefficients and three function parameters for a total of six parameters to be estimated. Many studies (Motamedi and Liang, 2014;Li et al., 2018;Wang et al., 2020) have used full maximum likelihood estimation (Nelsen, 2000) or the two-step maximum likelihood estimation (Genest et al., 1995) to estimate the model parameters. Both methods required a prior assumption of the type of marginal distribution, which directly affected the structure of the parameter estimation if the marginal distribution was incorrectly selected. In this study, we used the semi-parametric estimation method proposed by Genest (Genest, 1987) and Kim (Kim et al., 2007), that is, the pseudo-MLE method, to estimate the parameters of the Copula model. The parameter estimates are presented in Tables 2, 3. In all of the tables and figures that appear in this paper, the explanation of all the factors and variables are shown as follows ( Table 1).
The second step is to select the Copula model with the best fit for the four Copula models according to the AIC and BIC (Wang and Liu, 2006). The calculation results presented in Tables 3, 4 and 5 show that the M-Copula model had the best fit for the data source and thus was chosen for use in the following correlation study.
We used the K-S test to verify whether the partial derivatives of the M-Copula model obeyed a (0, 1) uniform distribution as a way to check whether the M-Copula model reflected the structure of the correlation between the variables well.
The original hypothesis was that the first-order partial derivatives of the M-Copula model would obey a uniform distribution within (0, 1). We tested the M-Copula function by constructing a K-S test statistic for the M-Copula model. As shown in Table 6, the p-value of the K-S test of the M-Copula model was greater than 0.05 for all six data sets, and the original hypothesis was accepted at the 95% confidence level. Thus, the first-order derivatives of the M-Copula model obeyed a uniform distribution within (0, 1), indicating that the M-Copula model reflected the correlation structure between the variables well.

Correlation Analysis Based on the M-Copula Model
We conducted a correlation analysis of the rainfall and reservoir water-level fluctuations and the landslide deformation using the M-Copula model. Figure 5 shows the spatial density distribution map and density contour metric of the M-Copula model. In all of the figures in this paper, the coordinate axis U represents each influencing factor and V represents the monthly displacement increment of the landslide. Table 7 presents the correlation metric parameters returned by the M-Copula function. Kendall's tau is the Kendall correlation coefficient and was used to show the overall correlation between two sets of random variables.  Table 6 means the test statistics of the K-S test of the M-Copula model. "P" in Table 6 means the p-value of the K-S test of the M-Copula model.
Frontiers in Earth Science | www.frontiersin.org April 2021 | Volume 9 | Article 637041 8 FIGURE 5 | (A), Spatial joint density distribution and contour plot between the 2 factors. "Reservoir water level rises in 1 month" and "Monthly displacement increment of ZG110"; (B), Spatial joint density distribution and contour plot between the 2 factors. "Reservoir water level drops in 1 month" and "Monthly displacement increment of ZG110"; (C), Spatial joint density distribution and contour plot between the 2 factors."Reservoir water level rises in 2 months" and "Monthly displacement increment of ZG110"; (D), Spatial joint density distribution and contour plot between the 2 factors. "Reservoir water level drops in 2 months" and "Monthly displacement increment of ZG110"; (E), Spatial joint density distribution and contour plot between the 2 factors. "Rainfall in 1 month" and "Monthly displacement increment of ZG110"; (F), Spatial joint density distribution and contour plot between the 2 factors. "Rainfall in 2 months" and "Monthly displacement increment of ZG110".
Frontiers in Earth Science | www.frontiersin.org April 2021 | Volume 9 | Article 637041 9 We used the upper-tail and lower-tail coefficients to measure the correlation between the two sets of random variables when the variable became larger or smaller.
The correlation parameters calculation results showed that the monthly displacement increment of the landslide was weakly negatively correlated with the one-month change Frontiers in Earth Science | www.frontiersin.org April 2021 | Volume 9 | Article 637041 10 in the reservoir's water level and the two-month change in the reservoir's water level, and it was weakly positively correlated with the other influencing factors. Note that the change in the reservoir's water level in the declining stage was more strongly correlated with the monthly displacement increment of the landslide. The upper-tail correlation coefficient of each group of variables was smaller than the lower-tail correlation coefficient, which is also shown in Figure 6. As the figure shows, the joint probability distribution was thicker in the lower tail and the lower-tail correlation was higher than the upper-tail correlation, indicating that the deformation of the landslide was likely to be smaller when the reservoir water-level fluctuations or rainfall were smaller.

Conditional Probability and Return Period of the Landslide Deformation
The conditional probability distribution diagrams and the return period diagrams of the random variables are shown in Figures 6,  7, respectively. When the amount by which the reservoir's water level changes was given, the conditional probability of the monthly displacement increment of the landslide gradually increased as the monthly displacement increment increased. When the amount of rainfall was given, the conditional probability of the monthly displacement increment gradually increased as the monthly displacement increment increased. In terms of the monthly displacement increment, the greater the change in the reservoir's water or the amount of rainfall, the smaller the probability value.
The return period of the monthly displacement increment of the landslide generally was characterized by a small return period when the change in the reservoir's water level was large and the monthly displacement increment was small. The return period was small when the amount of rainfall was large and the monthly displacement increment was small, and vice versa. When the reservoir water-level change or the rainfall was given, the return period became progressively larger as the monthly displacement increment increased. When the monthly displacement increment was given, it decreased as the reservoir-level change or rainfall level increased.  Frontiers in Earth Science | www.frontiersin.org April 2021 | Volume 9 | Article 637041 We arranged the Bazimen landslide displacement monitoring data in order of magnitude to create a new displacement sequence. The upper quartile of the displacement sequence of the monthly displacement increment of monitoring point ZG110 was 22.42 mm. If the monthly displacement increment was greater than 22.42 mm, we determined that the landslide had a large deformation in this month, and if it was less than 22.42 mm, we considered it to have a small deformation.
On the basis of this division, we reviewed the probability distribution and the return period of the monthly displacement increment for the selected −10 and 20 m monthly reservoir waterlevel changes, the −20 and 20 m two-month reservoir water-level changes, the 100 and 300 mm monthly rainfall levels, and the 200 and 500 mm two-month cumulative rainfall levels.
As shown in Figure 8A, the probability of a large landslide deformation in a month when the reservoir level dropped by 10 m was the highest (0.965) under various conditions. The probabilities for other changes in the reservoir's water level are shown in Table 8. Large landslide deformation was most likely to occur when the reservoir's water level fluctuated significantly, and the probability of a large landslide deformation occurring in a month when the reservoir's water level dropped was 0.965.
We obtained the following analysis results: When the reservoir's water level fluctuated significantly, the landslide was likely to experience deformation. In this case, the probability of a large landslide deformation occurring in a month in which the reservoir's water level dropped was higher than the probability of it occurring in a month in which the reservoir's water level rose. When the reservoir's water level remained unchanged, the probability of large landslide deformation occurring was significantly smaller than the probability of reservoir water storage or flooding. A reasonable explanation for this phenomenon was that most of the Bazimen landslides were below 145 m, and a rapid decrease in the reservoir's water level and the poor permeability of the landslide's geotechnical body caused the groundwater level within the slope to lag behind the reservoir's water level. This resulted in a positive difference between the groundwater level and the reservoir's water level. The outward infiltration of the groundwater from the landslide body caused the osmotic pressure to be directed toward the outside of the slope body. As a result, a significant deformation of the slope body was likely to occur during a sudden drop in the reservoir's water level.
As shown in Figure 8B, under different rainfall conditions in the reservoir area, the probability of a large landslide deformation occurring was the highest (0.900) when the cumulative rainfall for two months was 500 mm. The probabilities for the other rainfall conditions are shown in Table 9. This analysis led us to the conclusion that the greater the cumulative rainfall for one month or two months in the reservoir area, the more likely it was that the landslide would experience a significant deformation. This phenomenon may have been due to the fact that rainfall that infiltrated the slope increased the weight of the landslide and created pore penetration pressure. Thus, the slip zone soil also was softened by the water, which in turn would reduce its shear strength and cause a significant deformation.
By comparing the degrees of the influences of the reservoir water-level fluctuations and the rainfall on the landslide displacement in the reservoir, it was evident that changes in the reservoir's water level, especially when the reservoir's water level dropped, were more likely to cause large amounts of landslide deformation compared with rainfall for both the one-and two-month conditions.
We analyzed the pattern of the monthly displacement increment return period of the Bazimen landslide under specific conditions. As is shown in Figures 8C,E, the return period of the large landslide Frontiers in Earth Science | www.frontiersin.org April 2021 | Volume 9 | Article 637041 deformation (≥22.42 mm) when the reservoir's water level fluctuated greatly (up or down) was much smaller than the return period when the reservoir's water level remained unchanged or fluctuated less. Note that the return period for large landslide deformations was smaller when the reservoir's water level decreased than when the reservoir's water level increased. This result indicated that the probability of large landslide deformations was greatly increased in the case of sudden decreases in the reservoir's water level.
As the monthly displacement increment of the landslide increased from 22.42 mm to a maximum of 72.44 mm, the return period of the monthly displacement increment of the landslide increased by less than three months in the case of large fluctuations in the reservoir's water level. When the reservoir's water level remained basically unchanged, the return period increased by nearly 10 years. This result indicated that the landslide was more likely to undergo deformation destabilization and cause serious geological hazards when large fluctuations in the reservoir's water level occurred. Figures 8D,F, the return period of large landslide deformation (≥22.42 mm) was smaller when there was heavy rainfall in the reservoir area than when there was less rainfall. In particular, the return period of extremely large landslide deformation events (close to 75 mm) was much smaller for heavy rainfall than for smaller rainfall, with almost an order of magnitude difference. This indicated that the probability of large landslide deformation increased considerably under the heavy-rainfall scenario. As the monthly displacement increment of the landslide increased from 22.42 mm to a maximum of 72.44 mm, the return period of the monthly displacement increment of the landslide increased by only about seven months in the case of heavy rainfall. In the case of light rainfall, however, it increased by nearly four years. This result indicated that in the case of heavy rainfall, the landslide was also more prone to deformation destabilization and caused serious geological disasters. Table 10 shows the VaRs, TVaRs and their corresponding displacement increments of all the influencing factors based Frontiers in Earth Science | www.frontiersin.org April 2021 | Volume 9 | Article 637041 on the calculations using Eqs 8Eqs 9. On the whole, the displacement increment corresponding to the VaRs and TVaRs of rainfall is smaller than the displacement increment corresponding to the VaRs and TVaRs of reservoir level change, indicating that the response of landslide deformation to reservoir level change is higher than the response to rainfall. Moreover, the increment of displacement corresponding to VaRs and TVaRs when the reservoir level is falling is greater than that when it is rising. It can be assumed that landslide deformation responds more to a fall than to a rise in reservoir water level. When the four influencing factors reached or exceeded their respective VaR, the possibility of a large amount of landslide deformation occurring was greater. At that moment, it would be necessary to pay more attention to the landslide deformation and corresponding landslide geohazard emergency measures should be taken.

Case Validation: Baishuihe Landslide
The Baishuihe landslide is another typical landslide in the Three Gorges Reservoir area, located on the south bank of the Yangtze River in Zigui County, and its deformation development is mainly influenced by a combination of factors such as reservoir water level changes and rainfall. Monitoring data from the Baishuihe landslide from 2009 to 2018 were selected to verify the effectiveness of the proposed pseudo MLE-M-Copula approach. First, the monitoring data of monthly displacement increment, reservoir level changes and rainfall at monitoring point ZG118 were tested for normality using A-D test, J-B test and KS test. The p values of all variables were less than 0.05, indicating that none of the variables obeyed normal distribution and there is tail correlation between two variables. The M-Copula model proposed in this paper was used to fit the data, and the best fit was obtained on all candidate models based on three evaluation indicators: AIC, BIC and RMSE. The Kendall's τ and the tail correlation coefficients between each influencing factor and the monthly displacement increment were then calculated (Table 11), and it can be concluded that the monthly displacement increment is weakly negatively correlated with the reservoir level change (rise or fall) and weakly positively correlated with rainfall, which is basically the same as that of the Bazimen landslide. And the decrease in reservoir level has a stronger correlation with the incremental displacement. In addition, the upper tail correlation coefficient of variables is greater than the lower tail correlation coefficient. Unlike the Bazhimen landslide, the upper tail correlation of the displacement increment of the Baishuihe landslide with the influencing factors is higher than the lower tail correlation, indicating that when the reservoir water level changes or the rainfall is greater, the landslide is more likely to have large deformation. Finally, the VaRs and TVaRs of displacement increments and each influencing factor were calculated (Table 12), the results were roughly the same as the Bazhimen landslide, but the displacement increments were more responsive to rainfall. In summary, there are good reasons to believe that the M-Copula model based on pseudo-MLE proposed in this paper can effectively and accurately evaluate the correlation between the influencing factors and deformation of hydrodynamic landslides in the TGRA.

DISCUSSION
In this paper, we proposed a new M-Copula method based on pseudo-MLE to systematically investigate the correlation between landslide influencing factors and deformation of the Bazimen landslide in the TGRA. The method offered several advantages. First, the Copula method is based on nonlinear theory and can be used to completely characterize the correlation structure between multidimensional random variables. Second, the hybrid M-Copula method combined the advantages of three different Archimedean Copula functions to provide a more comprehensive and accurate description of the complex correlations. In addition, the pseudo-MLE method used to estimate the parameters of the model was able to avoid the adverse effects on the analysis results caused by the incorrect selection of the marginal distribution. Third, the introduction of conditional probability distributions and return periods allowed for a more accurate description of the effect of one variable on landslide deformation when one variable was fixed. This method also had limitations. First, the selection of influencing factors was largely dependent on empirical and qualitative analysis, and lacked concrete and computational support. Second, the study examined only the correlation between the monthly displacement increments of single monitoring point data and the influencing factor, and failed to consider the spatial characteristics of the landslide. Third, data from several monitoring points at different locations should be added to analyze the correlation structure between the three-dimensional spatial characteristics of the landslide and the landslide deformation, given that this study did not consider the changes in landslide displacement under the joint action of multiple influencing factors. The M-Copula model was firstly developed using a semiparametric estimation method (pseudo MLE) based on the monitoring data and compared with parametric and nonparametric estimation methods. In order to select the best-fit model on the monitoring data and to confirm its validity and accuracy, the M-Copula model was then compared with the Frank, Gumbel and Clayton copula models on the monitoring data. Three evaluation metrics, AIC, BIC and RMSE, were chosen to assess the goodness-of-fit, with the M-Copula model achieving the best results. These models were tested for statistical hypotheses using the A-D test, J-B test and KStest. Through correlation analysis, including overall Note: "+" means the reservoir water level is rising; "-" means the reservoir water level is falling.  correlation analysis and tail correlation analysis, it was found that the tail correlation was greater than the overall correlation and that the lower tail correlation was greater than the upper tail correlation. This result indicated that the possibility of a decrease in the displacement increment was significantly enhanced when there was a large decrease in the reservoir's water level or the amount of rainfall in the reservoir area. VaRs and TVaRs were used to calculate the threshold values for each influencing factor and its corresponding landslide displacement.
The results of the computational study validate the ability of the pseudo MLE-M-Copula model to analyze landslide deformation correlations and it can be well applied to other landslides in the same reservoir area. According to the results of the correlation calculations (tail correlation, conditional probability and return period, VaR and TVaR), there is a significant correlation between the landslide deformation, i.e., the landslide displacement increment, and the sudden drop of the reservoir water level, heavy rainfall. Therefore, the prevention of landslide hazards can be predicted by enhancing the monitoring of reservoir level changes and rainfall.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
1) The contribution of RW is the idea, introduction, theory and methodology, overview of the Bazimen Landslide, monitoring data analysis, conclusions. 2) The contribution of KZ is the relevance analysis.
3) The contribution of JQ is the programming of relevance analysis. 4) The contribution of WX is the displacement mechanism manof hydrodynamic landslide. 5) The contribution of Huanling Wang is the theory and methodology of Copula 5) The contribution of HH is the monitoring data analysis of Bazimen Landslide.