Dryness–Wetness Encounter Probabilities’ Analysis for Lake Ecological Water Replenishment Considering Non-Stationarity Effects

Ecological water replenishment (EWR) via interbasin water transfer projects has been regarded as a critical solution to reducing the risk of lake shrinkage and wetland degradation. The hydrological conditions of EWR water sources do not change synchronously, which may have an impact on the transferable water. Based on the GAMLSS model and the multivariate Copula model, this work presents a research approach for EWR via interbasin water transfer projects that can capture the non-stationarity of the runoff series and the frequency of dryness–wetness encounters, as well as speculates on various scenarios throughout the project operation phase. We present a case study on the Baiyangdian Lake, acting as the largest freshwater wetland in North China, which has suffered from severe degradation during the past decades and deserves thorough ecological restoration. The GAMLSS model was used to examine the non-stationarity characteristics of EWR water sources including the Danjiangkou Reservoir (DJK), the Huayuankou reach of the Yellow River (HYK), and upstream reservoirs (UR). The multivariate Copula model was implemented to evaluate the synchronous–asynchronous characteristics for hydrological probabilities for the multiple water sources. Results show that 1) significant non-stationarity has been detected for all water sources. Particularly, a significant decreasing trend has been found in UR and HYK. 2) The non-stationary model with time as the explanatory variable is more suitable for the runoff series of DJK, HYK, and UR. Under the non-stationary framework, the wet–dry classification of runoff series is completely changed. 3) Whether the bivariate or trivariate combination types, the asynchronous probability among the three water resources is over 0.6 except DJK-HYK, which indicates the complementary relationship. Multiple water resources are necessary for EWR. What is more, during a dry year of UR, the conditional probability that both DJK and HYK are in a dry year is 0.234. To alleviate the problem of not having enough water, some additional water resources and an acceptable EWR plan are required.


