Study on the streamflow compensation characteristics of reservoirs in Luanhe River Basin based on Copula function

With the rapid development of economy and society, the contradiction between supply and demand of water resources is increasing. Efficient utilization and allocation of limited water resources are one of the main means to solve the above contradictions. In this paper, the multidimensional joint distribution of natural streamflow series in reservoirs is constructed by introducing the mixed Copula function, and the probability of wet and dry encounters between natural streamflow is analyzed. Luan River is located in the northeastern part of Hebei Province, China, taking the group of Panjiakou Reservoir, Douhe Reservoir and Yuqiao Reservoir in the downstream of Luan River Basin as an example, the probabilities of synchronous and asynchronous abundance and depletion of inflow from the reservoirs are calculated. The results show that the probability of natural streamflow series between reservoirs is 61.14% for wetness and dryness asynchronous, which has certain mutual compensation ability. Therefore, it is necessary to minimize the risk of water supply security in Tianjin, Tangshan and other cities, and strengthen the optimal joint water supply scheduling of reservoirs. The research results are reasonable and reliable, which can provide reference for water supply operation of other basins.


INTRODUCTION
Study on the characteristics of wet and dry compensation in different regions is actually a joint probability distribution problem between variables. At present, the commonly used two-variable joint probability distribution models include two-variable normal distribution model, mixed Gumbel model, two-variable Gumbel-logistic model, exponential distribution model, Pearson type III model, Farlier-Gumbel-Morgenstern model, etc. These models are established based on the linear correlation between random variables, and the correlation between variables is measured by the linear correlation coefficient, which will lead to wrong conclusions when used to describe the nonlinear correlation problem. However, various random variables in the field of hydrology science often show complex linear and nonlinear correlations, so the above linear correlation variable distribution model cannot accurately describe the joint distribution between variables. Therefore, Copula function is introduced to describe the correlation between variables. Copula function can connect the joint distribution function and its respective edge distribution function. Based on the nonlinear correlation between variables, Copula function can describe the nonlinear, asymmetric and symmetric correlation between variables, which is widely used in hydrology.
Copula theory was first introduced by Sklar in 1959, and subsequently Nelsen provided a detailed introduction to Copula theory (Nelsen 2000). In recent years Copula theory has been widely used in the field of hydrology to solve the problem of multivariate joint distribution. Chen et al. (2012) established two four-dimensional Copula functions according to the joint distribution of flood magnitude and occurrence date, and calculated the coincidence probability of flood magnitude and occurrence date. Kuai et al. (2014) used Copula theory to construct the two-dimensional joint distribution function of the natural inflow of the two reservoirs in the whole year, flood season and non-flood season, and analyzed the streamflow encounter characteristics of the two reservoirs at three time scales. Yu et al. (2017) analyzed the influence of different Copula function parameter estimation methods on the selection of function types and the encounter of annual streamflow, which provided a theoretical basis for water resources development and water supply dispatching management in Xi' an. Sajjad Abdollahi et al. (2019) significantly optimized the modeling of complex hydrological phenomena using Copula function for multivariate probabilistic analysis of hydrological elements. Based on the long series of streamflow data of Dahuofang Reservoir, the Copula function and mathematical statistics method are used to analyze the streamflow compensation characteristics of Dahuofang Reservoir, and the diversion dispatching line of streamflow compensation conditions in wet, moderate and dry years is formulated (Zhang et al. 2020).

DISTRIBUTION MODEL BASED ON MIXED COPULA FUNCTION
Overview of the study area Panjiakou Reservoir is located at the junction of Tangshan City and Chengde City in Hebei Province. It is the source of the Luan River Diversion Project, which brings distributable water to Tianjin and Tangshan through the Luan River Diversion Project; Yuqiao Reservoir is an important hub project in Tianjin, which undertakes the task of water supply from Luan River Diversion Project to Tianjin; Douhe Reservoir is a comprehensive water conservancy hub for domestic water and industrial and agricultural production water in Tangshan City. In recent years, with the development of social economy, the shortage of water resources in Tianjin and Tangshan is becoming more and more serious, and the inflow of each reservoir shows a downward trend. Therefore, it has great significance to analyze the compensation characteristics between Panjiakou, Yuqiao and Douhe reservoirs for the joint scheduling of reservoirs and the rational allocation of water resources.
The joint scheduling between Panjiakou Reservoir, Daheiting Reservoir and Taolinkou Reservoir is mainly reflected in the agricultural irrigation water supply to the downstream of Luan River: (1) The Taolinkou Reservoir uses the remaining indicator water volume of Qinhuangdao City to increase the water supply to the downstream agriculture of the Luan River through the Qinglong River channel, and reduces the water supply to the downstream agriculture of the Luan River from the Panjiakou Reservoir and Daheiting Reservoir, thus increasing the water supply to Tianjin City or Tangshan City. (2) Using the remaining index water volume in Tianjin or Tangshan City, increase the agricultural water supply from Panjiakou and Daheiting Reservoirs to the lower reaches of the Luan River, The reduced agricultural water supply from Taolinkou Reservoir to the lower reaches of Luanhe River is supplied to Qinhuangdao City through the water diversion project from Qing to Qin.
Since Daheiting Reservoir is located 30 km downstream of the main dam of Panjiakou Reservoir on the main stream of the Luan River, it is close to the geographical location and climatic conditions, so the flow of the two reservoirs is basically the same abundance and depletion. Therefore, the compensation characteristics of inflow streamflow between Panjiakou Reservoir and Taolinkou Reservoir are analyzed.The hydraulic relationship between reservoirs is shown in Figure 2.