INTRODUCTION
With the rapid increase in the human population and social progress, development has generated significant economic and social benefits, but these benefits often come at high costs. The overexploitation of freshwater resources threatens the ecological environment and the overall well-being of humankind in many parts of the world (Kummu et al., 2016). At present, there are a large number of lakes shrinking and wetlands degrading around the world (Acreman et al., 2007;Chen et al., 2013;Mei et al., 2015;Liu et al., 2019;Ussenaliyeva and Aizhan, 2020;Jones and Fleck, 2020;Stone, 2021). In China, for example, 60.0% of lakes and 28.0% of marshes were threatened by the overuse of water resources, while 43.3% of lakes were threatened by sediments . Several wetland restoration measures, such as ecological water replenishment, restoring natural waterways, recreating the natural river, establishing habitats, and returning cropland to wetland, are recognized (Xin, 2014;Wu et al., 2020). Among these measures, the ecological water replenishment (EWR) has been widely applied in restoring ecology and hydrological conditions across various climatic and geophysical regions (Onuoha, 2008;Weigang et al., 2018), such as the Baiyangdian wetland (Ding et al., 2019), the Boluo Lake (Huang et al., 2021), and the Chagan Lake . Understanding the hydrological synchronization of numerous water sources is critical for the EWR, as the EWR measures are heavily reliant on the hydrological conditions of water sources. The water source area and intake area are usually geographically far apart, which leads to the temporal variation of runoff in the water source area and the intake area being not always synchronous. The asynchrony of the hydrological conditions bears directly on the transferable water quantity and the elapsed time for the water replenishment project. Furthermore, when water sources for EWR become diverse, the hydrological conditions of the water sources can hardly have a consistent change and made the EWR more complicated. Hence, a comprehensive analysis on the hydrological conditions of multiple water sources is suggested in aiding the formation of water diversion schemes Yan et al., 2018).
The basic assumption of traditional hydrological frequency analysis is the assumption of stationarity . However, hydrological series might become non-stationary due to changing climate and underlying surfaces, which leads to the invalidity of the results of hydrological frequency analysis (Lu et al., 2013;Jiang et al., 2015). Therefore, the frequency analysis of non-stationary hydrological series becomes a research focus in the past 20 years. The non-stationarity correction method (Ping et al., 2009) and the time-varying moment method Strupczewski and Kaczmarek, 2001) are widely used methods of non-stationary hydrological frequency analysis. The non-stationarity correction method attempted to accurately detect and decompose the abrupt and trend changes in hydrological time series and then compose these components (Xie et al., 2018a;2018b). The time-varying moment method assumes that hydrologic variables follow some distributions, including GEV and GAMLSS, in which one or more parameters are allowed to vary in time to reflect the nonstationarities of the hydrological series. One of the commonly used tools to examine the hydrological non-stationarity is the GAMLSS model Ahn and Palmer, 2016;Li et al., 2018). GAMLSS was proposed for fitting regression type models, where the distribution of the response variable does not have to belong to the exponential family, and includes highly skew, kurtotic continuous and discrete distribution (Rigby and Stasinopoulos, 2005). GAMLSS allows all the parameters of the distribution of the response variable to be modeled as linear/nonlinear or smooth functions of the explanatory variables (Stasinopoulos and Rigby, 2007). Because of these advantages, GAMLSS is widely used in the non-stationary hydrological frequency analysis (Villarini et al., 2009b;Yin et al., 2018;Zheng et al., 2018;Rashid and Beecham, 2019).
Univariate hydrological frequency analysis often fails to accurately characterize the hydrological conditions of multiple water sources. Therefore, copula functions were proposed (Sklar, 1959), which could quantify the dependence structure among correlated variables (Genest and Favre, 2007;Ariff et al., 2012), to determine the multivariate probability distribution. A copula is described as a function that links a multidimensional probability distribution function to its one-dimensional margins. Copula function is widely used in hydrology, including the joint frequency analysis of precipitation, drought, flood, and other extreme events (Grimaldi and Serinaldi, 2006;Shiau and Modarres, 2009;Xu et al., 2015). The studies mentioned above did not consider non-stationarity in the multivariate frequency analysis. In recent years, the non-stationarity in multivariate hydrological series has just begun to attract some attention only recently . Some studies have introduced the non-stationarity of marginal distribution into the joint frequency analysis based on the copula function (Kwon and Lall, 2016;Wu et al., 2020). Based on these, Jiang et al. believe that the changing environments have altered not only the statistical characteristics of some single random variables but also the dependence (i.e., statistical correlation) structure between different individual random variables. They employed a time-varying copula to analyze the effect of the time variation in the joint distribution on joint return periods of low flows . Some studies followed a similar method (Ahn and Palmer, 2016;Li et al., 2018;Vinnarasi and Dhanya, 2019;Wen et al., 2019).
However, the above studies about the non-stationarity in multivariate hydrological series focused on the joint probability distribution and return period of hydrological extreme events. Current studies mainly analyzed the dryness-wetness encounter probabilities of flow series under stationary conditions (Feng et al., 2010;Liu et al., 2015;Wang et al., 2017). The non-stationarity in multivariate series has not been fully considered when analyzing the dryness-wetness coupling of multiple water sources, particularly in formulating EWR schemes. Therefore, we recognize the necessity of a coupled analysis on the hydrological non-stationarity for multiple water sources in practice. The analysis of the encounter probability between water resources of water transfer projects enables decision-makers to make a reasonable lake water Frontiers in Environmental Science | www.frontiersin.org March 2022 | Volume 10 | Article 806794 replenishment plan. It motivates this study. In this study, we present a case study on the Baiyangdian Lake, acting as the largest freshwater wetland in North China, which has suffered from severe degradation during the past decades and deserves thorough ecological restoration. The GAMLSS model was employed to examine the hydrological non-stationarity of discharge for multiple water sources. Discharges considered in this study include those in the source regions of the South-to-North Water Diversion Project (SNWDP), the discharge at the Huayuankou section of the Yellow River, which represents the water diversion for ecological water replenishment in the Baiyangdian Lake from the Yellow River and the incoming discharge of the upstream reservoirs, that is, the Wangkuai Reservoir, the Xidayang Reservoir, and the Angezhuang Reservoir. After that, the copula function was employed to analyze the joint probability of the three water sources. Hence, the dryness-wetness encounter probabilities of the three water sources were calculated. The rest of this article is organized as follows. Section 2 introduces the main methods used in the present study briefly, including GAMLSS and copula. Details of the study area and data are presented in Section 3. Section 4 presents the parameter estimation and selection of marginal distribution and copula function along with the dryness-wetness encounter probabilities of the three water sources. Section 5 gives the conclusions of this study.

Trend and Change Points Analysis Methods
The trend component identification method used in this article is the Mann-Kendall trend test (Mccuen, 1994;Libiseller and Grimvall, 2010) and Sen's slope test (Mccuen, 1994;Gocic and Trajkovic, 2013). The two test methods are non-parametric tests, which make only mild assumptions about the data, and are appropriate when the distribution of the data is non-normal. Pettitt's test (Pettitt, 1979), the standard normal homogeneity test (SNHT) (Alexandersson, 1986), and Lanzante's test (Pettitt, 1979) are employed to identify change points of annual runoff time series. These methods are widely used, so the principles and calculation processes will not be reiterated here. The analyses were finished in R (package = "trend").

Marginal Distribution Using GAMLSS
To construct the dependence structure of hydrological variables by copulas, the marginal distribution of each variable should be determined first. Under the changing environments, the annual runoff series of many watersheds have been found to exhibit the so-called non-stationarity due to the effects of both climate change and human activates. As a result, the traditional method for runoff frequency analysis, which is based on the stationary assumption that the hydrological series should be independently and identically distributed, maybe no longer valid. Hence, we employed the time-varying moment model that expresses the distribution parameters as functions of time explanatory variable to capture the non-stationary characteristics of univariate runoff series by GAMLSS packages in R.
GAMLSS are univariate distributional regression models, where all the parameters of the assumed distribution for the response can be modeled as additive functions of the explanatory variables. GAMLSS provides over 100 continuous, discrete, and mixed distributions for modeling the response variable.
A GAMLSS model assumes independent observations y i for i 1, 2, . . . , n with probability (density) function f(y is a vector of 4 parameters, each of which is related to the explanatory variables. The first two population parameters μ and σ are usually characterized as location and scale parameters, respectively, while the remaining parameter(s), if any, are characterized as shape parameters, although the model may be applied more generally to the parameters of any population distribution.
where g k (·) is the known monotonic link function relating the distribution parameters to explanatory variables by semiparametric additive models. θ k and η k are vectors of length n. β T k {β 1k , β 2k , . . . , β J k k } is the parameter vector with sample size of J k ; X k is the matrix of the explanatory variable with length of n × J k ; Z jk is the fixed design matrix of n × q jk ; γ jk is the variable following the standard normal distribution. Without considering the impact of random effects on the distribution parameters, g k (θ k ) η k X k β k . For the location and scale parameters μ and σ, one or more frequently used variables need to be selected as the explanatory variables, including time (Villarini et al., 2009a(Villarini et al., , 2009bSun et al., 2018;Ting et al., 2018); climatic factors such as temperature, precipitation, and North Atlantic Oscillation (Villarini et al., 2012;Feng et al., 2020); factors relating to the underlying surface; and hydraulic engineering (Su and Chen, 2019;Wen et al., 2019). Previous research (Chen et al., 2007;Cong et al., 2009;Wang and Sun, 2021) has established that the runoff series of three study areas in this article exhibit significant non-stationarities with time, such as decreasing trends and change points, especially in the Yellow River Basin and the Haihe River Basin. Climate change and human activities have been identified as the two main reasons for the decrease in runoff. In the largest tributary of the Yellow River, the contributions of climate change and human activities to runoff decrease are 22-29% and 71-78%, respectively (Zhan et al., 2014). While in the Haihe River Basin, the contributions are 20-40% and 60-80% (Jia et al., 2012;Xu et al., 2014). The human activities here include land use, vegetation, and other land surface conditions, while climate change and climate variability are reflected in precipitation and temperature. It demonstrates that the decrease in runoff is the result of a mix of factors. Unfortunately, it is difficult to find appropriate indicators to fully characterize these factors in the absence of some key observation data, such as artificial water consumption. Hence, we take time T as the explanatory variable, which is effective and the most frequently used explanatory variable. A full-parameter model with time as the explanatory variable and without considering the impact of random effects is proposed as where μ and σ are vectors of length n, which indicate that the statistical parameters of non-stationary series change with time.

Copulas
Copula (Sklar, 1959) is described as a function that links a multidimensional probability distribution function to its onedimensional margins. The copula models are tools for studying the dependence structure of multivariate distributions. The usual joint distribution function comprises information on the marginal behavior of individual random variables as well as the dependency structure between the variables. The relationship between the correlated variables X i for i 1, 2, . . . , n is described by the following: where C is the distribution function and u i is the cumulative distribution of the X i variable and u i ∈ [0, 1]. A copula permits its marginal distributions to be evaluated by using different distributions. Among many families of copulas, the Archimedean copula has been most commonly applied in hydrology (Ahn and Palmer, 2016). Let function ϕ: [0, 1] → [0, ∞) be a strict Archimedean copula generator function and suppose its inverse ϕ −1 is completely monotonic on [0, ∞). A strict generator is a decreasing function ϕ: An Archimedean copula is defined as follows: The Archimedean copulas commonly applied are the Clayton copula, the Frank copula, and the Gumbel copula, as shown in Table 1.

Joint Probability Analysis Based on Copulas
A d-dimensional joint distribution probability can be defined as follows (Song, 2012): (6) where f X 1 X 2 ...X d (w 1 ,w 2 ,..,w d ) is the density function of X 1 , X 2 , . . . , X d . The marginal probability distributions of X 1 , X 2 , . . . , X d are denoted by Similar to the classification of wet-dry season, it is a common practice to classify runoff as wet, median, and dry periods according to the result of the frequency analysis in China to facilitate the management of watershed management. The common classification criterion is the cumulative probability distribution method. P(X < X i ) < P d : dry year, P d < P(X < X i ) < P w : median water year, and P(X < X i ) > P w : wet year.
where P d and P w are the thresholds of classification. P(X < X i ) is the cumulative probability distribution of runoff value X. Combined with the copula function, the dryness-wetness encounter probabilities of water replenishment resources are transformed into the multivariate joint probability.
For the bivariate joint probability, let X d and X w , Y d , and Y w be the thresholds of wet-dry classification of two runoff series, respectively. According to the probability inclusion-exclusion principle, the combination type can be calculated as follows: Synchronization probability is given as follows: Asynchronous probability is For the trivariate joint probability, there are 3 3 = 27 combination types based on the probability inclusion-exclusion principle. No more details are provided here, except the diagrammatic sketch ( Figure 1B), due to limited space.
Conditional probability is also employed in the analysis of the dryness-wetness encounter probability. Conditional probability is defined as the likelihood of an event or outcome occurring, based on the occurrence of a previous event.
Frontiers in Environmental Science | www.frontiersin.org March 2022 | Volume 10 | Article 806794 where P(A|B) is the probability of event A occurring given that event B has already occurred. P(AB) is the joint probability of events A and B. P(B) is the probability of event B.

STUDY AREA AND DATA
The BaiYangDian Lake (the BYD Lake) is located in the Haihe River Basin, which is the largest shallow lake/wetland in the North China plain. In the past 50 years, climate changes and human impacts have led to a sharp decrease in the amount of water entering the lake, and the water level in the lake has dropped significantly . According to the measured data of the Zaolinzhuang station at the outlet of the BYD Lake, the average water level was 7.67 m in 1950, and by 2018, it was only 6.88 m. The maximum water depth is 5~6m, and the average water depth is only 1~2 m. The decline of water level led to frequent drying up of the BYD Lake. In the 1970s, the lake dried up for 647 days in 1970-1973 and 1976. In the 1980s, the lake dried up for 6 years in a row from 1983 to 1988, for a total of 1845 days. In the 2000s, the lake dried up for a total of 1,488 days between 2000 and 2008. Low water-level events that occurred frequently and for a long time wreaked havoc on the lake's ecological health (Xu et al., 2011;Yang et al., 2016). The ecological restoration of the BYD Lake is nearing completion, and water replenishment via water transfer projects is one of the most efficient ways to meet the ecological water demand while also addressing the local water shortage. Three water sources of EWR are available. The upstream reservoirs (the Wangkuai Reservoir, the Xidayang Reservoir, and the Angezhuang Reservoir) are important water sources as local water resources. The total annual runoff of the three reservoirs is less than 0.1 billion m 3 . The targets are urban, agricultural, and ecological water consumption downstream of the BYD Basin. The locations of the three reservoirs and the BYD Lake are shown in Figures 2B. Except for the upstream reservoirs, two water transfer projects have been in operation for the EWR of the BYD Lake.

1) Danjiangkou Reservoir (DJK) in the Hanjiang River Basin.
The SNWDP is a large-scale water diversion project led by the Chinese government, which transfers water from DJK to northern China. The project was completed and put into operation in December 2014, and the BYD Lake received its first ecological water replenishment via the SNWDP in 2018.
The average annual inflow of the DJK reservoir is over 35 billion m 3 , and the average annual water diversion of the SNWDP is 8.54 billion m 3 . The main targets are 14 cities in northern China, including Beijing and Tianjin, that are suffering from water resources' shortage. It is worth mentioning that the SNWDP targets have substantial water use competition. The amount of water that can be used for ecological water replenishment is only 0.1-0.3 billion m 3 . 2) The mainstream of the Yellow River. The Yellow River Diversion Project is also a comprehensive large-scale water transfer project, which transfers water from the lower Yellow River to the BYD Lake. The targets of the project are the BYD Lake and five cities along the way. The project was completed and put into operation in November 2017. The annual average runoff at the Huayuankou section of the Yellow River (HYK) is 34.9 billion m 3 . HYK is a large hydrological station nearest to the water intake position of the project. The annual average water diversion of the project is 0.62 billion m 3 and the amount of water that can flow in the BYD Lake is about 0.1 billion m 3 .
In general, the annual average runoff is large, but the amount for the water transfer project is small and the amount flowing into the BYD Lake is tiny. The transferable water quantity of the water sources is limited during a dry year, and it is difficult to guarantee the water demand of the BYD Lake using a single water source. Therefore, EWR of the BYD Lake from multiple water sources and the consequent analyses of dryness-wetness encounter probabilities among multiple water sources are necessary and meaningful. The data used in this article include 1) the annual inflow of Danjiangkou Reservoir (DJK), 2) the annual runoff measuring at Huayuankou section of the Yellow River (HYK), and 3) the total annual inflow of the upstream three reservoirs, including the Wangkuai Reservoir, the Xidayang Reservoir, and the Angezhuang Reservoir (UR). The dataset with strict quality control covers 1961-2018 and is acquired from "Annual Hydrological Report P.R China" released by the Chinese government. It is worth noting that the EWR of the BYD Lake has only been operational since 2017. Long-term observation data after the EWR operation cannot be able to obtained. The nonstationarity of the dataset will be given full consideration in this study. What is more, the hydrometric stations observing this dataset are located upstream of the diversion ports (as shown in Figure 2), effectively avoiding the impact of water diversion projects' construction on data observation. Given the reasons, it is believed that the long-term datasets were of a good quality, which could fit the requirement of the dryness-wetness encounter probability analysis in this study.