Copula theory and Copula function distribution model
In 1959, Copula function theory was proposed, Sklar suggested that a joint distribution function could be decomposed into a Copula function and K marginal distribution functions, and the decomposed Copula function could be used to describe the correlation between variables, that is, a function that could relate the joint distribution to its respective marginal distribution function. The N-element Copula function has the following properties: (1) Nondecreasing for any variable u n [ [0, 1], n ¼ 1, 2, . . . , N, C(u 1 , u 2 , . . . , u N ) (2) C(u 1 , u 2 , . . . , 0, . . . , u N ) ¼ 0, C(1, . . . , 1, u n, , 1, . . . , 1) ¼ u n (3) For any variable u n , v n [ [0, 1], n ¼ 1, 2, . . . , N: . . , N is independent of each other, and C ? is used to represent the Copula function of independent variables, then: There are many kinds of commonly used Copula functions, This paper mainly introduces multivariate normal Copula function, Archimedes Copula function and extreme value Copula function. Commonly used Archimedes Copula functions include Clayton Copula function, Gumbel Copula function and Frank Copula function.

Multivariate normal Copula function
The distribution function and density function of the N-element normal Copula function can be expressed as: where r is a symmetric positive definite matrix with element 1 on the diagonal; jrj represents the value of the determinant corresponding to matrix r; F r (Á, . . . , Á) denotes the standard normal distribution function with correlation coefficient matrix r; F À1 ( Á ) is the inverse function of standard multivariate normal function F( Á ); 6 ¼ (6 1 , 6 2 , . . . , 6 N ) 0 , 6 n ¼ F À1 (u n ), n ¼ 1, 2, . . . , N, I is a unit matrix. Because the multivariate normal Copula function has the characteristics of symmetry, it requires that the correlation structure between variables is the same, or the correlation coefficient between variables is very close, so it is difficult to fit the asymmetric correlation of variables.

Archimedes Copula function
The Archimedes Copula function can be expressed as: where w( Á ) is the generating element of Archimedes Copula function C(u 1 , u 2 , . . . , u N ), which needs to satisfy the following two conditions: Clayton Copula function, Gumbel Copula function and Frank Copula function are commonly used bivariate Archimedes Copula functions, which can be expanded as N-ary Archimedes Copula functions respectively: (1) The expression of N-element Gumbel Copula function is: (2) The expression of N-element Clayton Copula function is: (3) The expression of N-element Frank Copula function is: where a, u, l are relevant parameters.
The analysis of the correlation between variables by Gumbel Copula function, Clayton Copula function and Frank Copula function covers various situations of the change of correlation structure, which can basically comprehensively describe the complex correlation between variables. By calculating their relevant parameter values, it is easy to get the correlation measure values that people pay attention to, especially the tail correlation measure values, so they are most widely used in reality.

Extreme value Copula function
The expression of extreme value Copula(EVC) function is: Extreme Copula function is suitable for analyzing the correlation between variables caused by extreme events, so extreme Copula functions are increasingly used in the measurement of extreme event correlation.

CONSTRUCTION AND CORRELATION ANALYSIS OF MIXED COPULA FUNCTION Comparative analysis of different types of Copula functions
In general, multivariate normal Copula function usually describes the correlation between variables, but due to the symmetry of this function, it is difficult to fit the asymmetric correlation of variables.
The Gumbel Copula function and Clayton Copula function have the characteristics of asymmetry. The Gumbel Copula function is 'J' shaped distribution, the lower tail is low and the upper tail is high. It is not very sensitive to the change of the lower tail of the hydrological variables, but sensitive to the distribution of the upper tail. Therefore, it is difficult to describe the relevant changes of the lower tail. When one hydrological variable appears to be a maximum, the probability of the other two hydrological variables also appearing to be a maximum increases.
The Clayton Copula function is a 'L' shaped distribution, and the lower tail is high and the upper tail is low. It is sensitive to the change of the lower tail of the hydrological but not sensitive to the change of the upper tail. Therefore, it is difficult to describe the related changes of the upper tail. That is, when a hydrological variable appears to be a minimum, the probability of the other two hydrological variables also appearing to be a minimum increases.
Frank Copula function is a 'U' shape distribution, it has the characteristics of symmetry, so it is difficult to describe the asymmetric relationship between hydrological variables. The function is only suitable for describing the correlation between variables with symmetric correlation structure, that is, the maximum correlation and the minimum correlation between each hydrological variable are symmetric growth, but the distribution of the tail variable is relatively independent, so Frank Copula function is not sensitive in describing the correlation between the upper tail and the lower tail, so it cannot describe the correlation of the tail change.

Construction and correlation analysis of mixed Copula function
Gumbel Copula function, Clayton Copula function and Frank Copula function are three kinds of commonly used Archimedes functions that can capture the tail dependence: the correlation of the upper tail, the correlation of the lower tail, and the symmetric correlation of the upper tail and the lower tail. These Copula functions can describe the relationship between the various hydrological variables, especially the tail correlation. The relationship between the variables in the hydrological system is complex and changeable, not confined to a specific relationship. They can only reflect a certain aspect of the correlation between hydrological variables, so it is difficult to use a simple function to comprehensively describe the correlation pattern between variables in the hydrological system. Therefore, a more flexible function needs to be constructed to describe the relationship between various hydrological variable patterns.
The mixed Copula function M-Copula expression is: where C G , C F , C cl are Gumbel Copula, Frank Copula, Clayton Copula functions respectively; w G , w F , w cl are the weight coefficients of the corresponding Copula functions. From Equations (6)- (8), it can be seen that six parameters are included in MC 3 , parameter vector (a, l, u) is used to describe the degree of correlation between variables, the pattern of correlation between variables is represented by a linear weight parameter (w G , w F , w cl ) vector.

PARAMETER ESTIMATION, TEST AND EVALUATION OF COPULA MODEL Parameter estimation of Copula model
There are many kinds of parameter estimation in the Copula function model, and the maximum likelihood estimation and moment estimation are generally used. Among them, the maximum likelihood estimation is a commonly used parameter estimation method for the Copula model. The density function of its joint distribution function is: Among them: where u c is the 1 Â m c dimensional parameter vector of the Copula function; F n (x n ; u n ) is the edge distribution function; u n is a 1 Â m n -dimension parameter vector of edge distribution function F n (x n ; u n ); u ¼ (u 1 , u 2 , . . . , u N , u c ) 0 ; n ¼ 1, 2, . . . , N. Therefore, the log-likelihood function of the sample (x 1t , x 2t , . . . , x Nt ), t ¼ 1, 2, . . . , T can be obtained as: The u that makes the maximum value taken by the likelihood function is the maximum likelihood estimate.

Test and evaluation of Copula model
Whether the applied Copula distribution function can well fit the correlation structure and distribution between variables, so the test and goodness of fit evaluation of Copula function need to be established. Kolmogorov Smirnov (K-S) test is used to test whether the samples obey the same distribution, so it is used to test Copula distribution function; The evaluation of the fit of the Copula function is calculated using the root mean square error RSME minimum criterion, its definition is as formula (14): where N is sample size; i is the sample serial number; p c is the theoretical frequency calculated by the model; p 0 is the empirical frequency of joint distribution. Where the statistic D formula for the K-S test is shown in formula (15): where N is sample size; C k is the Copula value of sample x k ¼ (x 1k , x 2k , x 3k ); u k is the number of samples satisfying x x k condition, That is, satisfies :

CALCULATION METHOD OF FREQUENCY CURVE
The P-III curve is a commonly used frequency curve in current hydrological calculation, also known as gamma distribution. The probability density function is: where G(d) is the gamma distribution of d; d, a 0 , b are the shape, location and scale parameters of the P-III curve distribution, and d . 0, b . 0. Therefore, once the parameters d, a 0 , b are determined, the density function is also determined. It is concluded that these three parameters have the following relationship with the mean value x, skewness coefficient C s and variation coefficient C v of the three overall statistical parameters: The moment method is used to estimate the parameters of P-III frequency curve: where n is the sample size; k i is modulus ratio coefficient.

Distribution of inflow streamflow of each reservoir
According to the inflow streamflow data of Panjiakou, Douhe and Yuqiao reservoirs from 1956 to 2000, the annual average annual distribution of annual streamflow of each reservoir is analyzed, as shown in Figures 2-4.
From Figures 3-5, it can be seen that the streamflow distribution of Panjiakou, Douhe and Yuqiao reservoirs is decreasing year by year. Due to the rapid development of industry and agriculture, the water consumption increases annually, which makes the joint operation of reservoir water supply imminent.
The x, C s and C v parameters of streamflow samples from Panjiakou, Douhe and Yuqiao reservoirs are shown in Table 1.

The calculation, test and evaluation results of Copula model
According to the calculation results in Table 1 and formulas (11)-(15), each correlation parameter and test evaluation value of Gumbel-Copula, Clayton-Copula, Frank-Copula and M-Copula functions were calculated for Panjiakou, Steep River and Yuqiao reservoirs, respectively. The calculation, test and evaluation results of the function are shown in Table 2 (K-S test is significant at 0.05 level). Table 2 shows that Gumbel-Copula, Clayton-Copula, Frank-Copula and M-Copula functions can fit well the edge distribution of reservoir streamflow series by K-S test. According to the minimum criterion of root mean square error RSME, M-Copula function is selected as the connection function to fit the streamflow joint distribution of the three reservoirs.  The wet and dry indexes of the reservoir are divided into three levels : wet, moderate and dry, and the guarantee rate of the corresponding reservoir is divided by the transcendental probability. The wet and dry indexes are respectively p f ¼ 37:5% and p k ¼ 62:5%. When the M-Copula function is used to calculate the encounter of wet and dry, the above indexes are transformed into the cumulative probability, namely, P(X , x f ) ¼ 0:625, P(X , x k ) ¼ 0:375. The specific division of the three reservoirs' flow wet and dry index values are shown in Table 3.   Uncorrected Proof Based on the above classification criteria, the M-Copula connection function is used to analyze the streamflow series wet and dry compensation of Panjiakou, Douhe and Yuqiao Reservoirs. There are 27 kinds of wet and dry compensation encounter forms. Equations (21)-(23) lists three types of wet-dry encounter calculation formulas for Panjiakou, Douhe and Yuqiao reservoirs respectively, for wet-moderate-dry, dry-wet-wet and moderate-wet-moderate. The same applies to other cases, and the specific forms and calculation results of each are shown in Table 4.
Situation of wet-moderate-dry: Situation of dry-wet-wet: Situation of moderate-wet-moderate: From the analysis of 27 kinds of wet and dry encounters of streamflow series of three reservoirs in Table 4, it is known that : The probabilities of Panjiakou, Douhe and Yuqiao reservoirs being in the same wet, moderate and dry water years are 18.25%, 4.08% and 16.53%, respectively, which means the probability of synchronous of wetness-dryness is 38.86%, and the probability of asynchronous wetness-dryness is 61.14%. The high probability of asynchronous wetness-dryness relative to synchronous wetness-dryness indicates that the three reservoirs have certain mutual compensation ability, which provides more favorable conditions for joint water supply scheduling of the reservoirs.
When Panjiakou Reservoir is a wet year, and at least one of Douhe Reservoir and Yuqiao Reservoir is a combination of moderate or dry years, Panjiakou Reservoir has a certain compensation capacity for Douhe Reservoir and Yuqiao Reservoir, namely: wet-moderate-moderate, wet-moderate-dry, wet-moderate-wet, wet-wet-moderate, wet-wet-dry, wet-dry-moderate, wet-dry-wet, wet-dry-dry, its probability is 20.77%. When the inflow of Panjiakou Reservoir is a moderate year, and there is at least one combination of moderate or dry years in Douhe Reservoir and Yuqiao Reservoir, Panjiakou Reservoir has a certain compensation ability for Douhe Reservoir and Yuqiao Reservoir, that is: moderate-moderate-moderate,moderatemoderate-dry, moderate-moderate-wet, moderate-wet-moderate, moderate-wet-dry, moderate-dry-moderate, moderate-drywet, moderate-dry-dry, its probability is 19.1%.
When the inflow of Panjiakou Reservoir is dry year, and one of Douhe Reservoir and Yuqiao Reservoir appears dry year or two reservoirs are both dry year, Panjiakou Reservoir has no compensation ability for Douhe Reservoir and Yuqiao Reservoir. There are several combinations as follows: dry-dry-dry, dry-moderate-dry, dry-wet-dry, dry-dry-wet, dry-dry-moderate, its probability is 25.77%.