Trend and Change Points Analysis
Hydrological series are generally composed of deterministic aperiodic components, deterministic periodic components, and stochastic components (Xinan and Ping, 2012). The deterministic aperiodic components include transient components such as trend and change point, which are often superimposed on other components. When the hydrological time series has a significant trend and change point, it is no longer consistent. The traditional frequency analysis is based on the consistency of hydrological series, so the trend and mutation point should be tested before fitting the marginal distribution.
The trend test results are shown in Table 2. The annual runoff series of DJK, HYK, and UR have decreasing trends. The annual changing rates in HYK and UR are −5.4042 and −0.1877 billion m 3 /10a, respectively. HYK's annual runoff was approximately 34.9 billion m 3 in 1961-2018 and 28 billion m 3 in 2009-2018, while that of UR was 0.959 and 0.590 billion m 3 , respectively. The result indicates that the decreasing trends in HYK and UR are significant (p value < 0.01). For DJK, the significance of the decreasing trend does not reach up to 0.05 level.
Pettitt, SNHT, and Mann-Kendall tests were employed to analyze the change point of the annual runoff series in DJK, HYK, and UR ( Table 3). The change point is 1990 in DJK. Some studies (Dongfei et al., 2016) suggested that the runoff decreased sharply because of the great climate change around 1990, such as the decrease of precipitation and the increase of temperature. The change points are 1985 and 1990 in HYK. Some studies (Jun et al., 2014) found that the atmospheric circulation in this area was abnormal from 1985 to 1990. The Mongolian low pressure significantly weakened and the summer monsoon was weak, resulting in a decrease in the precipitation and then in the runoff. The change point is 1996 in UR. Some studies (Lingling and Si-rui, 2016) believed that the abnormal decrease of runoff was caused by the change of the runoff generation mechanism. Linear regression curve and mean values before/ after change points are shown in Figure 3.
From the results of trend and change point analysis, it is clear that the runoff in the regions of the three water sources showed significant non-stationarity in the past decades, which is in accordance with the previous studies. The annual runoff (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018) of three hydrological stations in the middle reaches of the Yellow River decreased by 14.33%, 32.16%, 40.79%, and 44.32% compared with the maximum annual runoff (1959)(1960)(1961)(1962)(1963)(1964)(1965)(1966)(1967)(1968) at each station (Wang and Sun, 2021). The results of the other two water sources with decreased trends are also supported by previous studies (Chen et al., 2007;. In short, before conducting a hydrological frequency analysis, the non-stationarity of the runoff series must be considered in the planning and construction of water infrastructures including interbasin water transfer projects, reservoirs, and groundwater projects (Dong and Zhang, 2014;Isensee et al., 2021).

Marginal Distribution Fitting Based on the GAMLSS Model
Annual runoff series were non-stationary, according to the above results, due to change points and trends. Based on the result of the consistency test in Section 4.1, the GAMLSS model was employed to fit the marginal distributions of annual runoff in DJK, HY, and UR. Gumbel (GU), Weibull (WEI), gamma (GA), logistic (LO), log normal (LOGNO), and normal distribution (NO) were the candidate distributions. There were two types of fitting: (A) stationary marginal distributions and (B) time-varying marginal distributions with time as the explanatory variable. The position parameter μ and scale parameter σ were regarded as time-varying parameters to avoid overly complex regression equation, while shape parameter ν was assumed as a constant; that is, 1) both μ and σ are constants; 2) μ is a polynomial function of time t, while σ is a constant; 3) σ is a polynomial function of time t, while μ is a constant; and 4) both μ and σ are polynomial functions of time t. Only the linear and quadratic polynomial functions are included. The Akaike information criterion (AIC) was used to evaluate the goodness of fit of the distribution models and functions. The worm plot was used for visualizing the fitting performance of candidate models.  The complete results of marginal distribution fitting are given in Appendix Table A1. Table 4 shows the distribution name, distribution parameters, and AIC of the selected optimal marginal distribution. The marginal distribution fitting results of stationary type, namely, both μ and σ are constants, are also listed in the table for comparison. As shown in Table 4, for the annual runoff series with significant change points (like HJK and UR), non-stationary and stationary models had distinctly different AIC values. However, for the modeling of annual runoff series with the narrow changing range before and after their change points, non-stationary and stationary models had similar AIC values, which indicates that time-varying marginal distributions are suitable to capture the non-stationary characteristics. According to the AIC values, non-stationary models perform better than stationary models. The position parameter of the log-normal distribution for describing runoff series at DJK, HYK, and UR relates to the time (explanatory variable) negatively, whereas the scale parameter is constant. Therefore, the log-normal distribution with time-varying parameters is the selected optimal marginal distribution for the annual runoff series of DJK, HYK, and UR.
Analysis of residuals is necessary to evaluate the performance of the selected models. Figure 4 demonstrates the worm plots of the residuals by GAMLSS for the selected optimal marginal distributions in Table 4. It can be observed from the figure that the sample points of the annual runoff series follow the red solid curves fluctuating between two dashed curves, implying that the selected models have a quite good fitting quality at the 95% confidence level. Figure 5 shows the comparison of the cumulative probability distribution of the selected models. CDF_s and CDF_n mean the cumulative probability distribution function of the stationary model and the non-stationary model, respectively. The bars are the difference between the two (CDF_s-CDF_n). The green vertical lines are change points of the annual runoff series. The CDF of the selected non-stationary model is different from that of the stationary model. Before the change point, the CDF is greater than that of the stationary model (CDF_s > CDF_n), while it is reversed after the change point (CDF_s < CDF_n). Considering the decreasing trend of annual runoff in DJK, HYK, and UR, it is easy to be understood.
In China, it is a common practice to classify runoff as wet, median, or dry periods based on the results of frequency analysis to facilitate water resource management. Typically, P d and P w are 37.5% and 62.5%, respectively. According to the above, the choice of the stationary or non-stationary model will change the CDF of the annual runoff, which in turn changes classification results ( Figure 6). The transfer matrix from the stationary to nonstationary model in DJK indicated shows that 8 out of 58 years in the classification result changed. For HYK and UR, the amount of change years is 29 and 19, respectively. This also indicates that the traditional runoff classification without  considering non-stationary is not accurate enough. Based on this fact, we further calculated the dryness-wetness encounter probability of annual runoff in DJK, HYK, and UR, using the non-stationary model.

The Dryness-Wetness Encounter Probability Based on the Copula Function
Clayton copula, Gumbel-Hougaard copula, and Frank copula were selected as candidate distributions to construct bivariate and trivariate joint distribution functions. The maximum likelihood method is used to estimate the parameters of the copula function. The Akaike information criterion (AIC) and the Cramér-von Mises test (Genest et al., 2009) were used to evaluate the goodness of fit of the copula function. Table 5 summarizes the results of parameter estimation and goodness-of-fit testing for three candidate copulas. Bold text means the selected copulas in this study. Since the correlation between DJK and UR is negative, the Gumbel-Hougaard copula that allows only for positive dependence variables was excluded from the candidate distributions. As shown in Table 5, the selected optimal copulas for DJK-HYK, DJK-UR, and HYK-UR were Gumbel copula, Frank copula, and Clayton copula respectively. And Clayton copula was the optimal copula for DJK-HYK-UR.
According to the selected optimal copula and the parameter estimated by the maximum likelihood method, the bivariate and trivariate joint distributions are the probability of the annual runoff series. As mentioned earlier, P d and P w are 37.5% and 62.5%, respectively, in the runoff classification analysis. We classified runoff as wet, median, and dry periods according to Eq. 7, and then the encounter probability of bivariate combination types was calculated using the probability inclusion-exclusion principle ( Table 6).
For the encounter probability of bivariate combination types, the asynchronous probabilities of DJK-UR and HYK-UR are over 0.6. It seems reasonable as HYK is to the south of UR and DJK is  Frontiers in Environmental Science | www.frontiersin.org March 2022 | Volume 10 | Article 806794 9 more south geographically ( Figure 2). Different geographical locations lead to differences in factors such as climate and vegetation, and further lead to a higher asynchronous probability of runoff series. The result also indicates that the interbasin water transfer project is a reasonable way for EWR of the BYD Lake.
For the encounter probability of trivariate combination types, the result of 27 types is presented in Table 7. The synchronous probability of runoff series of three water sources, namely, the synchronous case (P S P(D&D&D) + P(M&M&M) + P(W&W&W)) is only 0.180. The asynchronous probability is higher than the bivariate combination types. But it is also noted that the asynchronous probability of DJK and HYK is 0.602, which indicates that there is no substitution between the two water sources. DJK and HYK are both important and necessary sources of EWP for the BYD Lake. A single water source is difficult to guarantee the ecological water demand of the BYD Lake. We suggest that the designed annual replenishment water amount of the BYD Lake should be adjusted to adapt to the wetness-dryness situation of DJK and HYK.
Because the BYD Lake and the three upstream reservoirs (UR) are in the same watershed, the dryness-wetness situation of UR represents, to some extent, the ecological water replenishment demand of the BYD Lake. Table 8 shows the result of conditional encounter probability. During a dry year of UR, the conditional probability that neither DJK nor HYK is in a dry year is 0.465. In this scenario, water from other watersheds partially meets BYR's water demand. But the conditional probability that both DJK and HYK are in a dry year is 0.234. Thus, to solve the problem of not having enough water, some additional water resources and a reasonable EWR plan are required.
The grim situation could be avoided by reservoir regulation, activating emergency water source and diversifying water sources (such as recycled water).   It is worth noting that 1) the probability that the water supply situation is tough is only calculated by the runoff series classification of each water source, not a direct comparison of water quantity. Perhaps when DJK or HYK is in a dry year, the amount of water that can be transferred into the BYD Lake is enough to meet the ecological water replenishment demand. In this case, the available amount of ecological water replenishment is only limited by the design amount of the water diversion projects. In other words, we only examined the water supply situation from the perspective of the water source, without considering the quantity of water diversion and water supply available from the water transfer projects. 2) DJK and UR are the reservoirs with regulating water storage capacity. The data used in this article are the inflow runoff series of reservoirs. The available water quantity of the water transfer projects is ultimately determined by the water level in the downstream area of the reservoirs.
Despite the limitations mentioned earlier, the results of encounter probability analysis could assist EWR project decision-makers in fully understanding potential challenges and devising appropriate countermeasures, and thus ensure the effectiveness of ecological replenishment projects. We believe that this article provides a research approach for EWR via interbasin water transfer projects that can capture the non-stationarity of runoff series and speculate on possible scenarios during the project operation phase. Current studies about EWR focus on targets' water demand and give simple constraints for water diversion and supplementation amount in their water replenishment optimization model (Huang et al., 2021). If more attention was paid to analyzing the non-stationarity and the dryness-wetness encounter probability, the optimization model would be more reasonable. Indeed, a large number of interbasin water transfer projects have been or are being built for the purpose of overcoming water scarcity and alleviating environmental problems (Guo et al., 2012;Akron et al., 2017;Roozbahani et al., 2020;Lei et al., 2021;Sun et al., 2021), which provides opportunities for the approach's application.

CONCLUSION
This article aimed to investigate the dryness-wetness encounter probability of runoff series between the lake replenishment water sources. The dryness-wetness bivariate and trivariate encounter probability analyses under non-stationary conditions is investigated in this study. Non-stationary frequency analysis is modeled with time as the explanatory variable using GAMLSS. The joint probability of runoff series using copula is calculated to assess the dryness-wetness encounter probability. The main results of this study are as follows: • Significant non-stationarity has been detected for DJK, HYK, and UR. Using MK and Sen's slope trend tests, there was a gradual fall in the runoff series of DJK and sharp drops in HYR and UR. Using Pettitt, SNHT, and Mann-Kendall tests, 1990, 1985, and 1990, and 1996 are considered as the change points of runoff in DJK, HYR, and UR, respectively. • The non-stationary model with time as the explanatory variable is more suitable for the runoff series of DJK, HYK, and UR. Under the non-stationary framework, the wet-dry classification of runoff series is completely changed. The transfer matrix from the stationary to non-stationary model in DJK indicated shows that eight out of 58 years in the classification result changed. For HYK and UR, the amount of change years is 29 and 19, respectively. • For the encounter probability of bivariate combination types, the asynchronous probabilities of DJK-UR and HYK-UR are over 0.6, which indicates that the interbasin water transfer project is a reasonable and scientific way for the BYD Lake ecological water replenishment. • For the encounter probability of trivariate combination types, the asynchronous probability of runoff series of three water sources (P A ) is 0.82, and the asynchronous probability of DJK and HYK is 0.602, which indicates that there is no substitution between DJK and HYK. A single water source is difficult to guarantee the ecological water demand of the BYD Lake. • During a dry year of UR, the conditional probability that both DJK and HYK are in a dry year is 0.234. Thus, some additional water resources and a reasonable EWR plan are required to solve the case when no water is available.