# Water Supply Risk Analysis Based on Runoff Sequence Simulation with Change Point under Changing Environment

Authors: Yafang Cheng; Ping Feng; Jianzhu Li; Yuangang Guo; Peizhen Ren
Journal: Advances in Meteorology
 (2019)

---

## Abstract

This study investigates the water supply risk of Panjiakou reservoir in the Luanhe River basin of China during 1956–2016 under environmental change. Since the monthly runoff series during 1956–2016 is a time series with change points, it is necessary to find a new stochastic streamflow series generation approach to preserve the statistical characteristics of the original series and to refine the reliability of water supply risk analysis. This paper improves a known stochastic streamflow simulation method of previous research to better reflect the characteristics of series with change points. And this paper simulates the monthly runoff series with change point of Panjiakou reservoir during 1956–2016 by three different methods, including Thomas–Fiering model, copula function stochastic simulation method, and copula function stochastic simulation method with the mixed distribution model. Among the three methods, the copula function stochastic simulation method with the mixed distribution model which is improved on the basis of copula function stochastic simulation method in this study performs best in simulating the observed monthly runoff series during 1956–2016, and the water supply risk indices including reliability (time-based and volume-based), resilience, vulnerability, drought risk index (DRI), and sustainability index (SUI) are evaluated for Panjiakou reservoir and analyzed by using the stochastic simulation results. By comparing with the previous studies, all indicators are between the corresponding results of 1956–1979 and 1980–2016 with stationary inflows; it can be seen that change point seriously affects the water supply risk of Panjiakou reservoir. These results make it easy to formulate water supply strategies and schemes in changing environment for water resources managers.


---

## Body

## 1. Introduction

Dams and reservoirs play a pivotal role in the development of the country. It is the most expensive element in the multipurpose river basin. They require very careful planning, design, construction, and operation. A reservoir is a manmade lake or large freshwater body which is used for multiple purposes such as for hydropower generation and water supply system, and the main function is to regulate surface runoff to satisfy certain requirements. The problem of the water supply system and the uncertainty factors such as randomness of rainfall and runoff make the great uncertainties in water supply and demand, resulting in some certain risks in the water supply system [1].According to the reliable scientific studies, research in the area of water supply risk is going momentum in order to mitigate the risk in the field of hydrological science [2, 3] and how to calculate the water supply risk becomes predominantly necessary. Three classical indicators (reliability, resilience, and vulnerability) were proposed to calculate the water supply risk of reservoirs and were defined mathematically [4]. Progressively, some composite indicators have been planned in order to better assess supply risk of the water resource system. Based on the three single indicators proposed by Hashimoto et al. [4], Xu et al. [1] put forward two composite indicators, including drought risk index (DRI) and drought damage index (DDI), and conducted a risk analysis of the water supply system in Fukuoka. This was in the first time that the three indicators of reliability, resilience, and vulnerability were weighted separately and collaboratively and drought risk index (DRI) is established to describe the probability and loss degree of water supply accidents. The weight index can reflect the contradiction between normal and failed water supply. However, it is very difficult to determine the weight, so it is necessary to make a comprehensive analysis of the uncertain factors affecting system security. Besides, Loucks [5] presented the sustainability index (SUI), which is evaluated by three single indicators through a certain mathematical formula, and the value is between 0 and 1.If the reliability or resilience index value is zero, then the SUI is close to one, indicating that the water resource system is relatively stable. Basis on the fifth assessment report of the Intergovernmental Panel on Climate Change (IPCC), the frequency and intensity of extreme weather events such as droughts and floods will increase significantly [6].With the changing of global climate and increasing demand of water, drought becomes a growing serious problem and attracts worldwide. Obviously, the shortage of water resources seriously affects urban and economic development. Thus, assessing the impact of drought on water supply system performance has become a key research topic in recent years [7–9]. Feng [10] applied the water supply risk analysis method to a water resource system during drought periods and calculated specific risk indicators such as reliability, resilience, and vulnerability for a reservoir. Due to the rapid development of socioeconomic activities such as water conservancy construction, urbanization, and industrialization, the “stationarity” hypothesis of the runoff series is no longer satisfied, which will bring a difficulty to hydrological analysis, especially in the calculation of water supply risk [11–13]. The calculation of water supply risk is based on stochastic simulation, which is the basis of measured runoff series. Consequently, uncertainty is brought to the calculation of water supply risk when change point exists in the measured runoff series. How to reduce this uncertainty is a great challenge in water supply risk analysis.Mixed distribution (MD) was put forward by Singh and Sinclair [14] for the first time. Two kinds of distribution were mixed to one distribution, and it was found that the fitting result was better than single distribution through the study of flood series with a change point in several river basins. Alila and Mtiraoui [15] reached the same conclusion that the MD method fitted better than the traditional single distribution model in the Jila River basin. Besides, they pointed out that one key to ensuring the accuracy of the MD model is to keep the number of subdistribution minimum as far as possible, mainly because the increase of the number of subdistributions will increase the number of parameters and affect the accuracy of parameter estimation. In general, mixed distribution is a good method to simulate the measured runoff, and it can preserve the statistical characteristics of runoff series [16–20].The aims of the present study are the following: (i) generate the runoff series with change point by using three stochastic simulation methods including the Thomas–Fiering model (T-F model), copula function stochastic simulation method (C(1) model), and copula function stochastic simulation method with mixed distribution model (C(m) model); (ii) choose the best stochastic simulation method to assess water supply risk of Panjiakou reservoir during 1956–2016. The innovation of this study is that the C(m) model is improved by mixed distribution on the basis of C(1) model proposed in previous research, which can better simulate the runoff series stochastic simulation with change point. Moreover, the water supply risk of Panjiakou reservoir can be analyzed more accurately.


## 2. Study Area and Datasets

### 2.1. Study Area

The Luanhe river basin, located in the northeastern China, has a drainage area of approximately 44750 km2. The basin mainly includes three types of landforms: plateau, mountainous area, and plain. The terrain inclines from northwest to southeast. It is a typical temperate continental monsoon climate with an average annual precipitation of 400–700 mm. The spatial difference of precipitation distribution is significant, and the annual distribution is uneven. The precipitation concentrates in summer which is about 200–550 mm, accounting for 66–76% of the annual precipitation.Panjiakou reservoir is located at the junction of Tangshan and Chengde cities, Hebei Province, on the main Luanhe River (Figure1), which was built in 1979. The main function of the reservoir is water supply, as well as flood control and power generation. The elevation of the main dam of the Panjiakou reservoir is 230.5 m, the normal water level is 222.0 m, the check flood water level is 227.0 m, the design flood water level is 224.5 m, and the flood limit water level is 216.0 m. The storage capacity is 1.95 billion m3, the flood control storage is 970 million m3, and the total reservoir storage capacity is 2.93 billion m3. The area controlled by the Panjiakou reservoir is 33700 km2, accounting for 75.3% of the entire area of the Luanhe river basin. The long-term average annual runoff of Panjiakou reservoir is 2.45 billion m3 [21]. In recent years, the annual runoff in the Panjiakou reservoir shows a decreasing trend which has become more significant since 1979, and drought disaster has become increasingly prominent [22].Figure 1
Location of Panjiakou reservoir in the Luanhe river basin.


### 2.2. Datasets

In this paper, the monthly streamflow series observed at Panjiakou reservoir during 1956–2016 are utilized as the basis of stochastic runoff simulation. Three different water supply guarantee rates for water volumes in 75% (1.95 billion m3/year), 85% (1.5 billion m3/year), and 95% (1.1 billion m3/year) are treated as the water demand of Panjiakou reservoir, respectively. The actual distribution of the average annual water demand of 2000–2010 is referenced to disaggregate the water volumes into monthly water demand shown in Figure 2(a) [23]. Panjiakou reservoir evaporation loss refers to the conversion coefficient of water surface evaporation in the Luanhe river basin (Figure 2(b)). Ecological water requirements of downstream channel are proposed by Wang et al. [24] (Figure 2(c)).Figure 2
Water use in different aspects for Panjiakou reservoir: (a) three different water supply guarantee rates of demand during 2000–2010; (b) water surface evaporation; (c) three different levels of ecological water requirements.
(a)
(b)
(c)


## 2.1. Study Area

The Luanhe river basin, located in the northeastern China, has a drainage area of approximately 44750 km2. The basin mainly includes three types of landforms: plateau, mountainous area, and plain. The terrain inclines from northwest to southeast. It is a typical temperate continental monsoon climate with an average annual precipitation of 400–700 mm. The spatial difference of precipitation distribution is significant, and the annual distribution is uneven. The precipitation concentrates in summer which is about 200–550 mm, accounting for 66–76% of the annual precipitation.Panjiakou reservoir is located at the junction of Tangshan and Chengde cities, Hebei Province, on the main Luanhe River (Figure1), which was built in 1979. The main function of the reservoir is water supply, as well as flood control and power generation. The elevation of the main dam of the Panjiakou reservoir is 230.5 m, the normal water level is 222.0 m, the check flood water level is 227.0 m, the design flood water level is 224.5 m, and the flood limit water level is 216.0 m. The storage capacity is 1.95 billion m3, the flood control storage is 970 million m3, and the total reservoir storage capacity is 2.93 billion m3. The area controlled by the Panjiakou reservoir is 33700 km2, accounting for 75.3% of the entire area of the Luanhe river basin. The long-term average annual runoff of Panjiakou reservoir is 2.45 billion m3 [21]. In recent years, the annual runoff in the Panjiakou reservoir shows a decreasing trend which has become more significant since 1979, and drought disaster has become increasingly prominent [22].Figure 1
Location of Panjiakou reservoir in the Luanhe river basin.


## 2.2. Datasets

In this paper, the monthly streamflow series observed at Panjiakou reservoir during 1956–2016 are utilized as the basis of stochastic runoff simulation. Three different water supply guarantee rates for water volumes in 75% (1.95 billion m3/year), 85% (1.5 billion m3/year), and 95% (1.1 billion m3/year) are treated as the water demand of Panjiakou reservoir, respectively. The actual distribution of the average annual water demand of 2000–2010 is referenced to disaggregate the water volumes into monthly water demand shown in Figure 2(a) [23]. Panjiakou reservoir evaporation loss refers to the conversion coefficient of water surface evaporation in the Luanhe river basin (Figure 2(b)). Ecological water requirements of downstream channel are proposed by Wang et al. [24] (Figure 2(c)).Figure 2
Water use in different aspects for Panjiakou reservoir: (a) three different water supply guarantee rates of demand during 2000–2010; (b) water surface evaporation; (c) three different levels of ecological water requirements.
(a)
(b)
(c)


## 3. Methodology

### 3.1. Runoff Time Series Generation

#### 3.1.1. Thomas–Fiering Model

The Thomas–Fiering model was developed by Thomas and Fiering in 1962 for the sequential generation of nonhistoric monthly streamflow [25]. The original T-F model is essentially a cyclostationary version of the classic stationary linear autoregressive model, in order to account for systematic changes and nonstationarities of statistical characteristics across seasons [26]. The marginal distributions of many hydrometeorological processes are not Gaussian, which motivated Thomas and Fiering [27] propose to replace the Gaussian white noise with Gamma (G) or Pearson type-III (PIII) distributed white noise in order to account for the high skewness coefficient [28]. Besides, the common approaches of cyclostationary modification with log-normal distribution for T-F model were proposed by Matalas [29]. And Harms et al. found that the monthly flows provide a consistently better fit to the log-normal distribution and then proposed a method of an extension to the Thomas–Fiering model which can provide an authentic representation of streamflow [30]. This method can preserve the log-normal distribution of monthly flows and correlation between monthly flows. Therefore, it is applied into the stochastic simulation process of monthly runoff and participates in the comparison of simulation results. Let i be the serial number of the year and j be the serial number of the month in the year. Initial streamflow of annual sequence and successive flows are as follows:(1)Q1=Q′¯+z1S,Qi+1=Q′¯+RQi−Q′¯+z1S1−R2,where Qi is the generated annual flow in the ith year; Q′¯ is the average measured annual flow; zj is the random normal deviation with zero mean and unit variance; S is the standard deviation of measured annual flow; and R is the autocorrelation coefficient of annual flow.When the monthly runoff distribution function is log-normal distribution, the general formula of monthly runoff simulation is as follows:(2)logq1,10=logq1′¯+z1s1,logqi,j0=logqj¯+bjlogqi,j−10−logqj−1¯+zjsj1−rj,j−12,qi,j=365qi,j0Qi∑j=112ajqi,j0,where qi,j0 is the initial estimate for qi,j; qj¯ is the average measured monthly flow; qj′ is the measured monthly flow for month j; sj is the standard deviation of measured monthly flow; rj,j−1 is the correlation coefficient between month j − 1 and month j; bj is the regression coefficient estimating logqj′ from logq′j−1; zj is independent for each month and it is independent from qi,j−1; and aj is the corresponding number of days in the monthly sequences used.


#### 3.1.2. Copula Function Stochastic Simulation Method

The copula function stochastic simulation method is proposed by Borgomeo et al. [3], which is a copula-based approach to test the vulnerability of the water resource system. The steps of the copula function stochastic simulation method are as follows [3]: The joint probability distribution models of 12 pairs of measured runoff marginal probability distribution in adjacent months are established based on Clayton copula function and the parameters are obtained. The bootstrapping method is used to sample the measured monthly runoff series and generate the pseudo-time series of monthly runoff, and the corresponding probability series is obtained for the generated pseudo-time series of monthly runoff. The joint probability formula of copula function is used to generate probability series for the next month. The runoff series of the next month is calculated by combining the marginal probability distribution (log-normal distribution is adopted in this paper). By analogy, monthly runoff series can be generated.


#### 3.1.3. Copula Function Stochastic Simulation Method with Mixed Distribution Model

In this paper, the stochastic runoff time series generation method using mixed distribution is improved based on the copula function stochastic simulation method proposed by Borgomeo et al. [3], and the parameter θ value of Clayton copula function of Panjiakou reservoir calculated by Li et al. [31] is used to generate the stochastic series of monthly runoff by bootstrapping the observed monthly runoff data. Then, the monthly marginal probability distribution before and after the change point of each month is calculated, and the before and after the change point distribution will be mixed by the mixed distribution model. At last, monthly runoff series is generated by the joint probability distribution so as to calculate the water supply risk index value.(1) Clayton Copula Function. For commonly used Archimedean copula functions, different function forms are different in describing the correlation structure among variables [32]. The upper and lower tail correlation coefficients of Gumbel copula function are 2−21/θ and 0, respectively, which show that it is more suitable for describing hydrological variables with upper tail correlation, such as the joint distribution with the upper tail correlation between flood peak and period of time flood volumes and period of time flood volumes and flood duration [33]; The upper and lower tail correlation coefficients of Frank copula function are all 0; namely, it does not require the correlation among variables and is suitable for describing variables without tail correlation. The upper and lower tail correlation coefficients of Clayton copula function are 0 and 2−1/θ, respectively, so it is more suitable for describing hydrological variables with lower tail correlation, such as drought analysis. Among the three copula functions, Clayton copula function is the best for the stochastic simulation of the Panjiakou reservoir [31]. The form of two-dimensional symmetric Clayton copula function is as follows:(3)Cθu,v=CθF1x1,F2x2=u−θ+v−θ−1−1/θ,θ>0,where Cθ is the two-dimensional copula function; u and v are defined as realizations of the random variables U=F1x1 and V=F2x2; and θ is the parameter of copula function, which changes throughout the year for different pairs of consecutive months.(2) Mixed Distribution. In order to better simulate the monthly runoff series of Panjiakou reservoir with change point during 1956–2016, in this study, the mixed distribution method is used to optimize the random simulation method of measured runoff series to characterize the change point of runoff series. The mixed distribution method assumes that the series of extreme samples with different distributions are composed of several subseries distributions [30]. It can be expressed as follows:(4)Fx=α1F1x+α2F2x+…+αkFkx,where Fx is the cumulative distribution function composed of k subseries distribution and αk is the weight of each subseries distribution and satisfies with the equation of αk>0 and ∑i=1kαk=1.In order to reduce the complexity of parameter estimation, it is generally assumed that the mixed distribution is composed of two subseries distributions. For a hydrological seriesX with change point, if the sample size is n and the year of variation is τ, it is assumed that the series before the year of variation is X1, following the distribution of probability density function f1x and the sample length is n1=τ, while the series after the year of variation is X2, following the distribution with probability density function f2x and the sample length is n2=n−τ. Then, the whole series X follows the mixed distribution of the probability density function fx which can be expressed as(5)fx=αf1x+1−αf2x,where α is the weight.In this study, subseries all match basically to log-normal distribution and the probability density functions aref1x and f2x, respectively. The expressions are(6)f1x=1xσ12πexp−lnx−μ122σ12,f2x=1xσ22πexp−lnx−μ222σ22.Then, the theoretical frequency formula of the mixed distributionFx of series with change point is as follows:(7)Fx=α1−1xσ12πexp−lnx−μ122σ12dt+1−α1−1xσ22πexp−lnx−μ222σ22,where σi and μi are the standard deviation parameter and mean parameter of lnx, respectively. So, there will be 5 parameters σ1,σ2,μ1,μ2,α for each month when using a mixed distribution function.Previous studies have shown that the runoff sequence of Panjiakou reservoir is nonstationary, which means that sudden changing years called “change points” exist in the sequence [21]. In this study, the subseries X1 is the observed monthly runoff series before the year of variation and X2 is the observed monthly runoff series after the year of variation. The Mann–Kendall trend analysis method is selected to test the change points of monthly runoff series. It is found that from January to December, the change points appear around the year of 1979, which is consistent with the previous studies in this region [31]. The M-K method was used to test the change points of monthly runoff series. It was found that the change points from January to December were 1976, 1973, 1972, 1976, 1980, 1980, 1978, 1977, 1978, 1977, 1981, and 1977, respectively. Mixed distribution needs to integrate the runoff distribution before and after the change point of each month to form a mixed distribution function which can represent the change point so as to optimize the stochastic simulation method.(3) Simulated Annealing Algorithm. Simulated annealing algorithm is selected to estimate the parameters of the mixed distribution, and the objective function is to minimize the sum of squares of frequency deviation. The idea of simulated annealing was first proposed by Metropolis in 1953. In 1983, Kirkpatrick et al. [34] successfully introduced the idea of simulated annealing into the optimization field. The starting point of simulated annealing algorithm for optimization problem is based on the similarity between annealing process of solid material in physics and general optimization problem. The solution of system objective function is the lowest energy state reached by simulation [35]. At a certain temperature, the simulated annealing algorithm can find the global optimal solution of the objective function randomly in the solution space by combining the probabilistic jumping property with the gradual decrease of the temperature parameters. It has been proved theoretically that the simulated annealing algorithm can converge to the global optimal solution with probability 1 as long as the simulation process is adequate [36–38]. The algorithm can be applied universally and has no special requirements for initial conditions and objective functions.(4) Sum of Squares of Deviations. To compare the results of the three stochastic simulation models, root-mean-squared error (RMSE), a comprehensive evaluation index, is used to evaluate the accuracy of the models. The statistical variable is calculated as follows:(8)RMSE=1n∑i=1nδ0,j−δs,jδ0,j2,where δ0,j is the statistical parameters measured in the month j, δs,j is the statistical parameters simulated in the month j, and n is the number of month and n = 12. RMSE indicates the fitting degree between the simulated and measured values, which can reflect the simulation effect of the model to a certain extent. The smaller the RMSE value is, the better the simulation effect of the model is.


### 3.2. Water Supply Risk Analysis Model

#### 3.2.1. Reservoir Operation

In this study, Standard Operating Policy (SOP) is adopted in reservoir operation in water supply risk simulation. SOP means that the discharge of reservoir is simplified as a function of the storage of reservoir and the inflow. It is widely used in the simulation analysis of reservoir operation because of its simplicity. The scheduling target of SOP is to meet water demand as much as possible during the water supply period [4], and the execution process of its operation strategy is as follows: if the reservoir cannot meet the target water demand in the simulation period, the reservoir will be released to dead storage capacity in order to meet the water demand in the current period, while if the excessive water inflow causes the reservoir water storage to exceed the normal water level, the excess water should be released to keep the normal water level.


#### 3.2.2. Water Supply Simulation

The water supply simulation of Panjiakou reservoir is based on the water supply system planning, designing, and other operation constraints, and water supply risk indices are evaluated based on the process of simulation.(1)
Water balance of reservoir constraints:(9)Vt=Vt−1+It−Rt−Ot−Et−Wt，t=1,2,…,n,wheret is the time period; Vt and Vt−1 are the effective storage at the beginning and end of the t period, respectively; It is the runoff into reservoir in the t period; Rt is the ecological water demand of the downstream of reservoir; Ot is the discharge of reservoir; Et is the evaporation of reservoir; and Wt is the planned water supply of reservoir; and n is the simulation length of reservoir operation.(2)
Storage capacity constraints:(10)Vt≤C,t=1,2,…,n，Vn−V0≥0,whereC is the storage capacity of reservoir and V0 and Vt are the initial and end storage of reservoir operation.(3)
Relation between available water supply and planned water supply:(11)ht≤WtVt−1+It−Rt−Et,t=1,2,…,n，whereht is an indicator variable for the state of water supply. It is a 0–1 type variable. That the value is 0 means the normal water supply while that the value is 1 means the water supply is failed.(4)
Discharge of reservoir constraints:(12)bt≤Vt−1+It−Rt−Et−WtVc,t=1,2,…,n,Ot=Vt−1+It−Rt−Et−Wt−Vc,t=1,2,…,n,wherebt is a 0-1 type variable, indicating the discharge state of reservoir. That the value is 0 means there is no discharged water in the reservoir while that the value is 1 is on the contrary. Vc is the storage capacity; Ot is the value of discharged water: when the planned water supply exceeds Vc and the value of bt is 1, the reservoir has discarded water, otherwise the reservoir has no discarded water.(5)
Nonsimultaneous constraints of incident of failed water supply and discharge incident at any period:(13)ht+bt≤1,t=1,2,…,n,(6)
Model variables:Vt, Vt−1, Rt, Ot, Wt, and It are defined as the nonnegative value; htandbt are defined as the nonnegative integers.


#### 3.2.3. Water Supply Risk Indices

Reliability, resilience, and vulnerability are widely used in the evaluation of the water supply risk [4]. In this study, the three indices are utilized to quantify the capability of the reservoir to define performance indices and two composite indices are also used.LetDt be the target water requirement that needs to be met during the tth period and Xt be the state of a water supply system. S indicates the periods of satisfactory state while U indicates the periods of unsatisfactory state.(1)
Reliability can be defined as(14)rt=PXt∈S=1NS∑t=1NSht,rv=PXt∈S=1−∑t=1NSDt−Rt∑t=1NSDt,ht=1,ifXt⊂S,0,ifXt⊂U,wherert and rv, respectively, represent the reliability of a water supply system based on time and based on volume; Rt is the release from the reservoir storage during t; NS is the total number of time intervals; and ht is the number of intervals that the target water requirement has fully met.(2)
Resilience can be defined as(15)res=∑t=1T−1WtT−∑t=1Tht,Wt=1,ifXt−1⊂U,Xt⊂S,0,otherwise,whereres is the resilience of a water supply system; Wt is an indicator for the transition from unsatisfactory to satisfactory state; and T is the total length of the time series.(3)
Vulnerability and relative vulnerability are utilized to evaluate the significance of failure and defined as follows:(16)v=max∑t∈JiD−Xt,i=1,…,Nfs,vr=vDf,0<vr≤1,wherev and vr are the vulnerability and the dimensionless vulnerability metric, respectively. J1, … , JN represent periods of unsatisfactory states and Ji ≥ 1 month. fs is the number of individual continuous sequences of failure periods, and Df is the consequent target water requirement during unsatisfactory periods.(4)
Drought risk index (DRI) can be expressed as a linear combination of reliability, resilience, and relative vulnerability:(17)DRI=α11−rt+α21−res+α3vr,whereα1, α2, and α3 are the weights and α1+α2+α3=1 should be satisfied. The three weights need to be predetermined, but it is difficult to evaluate, and, in general, the three weights are treated as the equal value [1, 39].(5)
Sustainability of a water supply system can be characterized by many methods, such as fairness [40] and reversibility [41]. Therefore, one index named “sustainability” is the only measure which combines the reliability, resilience, and vulnerability [5] and defined as(18)SUI=rtres1−vr,whereSUI is the sustainability index and ranges from 0 to 1.4.


## 3.1. Runoff Time Series Generation

### 3.1.1. Thomas–Fiering Model

The Thomas–Fiering model was developed by Thomas and Fiering in 1962 for the sequential generation of nonhistoric monthly streamflow [25]. The original T-F model is essentially a cyclostationary version of the classic stationary linear autoregressive model, in order to account for systematic changes and nonstationarities of statistical characteristics across seasons [26]. The marginal distributions of many hydrometeorological processes are not Gaussian, which motivated Thomas and Fiering [27] propose to replace the Gaussian white noise with Gamma (G) or Pearson type-III (PIII) distributed white noise in order to account for the high skewness coefficient [28]. Besides, the common approaches of cyclostationary modification with log-normal distribution for T-F model were proposed by Matalas [29]. And Harms et al. found that the monthly flows provide a consistently better fit to the log-normal distribution and then proposed a method of an extension to the Thomas–Fiering model which can provide an authentic representation of streamflow [30]. This method can preserve the log-normal distribution of monthly flows and correlation between monthly flows. Therefore, it is applied into the stochastic simulation process of monthly runoff and participates in the comparison of simulation results. Let i be the serial number of the year and j be the serial number of the month in the year. Initial streamflow of annual sequence and successive flows are as follows:(1)Q1=Q′¯+z1S,Qi+1=Q′¯+RQi−Q′¯+z1S1−R2,where Qi is the generated annual flow in the ith year; Q′¯ is the average measured annual flow; zj is the random normal deviation with zero mean and unit variance; S is the standard deviation of measured annual flow; and R is the autocorrelation coefficient of annual flow.When the monthly runoff distribution function is log-normal distribution, the general formula of monthly runoff simulation is as follows:(2)logq1,10=logq1′¯+z1s1,logqi,j0=logqj¯+bjlogqi,j−10−logqj−1¯+zjsj1−rj,j−12,qi,j=365qi,j0Qi∑j=112ajqi,j0,where qi,j0 is the initial estimate for qi,j; qj¯ is the average measured monthly flow; qj′ is the measured monthly flow for month j; sj is the standard deviation of measured monthly flow; rj,j−1 is the correlation coefficient between month j − 1 and month j; bj is the regression coefficient estimating logqj′ from logq′j−1; zj is independent for each month and it is independent from qi,j−1; and aj is the corresponding number of days in the monthly sequences used.


### 3.1.2. Copula Function Stochastic Simulation Method

The copula function stochastic simulation method is proposed by Borgomeo et al. [3], which is a copula-based approach to test the vulnerability of the water resource system. The steps of the copula function stochastic simulation method are as follows [3]: The joint probability distribution models of 12 pairs of measured runoff marginal probability distribution in adjacent months are established based on Clayton copula function and the parameters are obtained. The bootstrapping method is used to sample the measured monthly runoff series and generate the pseudo-time series of monthly runoff, and the corresponding probability series is obtained for the generated pseudo-time series of monthly runoff. The joint probability formula of copula function is used to generate probability series for the next month. The runoff series of the next month is calculated by combining the marginal probability distribution (log-normal distribution is adopted in this paper). By analogy, monthly runoff series can be generated.


### 3.1.3. Copula Function Stochastic Simulation Method with Mixed Distribution Model

In this paper, the stochastic runoff time series generation method using mixed distribution is improved based on the copula function stochastic simulation method proposed by Borgomeo et al. [3], and the parameter θ value of Clayton copula function of Panjiakou reservoir calculated by Li et al. [31] is used to generate the stochastic series of monthly runoff by bootstrapping the observed monthly runoff data. Then, the monthly marginal probability distribution before and after the change point of each month is calculated, and the before and after the change point distribution will be mixed by the mixed distribution model. At last, monthly runoff series is generated by the joint probability distribution so as to calculate the water supply risk index value.(1) Clayton Copula Function. For commonly used Archimedean copula functions, different function forms are different in describing the correlation structure among variables [32]. The upper and lower tail correlation coefficients of Gumbel copula function are 2−21/θ and 0, respectively, which show that it is more suitable for describing hydrological variables with upper tail correlation, such as the joint distribution with the upper tail correlation between flood peak and period of time flood volumes and period of time flood volumes and flood duration [33]; The upper and lower tail correlation coefficients of Frank copula function are all 0; namely, it does not require the correlation among variables and is suitable for describing variables without tail correlation. The upper and lower tail correlation coefficients of Clayton copula function are 0 and 2−1/θ, respectively, so it is more suitable for describing hydrological variables with lower tail correlation, such as drought analysis. Among the three copula functions, Clayton copula function is the best for the stochastic simulation of the Panjiakou reservoir [31]. The form of two-dimensional symmetric Clayton copula function is as follows:(3)Cθu,v=CθF1x1,F2x2=u−θ+v−θ−1−1/θ,θ>0,where Cθ is the two-dimensional copula function; u and v are defined as realizations of the random variables U=F1x1 and V=F2x2; and θ is the parameter of copula function, which changes throughout the year for different pairs of consecutive months.(2) Mixed Distribution. In order to better simulate the monthly runoff series of Panjiakou reservoir with change point during 1956–2016, in this study, the mixed distribution method is used to optimize the random simulation method of measured runoff series to characterize the change point of runoff series. The mixed distribution method assumes that the series of extreme samples with different distributions are composed of several subseries distributions [30]. It can be expressed as follows:(4)Fx=α1F1x+α2F2x+…+αkFkx,where Fx is the cumulative distribution function composed of k subseries distribution and αk is the weight of each subseries distribution and satisfies with the equation of αk>0 and ∑i=1kαk=1.In order to reduce the complexity of parameter estimation, it is generally assumed that the mixed distribution is composed of two subseries distributions. For a hydrological seriesX with change point, if the sample size is n and the year of variation is τ, it is assumed that the series before the year of variation is X1, following the distribution of probability density function f1x and the sample length is n1=τ, while the series after the year of variation is X2, following the distribution with probability density function f2x and the sample length is n2=n−τ. Then, the whole series X follows the mixed distribution of the probability density function fx which can be expressed as(5)fx=αf1x+1−αf2x,where α is the weight.In this study, subseries all match basically to log-normal distribution and the probability density functions aref1x and f2x, respectively. The expressions are(6)f1x=1xσ12πexp−lnx−μ122σ12,f2x=1xσ22πexp−lnx−μ222σ22.Then, the theoretical frequency formula of the mixed distributionFx of series with change point is as follows:(7)Fx=α1−1xσ12πexp−lnx−μ122σ12dt+1−α1−1xσ22πexp−lnx−μ222σ22,where σi and μi are the standard deviation parameter and mean parameter of lnx, respectively. So, there will be 5 parameters σ1,σ2,μ1,μ2,α for each month when using a mixed distribution function.Previous studies have shown that the runoff sequence of Panjiakou reservoir is nonstationary, which means that sudden changing years called “change points” exist in the sequence [21]. In this study, the subseries X1 is the observed monthly runoff series before the year of variation and X2 is the observed monthly runoff series after the year of variation. The Mann–Kendall trend analysis method is selected to test the change points of monthly runoff series. It is found that from January to December, the change points appear around the year of 1979, which is consistent with the previous studies in this region [31]. The M-K method was used to test the change points of monthly runoff series. It was found that the change points from January to December were 1976, 1973, 1972, 1976, 1980, 1980, 1978, 1977, 1978, 1977, 1981, and 1977, respectively. Mixed distribution needs to integrate the runoff distribution before and after the change point of each month to form a mixed distribution function which can represent the change point so as to optimize the stochastic simulation method.(3) Simulated Annealing Algorithm. Simulated annealing algorithm is selected to estimate the parameters of the mixed distribution, and the objective function is to minimize the sum of squares of frequency deviation. The idea of simulated annealing was first proposed by Metropolis in 1953. In 1983, Kirkpatrick et al. [34] successfully introduced the idea of simulated annealing into the optimization field. The starting point of simulated annealing algorithm for optimization problem is based on the similarity between annealing process of solid material in physics and general optimization problem. The solution of system objective function is the lowest energy state reached by simulation [35]. At a certain temperature, the simulated annealing algorithm can find the global optimal solution of the objective function randomly in the solution space by combining the probabilistic jumping property with the gradual decrease of the temperature parameters. It has been proved theoretically that the simulated annealing algorithm can converge to the global optimal solution with probability 1 as long as the simulation process is adequate [36–38]. The algorithm can be applied universally and has no special requirements for initial conditions and objective functions.(4) Sum of Squares of Deviations. To compare the results of the three stochastic simulation models, root-mean-squared error (RMSE), a comprehensive evaluation index, is used to evaluate the accuracy of the models. The statistical variable is calculated as follows:(8)RMSE=1n∑i=1nδ0,j−δs,jδ0,j2,where δ0,j is the statistical parameters measured in the month j, δs,j is the statistical parameters simulated in the month j, and n is the number of month and n = 12. RMSE indicates the fitting degree between the simulated and measured values, which can reflect the simulation effect of the model to a certain extent. The smaller the RMSE value is, the better the simulation effect of the model is.


## 3.1.1. Thomas–Fiering Model

The Thomas–Fiering model was developed by Thomas and Fiering in 1962 for the sequential generation of nonhistoric monthly streamflow [25]. The original T-F model is essentially a cyclostationary version of the classic stationary linear autoregressive model, in order to account for systematic changes and nonstationarities of statistical characteristics across seasons [26]. The marginal distributions of many hydrometeorological processes are not Gaussian, which motivated Thomas and Fiering [27] propose to replace the Gaussian white noise with Gamma (G) or Pearson type-III (PIII) distributed white noise in order to account for the high skewness coefficient [28]. Besides, the common approaches of cyclostationary modification with log-normal distribution for T-F model were proposed by Matalas [29]. And Harms et al. found that the monthly flows provide a consistently better fit to the log-normal distribution and then proposed a method of an extension to the Thomas–Fiering model which can provide an authentic representation of streamflow [30]. This method can preserve the log-normal distribution of monthly flows and correlation between monthly flows. Therefore, it is applied into the stochastic simulation process of monthly runoff and participates in the comparison of simulation results. Let i be the serial number of the year and j be the serial number of the month in the year. Initial streamflow of annual sequence and successive flows are as follows:(1)Q1=Q′¯+z1S,Qi+1=Q′¯+RQi−Q′¯+z1S1−R2,where Qi is the generated annual flow in the ith year; Q′¯ is the average measured annual flow; zj is the random normal deviation with zero mean and unit variance; S is the standard deviation of measured annual flow; and R is the autocorrelation coefficient of annual flow.When the monthly runoff distribution function is log-normal distribution, the general formula of monthly runoff simulation is as follows:(2)logq1,10=logq1′¯+z1s1,logqi,j0=logqj¯+bjlogqi,j−10−logqj−1¯+zjsj1−rj,j−12,qi,j=365qi,j0Qi∑j=112ajqi,j0,where qi,j0 is the initial estimate for qi,j; qj¯ is the average measured monthly flow; qj′ is the measured monthly flow for month j; sj is the standard deviation of measured monthly flow; rj,j−1 is the correlation coefficient between month j − 1 and month j; bj is the regression coefficient estimating logqj′ from logq′j−1; zj is independent for each month and it is independent from qi,j−1; and aj is the corresponding number of days in the monthly sequences used.


## 3.1.2. Copula Function Stochastic Simulation Method

The copula function stochastic simulation method is proposed by Borgomeo et al. [3], which is a copula-based approach to test the vulnerability of the water resource system. The steps of the copula function stochastic simulation method are as follows [3]: The joint probability distribution models of 12 pairs of measured runoff marginal probability distribution in adjacent months are established based on Clayton copula function and the parameters are obtained. The bootstrapping method is used to sample the measured monthly runoff series and generate the pseudo-time series of monthly runoff, and the corresponding probability series is obtained for the generated pseudo-time series of monthly runoff. The joint probability formula of copula function is used to generate probability series for the next month. The runoff series of the next month is calculated by combining the marginal probability distribution (log-normal distribution is adopted in this paper). By analogy, monthly runoff series can be generated.


## 3.1.3. Copula Function Stochastic Simulation Method with Mixed Distribution Model

In this paper, the stochastic runoff time series generation method using mixed distribution is improved based on the copula function stochastic simulation method proposed by Borgomeo et al. [3], and the parameter θ value of Clayton copula function of Panjiakou reservoir calculated by Li et al. [31] is used to generate the stochastic series of monthly runoff by bootstrapping the observed monthly runoff data. Then, the monthly marginal probability distribution before and after the change point of each month is calculated, and the before and after the change point distribution will be mixed by the mixed distribution model. At last, monthly runoff series is generated by the joint probability distribution so as to calculate the water supply risk index value.(1) Clayton Copula Function. For commonly used Archimedean copula functions, different function forms are different in describing the correlation structure among variables [32]. The upper and lower tail correlation coefficients of Gumbel copula function are 2−21/θ and 0, respectively, which show that it is more suitable for describing hydrological variables with upper tail correlation, such as the joint distribution with the upper tail correlation between flood peak and period of time flood volumes and period of time flood volumes and flood duration [33]; The upper and lower tail correlation coefficients of Frank copula function are all 0; namely, it does not require the correlation among variables and is suitable for describing variables without tail correlation. The upper and lower tail correlation coefficients of Clayton copula function are 0 and 2−1/θ, respectively, so it is more suitable for describing hydrological variables with lower tail correlation, such as drought analysis. Among the three copula functions, Clayton copula function is the best for the stochastic simulation of the Panjiakou reservoir [31]. The form of two-dimensional symmetric Clayton copula function is as follows:(3)Cθu,v=CθF1x1,F2x2=u−θ+v−θ−1−1/θ,θ>0,where Cθ is the two-dimensional copula function; u and v are defined as realizations of the random variables U=F1x1 and V=F2x2; and θ is the parameter of copula function, which changes throughout the year for different pairs of consecutive months.(2) Mixed Distribution. In order to better simulate the monthly runoff series of Panjiakou reservoir with change point during 1956–2016, in this study, the mixed distribution method is used to optimize the random simulation method of measured runoff series to characterize the change point of runoff series. The mixed distribution method assumes that the series of extreme samples with different distributions are composed of several subseries distributions [30]. It can be expressed as follows:(4)Fx=α1F1x+α2F2x+…+αkFkx,where Fx is the cumulative distribution function composed of k subseries distribution and αk is the weight of each subseries distribution and satisfies with the equation of αk>0 and ∑i=1kαk=1.In order to reduce the complexity of parameter estimation, it is generally assumed that the mixed distribution is composed of two subseries distributions. For a hydrological seriesX with change point, if the sample size is n and the year of variation is τ, it is assumed that the series before the year of variation is X1, following the distribution of probability density function f1x and the sample length is n1=τ, while the series after the year of variation is X2, following the distribution with probability density function f2x and the sample length is n2=n−τ. Then, the whole series X follows the mixed distribution of the probability density function fx which can be expressed as(5)fx=αf1x+1−αf2x,where α is the weight.In this study, subseries all match basically to log-normal distribution and the probability density functions aref1x and f2x, respectively. The expressions are(6)f1x=1xσ12πexp−lnx−μ122σ12,f2x=1xσ22πexp−lnx−μ222σ22.Then, the theoretical frequency formula of the mixed distributionFx of series with change point is as follows:(7)Fx=α1−1xσ12πexp−lnx−μ122σ12dt+1−α1−1xσ22πexp−lnx−μ222σ22,where σi and μi are the standard deviation parameter and mean parameter of lnx, respectively. So, there will be 5 parameters σ1,σ2,μ1,μ2,α for each month when using a mixed distribution function.Previous studies have shown that the runoff sequence of Panjiakou reservoir is nonstationary, which means that sudden changing years called “change points” exist in the sequence [21]. In this study, the subseries X1 is the observed monthly runoff series before the year of variation and X2 is the observed monthly runoff series after the year of variation. The Mann–Kendall trend analysis method is selected to test the change points of monthly runoff series. It is found that from January to December, the change points appear around the year of 1979, which is consistent with the previous studies in this region [31]. The M-K method was used to test the change points of monthly runoff series. It was found that the change points from January to December were 1976, 1973, 1972, 1976, 1980, 1980, 1978, 1977, 1978, 1977, 1981, and 1977, respectively. Mixed distribution needs to integrate the runoff distribution before and after the change point of each month to form a mixed distribution function which can represent the change point so as to optimize the stochastic simulation method.(3) Simulated Annealing Algorithm. Simulated annealing algorithm is selected to estimate the parameters of the mixed distribution, and the objective function is to minimize the sum of squares of frequency deviation. The idea of simulated annealing was first proposed by Metropolis in 1953. In 1983, Kirkpatrick et al. [34] successfully introduced the idea of simulated annealing into the optimization field. The starting point of simulated annealing algorithm for optimization problem is based on the similarity between annealing process of solid material in physics and general optimization problem. The solution of system objective function is the lowest energy state reached by simulation [35]. At a certain temperature, the simulated annealing algorithm can find the global optimal solution of the objective function randomly in the solution space by combining the probabilistic jumping property with the gradual decrease of the temperature parameters. It has been proved theoretically that the simulated annealing algorithm can converge to the global optimal solution with probability 1 as long as the simulation process is adequate [36–38]. The algorithm can be applied universally and has no special requirements for initial conditions and objective functions.(4) Sum of Squares of Deviations. To compare the results of the three stochastic simulation models, root-mean-squared error (RMSE), a comprehensive evaluation index, is used to evaluate the accuracy of the models. The statistical variable is calculated as follows:(8)RMSE=1n∑i=1nδ0,j−δs,jδ0,j2,where δ0,j is the statistical parameters measured in the month j, δs,j is the statistical parameters simulated in the month j, and n is the number of month and n = 12. RMSE indicates the fitting degree between the simulated and measured values, which can reflect the simulation effect of the model to a certain extent. The smaller the RMSE value is, the better the simulation effect of the model is.


## 3.2. Water Supply Risk Analysis Model

### 3.2.1. Reservoir Operation

In this study, Standard Operating Policy (SOP) is adopted in reservoir operation in water supply risk simulation. SOP means that the discharge of reservoir is simplified as a function of the storage of reservoir and the inflow. It is widely used in the simulation analysis of reservoir operation because of its simplicity. The scheduling target of SOP is to meet water demand as much as possible during the water supply period [4], and the execution process of its operation strategy is as follows: if the reservoir cannot meet the target water demand in the simulation period, the reservoir will be released to dead storage capacity in order to meet the water demand in the current period, while if the excessive water inflow causes the reservoir water storage to exceed the normal water level, the excess water should be released to keep the normal water level.


### 3.2.2. Water Supply Simulation

The water supply simulation of Panjiakou reservoir is based on the water supply system planning, designing, and other operation constraints, and water supply risk indices are evaluated based on the process of simulation.(1)
Water balance of reservoir constraints:(9)Vt=Vt−1+It−Rt−Ot−Et−Wt，t=1,2,…,n,wheret is the time period; Vt and Vt−1 are the effective storage at the beginning and end of the t period, respectively; It is the runoff into reservoir in the t period; Rt is the ecological water demand of the downstream of reservoir; Ot is the discharge of reservoir; Et is the evaporation of reservoir; and Wt is the planned water supply of reservoir; and n is the simulation length of reservoir operation.(2)
Storage capacity constraints:(10)Vt≤C,t=1,2,…,n，Vn−V0≥0,whereC is the storage capacity of reservoir and V0 and Vt are the initial and end storage of reservoir operation.(3)
Relation between available water supply and planned water supply:(11)ht≤WtVt−1+It−Rt−Et,t=1,2,…,n，whereht is an indicator variable for the state of water supply. It is a 0–1 type variable. That the value is 0 means the normal water supply while that the value is 1 means the water supply is failed.(4)
Discharge of reservoir constraints:(12)bt≤Vt−1+It−Rt−Et−WtVc,t=1,2,…,n,Ot=Vt−1+It−Rt−Et−Wt−Vc,t=1,2,…,n,wherebt is a 0-1 type variable, indicating the discharge state of reservoir. That the value is 0 means there is no discharged water in the reservoir while that the value is 1 is on the contrary. Vc is the storage capacity; Ot is the value of discharged water: when the planned water supply exceeds Vc and the value of bt is 1, the reservoir has discarded water, otherwise the reservoir has no discarded water.(5)
Nonsimultaneous constraints of incident of failed water supply and discharge incident at any period:(13)ht+bt≤1,t=1,2,…,n,(6)
Model variables:Vt, Vt−1, Rt, Ot, Wt, and It are defined as the nonnegative value; htandbt are defined as the nonnegative integers.


### 3.2.3. Water Supply Risk Indices

Reliability, resilience, and vulnerability are widely used in the evaluation of the water supply risk [4]. In this study, the three indices are utilized to quantify the capability of the reservoir to define performance indices and two composite indices are also used.LetDt be the target water requirement that needs to be met during the tth period and Xt be the state of a water supply system. S indicates the periods of satisfactory state while U indicates the periods of unsatisfactory state.(1)
Reliability can be defined as(14)rt=PXt∈S=1NS∑t=1NSht,rv=PXt∈S=1−∑t=1NSDt−Rt∑t=1NSDt,ht=1,ifXt⊂S,0,ifXt⊂U,wherert and rv, respectively, represent the reliability of a water supply system based on time and based on volume; Rt is the release from the reservoir storage during t; NS is the total number of time intervals; and ht is the number of intervals that the target water requirement has fully met.(2)
Resilience can be defined as(15)res=∑t=1T−1WtT−∑t=1Tht,Wt=1,ifXt−1⊂U,Xt⊂S,0,otherwise,whereres is the resilience of a water supply system; Wt is an indicator for the transition from unsatisfactory to satisfactory state; and T is the total length of the time series.(3)
Vulnerability and relative vulnerability are utilized to evaluate the significance of failure and defined as follows:(16)v=max∑t∈JiD−Xt,i=1,…,Nfs,vr=vDf,0<vr≤1,wherev and vr are the vulnerability and the dimensionless vulnerability metric, respectively. J1, … , JN represent periods of unsatisfactory states and Ji ≥ 1 month. fs is the number of individual continuous sequences of failure periods, and Df is the consequent target water requirement during unsatisfactory periods.(4)
Drought risk index (DRI) can be expressed as a linear combination of reliability, resilience, and relative vulnerability:(17)DRI=α11−rt+α21−res+α3vr,whereα1, α2, and α3 are the weights and α1+α2+α3=1 should be satisfied. The three weights need to be predetermined, but it is difficult to evaluate, and, in general, the three weights are treated as the equal value [1, 39].(5)
Sustainability of a water supply system can be characterized by many methods, such as fairness [40] and reversibility [41]. Therefore, one index named “sustainability” is the only measure which combines the reliability, resilience, and vulnerability [5] and defined as(18)SUI=rtres1−vr,whereSUI is the sustainability index and ranges from 0 to 1.4.


## 3.2.1. Reservoir Operation

In this study, Standard Operating Policy (SOP) is adopted in reservoir operation in water supply risk simulation. SOP means that the discharge of reservoir is simplified as a function of the storage of reservoir and the inflow. It is widely used in the simulation analysis of reservoir operation because of its simplicity. The scheduling target of SOP is to meet water demand as much as possible during the water supply period [4], and the execution process of its operation strategy is as follows: if the reservoir cannot meet the target water demand in the simulation period, the reservoir will be released to dead storage capacity in order to meet the water demand in the current period, while if the excessive water inflow causes the reservoir water storage to exceed the normal water level, the excess water should be released to keep the normal water level.


## 3.2.2. Water Supply Simulation

The water supply simulation of Panjiakou reservoir is based on the water supply system planning, designing, and other operation constraints, and water supply risk indices are evaluated based on the process of simulation.(1)
Water balance of reservoir constraints:(9)Vt=Vt−1+It−Rt−Ot−Et−Wt，t=1,2,…,n,wheret is the time period; Vt and Vt−1 are the effective storage at the beginning and end of the t period, respectively; It is the runoff into reservoir in the t period; Rt is the ecological water demand of the downstream of reservoir; Ot is the discharge of reservoir; Et is the evaporation of reservoir; and Wt is the planned water supply of reservoir; and n is the simulation length of reservoir operation.(2)
Storage capacity constraints:(10)Vt≤C,t=1,2,…,n，Vn−V0≥0,whereC is the storage capacity of reservoir and V0 and Vt are the initial and end storage of reservoir operation.(3)
Relation between available water supply and planned water supply:(11)ht≤WtVt−1+It−Rt−Et,t=1,2,…,n，whereht is an indicator variable for the state of water supply. It is a 0–1 type variable. That the value is 0 means the normal water supply while that the value is 1 means the water supply is failed.(4)
Discharge of reservoir constraints:(12)bt≤Vt−1+It−Rt−Et−WtVc,t=1,2,…,n,Ot=Vt−1+It−Rt−Et−Wt−Vc,t=1,2,…,n,wherebt is a 0-1 type variable, indicating the discharge state of reservoir. That the value is 0 means there is no discharged water in the reservoir while that the value is 1 is on the contrary. Vc is the storage capacity; Ot is the value of discharged water: when the planned water supply exceeds Vc and the value of bt is 1, the reservoir has discarded water, otherwise the reservoir has no discarded water.(5)
Nonsimultaneous constraints of incident of failed water supply and discharge incident at any period:(13)ht+bt≤1,t=1,2,…,n,(6)
Model variables:Vt, Vt−1, Rt, Ot, Wt, and It are defined as the nonnegative value; htandbt are defined as the nonnegative integers.


## 3.2.3. Water Supply Risk Indices

Reliability, resilience, and vulnerability are widely used in the evaluation of the water supply risk [4]. In this study, the three indices are utilized to quantify the capability of the reservoir to define performance indices and two composite indices are also used.LetDt be the target water requirement that needs to be met during the tth period and Xt be the state of a water supply system. S indicates the periods of satisfactory state while U indicates the periods of unsatisfactory state.(1)
Reliability can be defined as(14)rt=PXt∈S=1NS∑t=1NSht,rv=PXt∈S=1−∑t=1NSDt−Rt∑t=1NSDt,ht=1,ifXt⊂S,0,ifXt⊂U,wherert and rv, respectively, represent the reliability of a water supply system based on time and based on volume; Rt is the release from the reservoir storage during t; NS is the total number of time intervals; and ht is the number of intervals that the target water requirement has fully met.(2)
Resilience can be defined as(15)res=∑t=1T−1WtT−∑t=1Tht,Wt=1,ifXt−1⊂U,Xt⊂S,0,otherwise,whereres is the resilience of a water supply system; Wt is an indicator for the transition from unsatisfactory to satisfactory state; and T is the total length of the time series.(3)
Vulnerability and relative vulnerability are utilized to evaluate the significance of failure and defined as follows:(16)v=max∑t∈JiD−Xt,i=1,…,Nfs,vr=vDf,0<vr≤1,wherev and vr are the vulnerability and the dimensionless vulnerability metric, respectively. J1, … , JN represent periods of unsatisfactory states and Ji ≥ 1 month. fs is the number of individual continuous sequences of failure periods, and Df is the consequent target water requirement during unsatisfactory periods.(4)
Drought risk index (DRI) can be expressed as a linear combination of reliability, resilience, and relative vulnerability:(17)DRI=α11−rt+α21−res+α3vr,whereα1, α2, and α3 are the weights and α1+α2+α3=1 should be satisfied. The three weights need to be predetermined, but it is difficult to evaluate, and, in general, the three weights are treated as the equal value [1, 39].(5)
Sustainability of a water supply system can be characterized by many methods, such as fairness [40] and reversibility [41]. Therefore, one index named “sustainability” is the only measure which combines the reliability, resilience, and vulnerability [5] and defined as(18)SUI=rtres1−vr,whereSUI is the sustainability index and ranges from 0 to 1.4.


## 4. Results

### 4.1. Runoff Time Series Generation

The monthly runoff series with change point of Panjiakou reservoir in 1956–2016 are simulated by the Thomas–Fiering model (T-F model), copula function stochastic simulation method (C(1) model), and copula function stochastic simulation method with mixed distribution model (C(m) model). All the three models generate monthly runoff series of 1000 years. The frequency distribution curves for measured data during 1956–1979 (before change point) and 1980–2016 (after change point) by log-normal distribution are shown in Figure 3. In Figure 3, it can be seen that the probability distribution before and after the change point is quite different, so mixed distribution is considered to utilize instead of single distribution. And the mixed-two log-normal frequency distribution curves of Panjiakou reservoir monthly runoff are shown in Figure 4.Figure 3
Frequency distribution curves of Panjiakou reservoir monthly runoff for different series (log-normal distribution). (a) Jan. (b) Feb. (c) Mar. (d) Apr. (e) May. (f) June. (g) July. (h) Aug. (i) Sept. (j) Oct. (k) Nov. (l) Dec.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)Figure 4
Mixed-two log-normal frequency distribution curves of Panjiakou reservoir monthly runoff. (a) Jan. (b) Feb. (c) Mar. (d) Apr. (e) May. (f) June. (g) July. (h) Aug. (i) Sept. (j) Oct. (k) Nov. (l) Dec.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)Then, the results of stochastic simulation are compared, which mean the statistical results of mean flow, standard deviation, variation coefficientCV, monthly flow skewness CS, and correlation coefficient r of adjacent months are compared. The comparison of statistical results of the three stochastic simulations is shown in Table 1 and Figure 5. By comparing the monthly measured statistical parameters of the T-F model, C(1) model, and C(m) model, it can be found that the mean flow simulated by the C(1) model and C(m) model is similar to those of the T-F model. The RMSE of the T-F model is 0.004, the RMSE of the C(1) model is 0.031, and the RMSE of the C(m) model is 0.005. It shows that three models are not significantly different in this index, as standard deviation, T-F model, and C(m) model have obvious advantages. The RMSE for CV of the T-F model is 0.068 and C(1) model is 0.111, but C(m) model is only 0.060. In correlation coefficient r, C(m) model also has some optimization on the basis of the C(1) model. Although the simulation results of T-F model (0.180) are slightly better than those of C(1) model (0.211) and C(m) model (0.181) in terms of correlation coefficient r of adjacent months, the advantage is not obvious. This is because T-F model is based on the linear correlation coefficient, while C(1) model and C(m) model are simulated on the basis of nonlinear correlation coefficient. In summary, by comparing the simulation results of the T-F model, C(1) model, and C(m) model, the C(m) model performs best in preserving and reflecting the original sequence of statistical parameters. This also proves that the C(m) model which uses mixed distribution can optimize the C(1) model to some extent. It can better simulate the statistical characteristics of the monthly runoff series with change point of Panjiakou reservoir by the C(m) model, so the C(m) model is used in the stochastic simulation to calculate water supply risk of Panjiakou reservoir.Table 1
Comparison of stochastic simulation results of monthly runoff for Panjiakou reservoir during 1956–2016.
Month
Mean flow
Standard deviation
C
V
r
Observed
T-F
C(1)
C(m)
Observed
T-F
C(1)
C(m)
Observed
T-F
C(1)
C(m)
Observed
T-F
C(1)
C(m)
1
10.803
10.750
10.509
10.776
2.105
1.903
2.068
2.025
0.195
0.177
0.197
0.190
0.948
0.811
0.677
0.778
2
11.118
11.105
11.033
11.15
2.22
2.050
2.175
2.146
0.2
0.185
0.197
0.191
0.873
0.710
0.721
0.731
3
19.91
19.981
19.831
20.006
3.024
3.434
3.094
3.146
0.152
0.172
0.156
0.157
0.662
0.747
0.574
0.53
4
30.44
30.565
30.251
30.297
3.812
3.634
3.856
3.748
0.125
0.119
0.127
0.120
0.771
0.615
0.568
0.552
5
18.926
19.099
18.474
18.937
3.212
3.520
4.249
3.103
0.17
0.184
0.230
0.167
0.351
0.372
0.45
0.369
6
34.636
34.420
32.772
34.143
5.913
6.040
4.942
4.914
0.171
0.175
0.151
0.144
0.34
0.403
0.405
0.414
7
155.712
155.597
158.216
155.435
12.588
12.782
11.797
12.29
0.081
0.082
0.075
0.079
0.489
0.682
0.66
0.623
8
204.781
204.750
207.772
204.698
14.02
14.145
14.539
14.414
0.068
0.069
0.070
0.070
0.841
0.843
0.702
0.73
9
92.999
92.707
87.153
92.654
10.313
9.532
9.825
9.955
0.111
0.103
0.113
0.106
0.922
0.854
0.78
0.762
10
59.389
59.320
57.842
59.49
6.073
5.667
5.8
5.638
0.102
0.096
0.100
0.094
0.909
0.810
0.77
0.766
11
35.206
35.179
34.716
35.218
4.057
3.826
4.082
4.04
0.115
0.109
0.118
0.116
0.845
0.683
0.725
0.739
12
15.248
15.179
14.54
15.183
2.639
2.699
2.544
2.447
0.173
0.178
0.175
0.162
0.637
0.532
0.609
0.602
RMSE
—
0.004
0.031
0.005
—
0.070
0.109
0.063
—
0.068
0.111
0.060
—
0.180
0.211
0.181Figure 5
Comparison of stochastic simulation results of monthly runoff for Panjiakou reservoir.
(a)
(b)
(c)
(d)It can be seen from the above comparison results that theC(m) model is the most suitable stochastic simulation, which is optimized by mixed distribution and simulated annealing algorithm based on copula function to simulate the monthly sequences with change point. The calculated results of mixed distribution parameters are shown in Table 2. Then, the monthly streamflow series of Panjiakou reservoir during 1956–2016 are generated based on the C(m) model.Table 2
Parameters of mixed distribution in the monthly observed streamflow series during 1956–2016.
Month
Before the change point
After the change point
α
Mean flow (m3/s)
Standard deviation (m3/s)
Mean flow (m3/s)
Standard deviation (m3/s)
1
2.401
0.246
2.248
0.450
0.555
2
2.427
0.277
2.266
0.486
0.591
3
3.080
0.325
2.805
0.516
0.747
4
3.569
0.354
3.175
0.466
0.949
5
3.034
0.444
2.512
1.053
0.438
6
3.559
0.555
3.072
0.753
0.601
7
5.120
0.698
4.305
0.987
0.868
8
5.444
0.706
4.567
1.036
0.869
9
4.727
0.710
3.799
0.774
0.983
10
4.278
0.404
3.727
0.545
0.969
11
3.692
0.396
3.276
0.457
0.914
12
2.757
0.321
2.559
0.470
0.801


### 4.2. Water Supply Risk of Panjiakou Reservoir

Based on the simulation results of theC(m) model and the water demand data of Panjiakou reservoir, the water supply risk during 1956–2016 is calculated by using the water supply risk analysis model. The results are shown in Table 3. With the simulation of C(m) model, the statistical characteristics of the original sequence are well preserved and the water supply risk during 1956–2016 is also evaluated.Table 3
Performance of generated monthly series for Panjiakou reservoir.
Water supply guarantee rate (%)
Ecological water level requirements
r
t
r
v
r
e
s
v
v
r
DRI
SUI
75
Min
0.490
0.704
0.075
2.053
0.699
0.711
0.011
Suitable
0.372
0.621
0.054
3.118
0.743
0.772
0.005
Ideal
0.216
0.493
0.040
5.315
0.778
0.841
0.002
85
Min
0.676
0.821
0.098
1.442
0.635
0.620
0.024
Suitable
0.504
0.718
0.075
2.240
0.678
0.699
0.012
Ideal
0.290
0.558
0.045
4.627
0.753
0.806
0.003
95
Min
0.858
0.929
0.125
0.924
0.571
0.529
0.046
Suitable
0.673
0.821
0.094
1.719
0.620
0.618
0.024
Ideal
0.379
0.629
0.058
3.746
0.704
0.756
0.006It can be found that the water supply risk indicators of Panjiakou reservoir during 1956–2016 are affected by series with change point. Under the water supply guarantee rate of 75% and the minimum ecological water level, the time-based reliabilityrt and the volume-based reliability rv are 0.490 and 0.704, respectively, which are far less than the water supply guarantee rate 75%, indicating that reliability of water supply is lower than expected. The resilience res is 0.075, which shows that the hydrological drought in the water supply system lasts for a long time, and it takes a relatively long time to return to normal water supply after the water supply system is destroyed. Besides, the vulnerability v is 2.053 and the relative vulnerability vr is 0.699, indicating that the water shortage in Panjiakou reservoir is serious and easy to occur. In addition, the two corresponding composite indices are as follows: the drought risk index DRI is 0.711, which indicates that the risk of drought is relatively high; the sustainability index SUI is only 0.011, which means that the ability of the water supply system to operate sustainably is very weak, and the normal water supply is basically unsustainable.


### 4.3. Comparison with Previous Research under Stationary Condition

In order to better illustrate the impact of change point on water supply risk, the results of water supply risk indices based on series with change point are compared with those of previous studies.Taking the water supply guarantee rate of 75% and the minimum ecological water level as an example, it can be found that the water supply risk indicators of Panjiakou reservoir during 1956–2016 are between the corresponding indicators of 1956–1979 and 1980–2016 [31]. As shown in Table 4, the time-based reliability rt is 0.490 (0.773 and 0.322), the volume-based reliability rv is 0.704 (0.884 and 0.554), the resilience res is 0.075 (0.137 and 0.073), the vulnerability v is 2.053 (1.700 and 2.382), the relative vulnerability vr is 0.699 (0.604 and 0.735), the drought risk index DRI is 0.711 (0.564 and 0.780), and the sustainability index SUI is 0.011 (0.042 and 0.006). From the above results, it can be seen that the water supply risk indices of Panjiakou reservoir during 1956–2016 are affected by series with change point. Compared with the water supply risk indices during 1956–1979, the reliability, resilience, and sustainability indices are reduced, while the vulnerability, relative vulnerability, and drought risk indices are increased. Overall, the water supply reliability of Panjiakou reservoir during 1956–2016 is still low, and the recovery ability after water supply destruction is still weak and vulnerable. In addition, the drought risk index DRI is still large, indicating that the risk of drought in Panjiakou reservoir is still high and the water supply system cannot operate sustainably for a long time. And under other water supply guarantee rates and other levels of ecological water demand conditions, similar conclusions are presented.Table 4
Water supply risk indices for Panjiakou reservoir based on stationary stochastic simulation.
Subseries
Water supply guarantee rate and ecological water level requirements
r
t
r
v
r
e
s
v
v
r
DRI
SUI
1956–1979
75% and min
0.773
0.884
0.137
1.700
0.604
0.564
0.042
1980–2016
75% and min
0.322
0.554
0.073
2.382
0.735
0.780
0.006


## 4.1. Runoff Time Series Generation

The monthly runoff series with change point of Panjiakou reservoir in 1956–2016 are simulated by the Thomas–Fiering model (T-F model), copula function stochastic simulation method (C(1) model), and copula function stochastic simulation method with mixed distribution model (C(m) model). All the three models generate monthly runoff series of 1000 years. The frequency distribution curves for measured data during 1956–1979 (before change point) and 1980–2016 (after change point) by log-normal distribution are shown in Figure 3. In Figure 3, it can be seen that the probability distribution before and after the change point is quite different, so mixed distribution is considered to utilize instead of single distribution. And the mixed-two log-normal frequency distribution curves of Panjiakou reservoir monthly runoff are shown in Figure 4.Figure 3
Frequency distribution curves of Panjiakou reservoir monthly runoff for different series (log-normal distribution). (a) Jan. (b) Feb. (c) Mar. (d) Apr. (e) May. (f) June. (g) July. (h) Aug. (i) Sept. (j) Oct. (k) Nov. (l) Dec.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)Figure 4
Mixed-two log-normal frequency distribution curves of Panjiakou reservoir monthly runoff. (a) Jan. (b) Feb. (c) Mar. (d) Apr. (e) May. (f) June. (g) July. (h) Aug. (i) Sept. (j) Oct. (k) Nov. (l) Dec.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)Then, the results of stochastic simulation are compared, which mean the statistical results of mean flow, standard deviation, variation coefficientCV, monthly flow skewness CS, and correlation coefficient r of adjacent months are compared. The comparison of statistical results of the three stochastic simulations is shown in Table 1 and Figure 5. By comparing the monthly measured statistical parameters of the T-F model, C(1) model, and C(m) model, it can be found that the mean flow simulated by the C(1) model and C(m) model is similar to those of the T-F model. The RMSE of the T-F model is 0.004, the RMSE of the C(1) model is 0.031, and the RMSE of the C(m) model is 0.005. It shows that three models are not significantly different in this index, as standard deviation, T-F model, and C(m) model have obvious advantages. The RMSE for CV of the T-F model is 0.068 and C(1) model is 0.111, but C(m) model is only 0.060. In correlation coefficient r, C(m) model also has some optimization on the basis of the C(1) model. Although the simulation results of T-F model (0.180) are slightly better than those of C(1) model (0.211) and C(m) model (0.181) in terms of correlation coefficient r of adjacent months, the advantage is not obvious. This is because T-F model is based on the linear correlation coefficient, while C(1) model and C(m) model are simulated on the basis of nonlinear correlation coefficient. In summary, by comparing the simulation results of the T-F model, C(1) model, and C(m) model, the C(m) model performs best in preserving and reflecting the original sequence of statistical parameters. This also proves that the C(m) model which uses mixed distribution can optimize the C(1) model to some extent. It can better simulate the statistical characteristics of the monthly runoff series with change point of Panjiakou reservoir by the C(m) model, so the C(m) model is used in the stochastic simulation to calculate water supply risk of Panjiakou reservoir.Table 1
Comparison of stochastic simulation results of monthly runoff for Panjiakou reservoir during 1956–2016.
Month
Mean flow
Standard deviation
C
V
r
Observed
T-F
C(1)
C(m)
Observed
T-F
C(1)
C(m)
Observed
T-F
C(1)
C(m)
Observed
T-F
C(1)
C(m)
1
10.803
10.750
10.509
10.776
2.105
1.903
2.068
2.025
0.195
0.177
0.197
0.190
0.948
0.811
0.677
0.778
2
11.118
11.105
11.033
11.15
2.22
2.050
2.175
2.146
0.2
0.185
0.197
0.191
0.873
0.710
0.721
0.731
3
19.91
19.981
19.831
20.006
3.024
3.434
3.094
3.146
0.152
0.172
0.156
0.157
0.662
0.747
0.574
0.53
4
30.44
30.565
30.251
30.297
3.812
3.634
3.856
3.748
0.125
0.119
0.127
0.120
0.771
0.615
0.568
0.552
5
18.926
19.099
18.474
18.937
3.212
3.520
4.249
3.103
0.17
0.184
0.230
0.167
0.351
0.372
0.45
0.369
6
34.636
34.420
32.772
34.143
5.913
6.040
4.942
4.914
0.171
0.175
0.151
0.144
0.34
0.403
0.405
0.414
7
155.712
155.597
158.216
155.435
12.588
12.782
11.797
12.29
0.081
0.082
0.075
0.079
0.489
0.682
0.66
0.623
8
204.781
204.750
207.772
204.698
14.02
14.145
14.539
14.414
0.068
0.069
0.070
0.070
0.841
0.843
0.702
0.73
9
92.999
92.707
87.153
92.654
10.313
9.532
9.825
9.955
0.111
0.103
0.113
0.106
0.922
0.854
0.78
0.762
10
59.389
59.320
57.842
59.49
6.073
5.667
5.8
5.638
0.102
0.096
0.100
0.094
0.909
0.810
0.77
0.766
11
35.206
35.179
34.716
35.218
4.057
3.826
4.082
4.04
0.115
0.109
0.118
0.116
0.845
0.683
0.725
0.739
12
15.248
15.179
14.54
15.183
2.639
2.699
2.544
2.447
0.173
0.178
0.175
0.162
0.637
0.532
0.609
0.602
RMSE
—
0.004
0.031
0.005
—
0.070
0.109
0.063
—
0.068
0.111
0.060
—
0.180
0.211
0.181Figure 5
Comparison of stochastic simulation results of monthly runoff for Panjiakou reservoir.
(a)
(b)
(c)
(d)It can be seen from the above comparison results that theC(m) model is the most suitable stochastic simulation, which is optimized by mixed distribution and simulated annealing algorithm based on copula function to simulate the monthly sequences with change point. The calculated results of mixed distribution parameters are shown in Table 2. Then, the monthly streamflow series of Panjiakou reservoir during 1956–2016 are generated based on the C(m) model.Table 2
Parameters of mixed distribution in the monthly observed streamflow series during 1956–2016.
Month
Before the change point
After the change point
α
Mean flow (m3/s)
Standard deviation (m3/s)
Mean flow (m3/s)
Standard deviation (m3/s)
1
2.401
0.246
2.248
0.450
0.555
2
2.427
0.277
2.266
0.486
0.591
3
3.080
0.325
2.805
0.516
0.747
4
3.569
0.354
3.175
0.466
0.949
5
3.034
0.444
2.512
1.053
0.438
6
3.559
0.555
3.072
0.753
0.601
7
5.120
0.698
4.305
0.987
0.868
8
5.444
0.706
4.567
1.036
0.869
9
4.727
0.710
3.799
0.774
0.983
10
4.278
0.404
3.727
0.545
0.969
11
3.692
0.396
3.276
0.457
0.914
12
2.757
0.321
2.559
0.470
0.801


## 4.2. Water Supply Risk of Panjiakou Reservoir

Based on the simulation results of theC(m) model and the water demand data of Panjiakou reservoir, the water supply risk during 1956–2016 is calculated by using the water supply risk analysis model. The results are shown in Table 3. With the simulation of C(m) model, the statistical characteristics of the original sequence are well preserved and the water supply risk during 1956–2016 is also evaluated.Table 3
Performance of generated monthly series for Panjiakou reservoir.
Water supply guarantee rate (%)
Ecological water level requirements
r
t
r
v
r
e
s
v
v
r
DRI
SUI
75
Min
0.490
0.704
0.075
2.053
0.699
0.711
0.011
Suitable
0.372
0.621
0.054
3.118
0.743
0.772
0.005
Ideal
0.216
0.493
0.040
5.315
0.778
0.841
0.002
85
Min
0.676
0.821
0.098
1.442
0.635
0.620
0.024
Suitable
0.504
0.718
0.075
2.240
0.678
0.699
0.012
Ideal
0.290
0.558
0.045
4.627
0.753
0.806
0.003
95
Min
0.858
0.929
0.125
0.924
0.571
0.529
0.046
Suitable
0.673
0.821
0.094
1.719
0.620
0.618
0.024
Ideal
0.379
0.629
0.058
3.746
0.704
0.756
0.006It can be found that the water supply risk indicators of Panjiakou reservoir during 1956–2016 are affected by series with change point. Under the water supply guarantee rate of 75% and the minimum ecological water level, the time-based reliabilityrt and the volume-based reliability rv are 0.490 and 0.704, respectively, which are far less than the water supply guarantee rate 75%, indicating that reliability of water supply is lower than expected. The resilience res is 0.075, which shows that the hydrological drought in the water supply system lasts for a long time, and it takes a relatively long time to return to normal water supply after the water supply system is destroyed. Besides, the vulnerability v is 2.053 and the relative vulnerability vr is 0.699, indicating that the water shortage in Panjiakou reservoir is serious and easy to occur. In addition, the two corresponding composite indices are as follows: the drought risk index DRI is 0.711, which indicates that the risk of drought is relatively high; the sustainability index SUI is only 0.011, which means that the ability of the water supply system to operate sustainably is very weak, and the normal water supply is basically unsustainable.


## 4.3. Comparison with Previous Research under Stationary Condition

In order to better illustrate the impact of change point on water supply risk, the results of water supply risk indices based on series with change point are compared with those of previous studies.Taking the water supply guarantee rate of 75% and the minimum ecological water level as an example, it can be found that the water supply risk indicators of Panjiakou reservoir during 1956–2016 are between the corresponding indicators of 1956–1979 and 1980–2016 [31]. As shown in Table 4, the time-based reliability rt is 0.490 (0.773 and 0.322), the volume-based reliability rv is 0.704 (0.884 and 0.554), the resilience res is 0.075 (0.137 and 0.073), the vulnerability v is 2.053 (1.700 and 2.382), the relative vulnerability vr is 0.699 (0.604 and 0.735), the drought risk index DRI is 0.711 (0.564 and 0.780), and the sustainability index SUI is 0.011 (0.042 and 0.006). From the above results, it can be seen that the water supply risk indices of Panjiakou reservoir during 1956–2016 are affected by series with change point. Compared with the water supply risk indices during 1956–1979, the reliability, resilience, and sustainability indices are reduced, while the vulnerability, relative vulnerability, and drought risk indices are increased. Overall, the water supply reliability of Panjiakou reservoir during 1956–2016 is still low, and the recovery ability after water supply destruction is still weak and vulnerable. In addition, the drought risk index DRI is still large, indicating that the risk of drought in Panjiakou reservoir is still high and the water supply system cannot operate sustainably for a long time. And under other water supply guarantee rates and other levels of ecological water demand conditions, similar conclusions are presented.Table 4
Water supply risk indices for Panjiakou reservoir based on stationary stochastic simulation.
Subseries
Water supply guarantee rate and ecological water level requirements
r
t
r
v
r
e
s
v
v
r
DRI
SUI
1956–1979
75% and min
0.773
0.884
0.137
1.700
0.604
0.564
0.042
1980–2016
75% and min
0.322
0.554
0.073
2.382
0.735
0.780
0.006


## 5. Discussion

At present, there are few water supply risk analysis based on runoff sequence simulation with change point, and there are still some problems in theoretical methods and practical applications, which need to be further explored and studied. This paper evaluates the water supply risk of Panjiakou reservoir which is a single reservoir. However, as for a region, it is often a combination of multiple water supply systems and often aimed at multiple water supply objectives [42], such as ensuring both the water quantity and quality, and even need to consider some social benefits. Under these conditions, water supply risk assessment will become very complex, which deserves further study.This study assumes that the reservoir is simulated under the standard operation strategy (SOP), but other operation strategies, such as dynamic programming strategy (DP) and stochastic dynamic programming strategy (SDP), are often used in the actual operation of the reservoir for actual demand [43–45] to supply water under different reservoir operation strategies. In different operation strategies, the change of water supply risk is a problem worth considering.At present, the most commonly used method to simulate the hydrological impact of global climate change is to use global climate models (GCMs) to simulate the changes of climate factors [46, 47], then downscale large-scale climate variables into small scale (basin scale) as input of hydrological model, and finally simulate the changes of climate factors under different scenarios. Therefore, how to analyze the risk of water supply with the global climate models (GCMs) is an interesting issue for further research to evaluate the impact of future climate conditions on water supply risk.A novel model proposed by Tsoukalas et al. shows that synthetic time series generation can simulate stationary univariate and multivariate processes with any-range dependence and arbitrary marginal distributions [48]; Gamma and Burr type-XII are also included in the target distributions. So, different dependencies and distributions can be tried for stochastic simulation to calculate water supply risk.


## 6. Conclusions

In this study, three different stochastic simulation methods were utilized to simulate the runoff series with change point. On this basis, the water supply risk indices of Panjiakou reservoir during 1956–2016 are calculated by the best stochastic simulation approach. The conclusions are as follows:(1)
The runoff series with change point of Panjiakou reservoir in 1956–2016 are simulated by the Thomas–Fiering model, copula function stochastic simulation method, and copula function stochastic simulation method with mixed distribution model, respectively. TheC(m) model is the optimized method on the basis of C(1) model of previous research. Compared with the other two methods, RMSE of the C(m) model is the lowest in standard deviation, variation coefficient CV, and monthly flow skewness CS. The copula function stochastic simulation method with mixed distribution model is the best method to preserve the statistical characteristics of runoff series with change point of Panjiakou reservoir during 1956–2016.(2)
Using the results of the Copula function stochastic simulation method with mixed distribution model and the water supply risk analysis model, the water supply risk indices (including reliability (time-based and volume-based), resilience, vulnerability, drought risk index (DRI), and sustainability index (SUI)) of Panjiakou reservoir during 1956–2016 are evaluated and analyzed. When water supply guarantee rate is 75% and the level of ecological water is minimum, the time-based reliabilityrt is 0.490, the volume-based reliability rv is 0.704, the resilience res is 0.075, the vulnerability v is 2.053, the relative vulnerability vr is 0.699, the drought risk index DRI is 0.711, and the sustainability index SUI is 0.011. The results of water supply risk indices reflect that failure of the water supply system lasts for a long time, and the time to return to normal water supply after the water supply system damaged is relatively long. Water shortage is serious, the ability of the water supply system to operate sustainably is very weak, and the normal water supply is unsustainable. Therefore, the risk of water supply of Panjiakou reservoir is relatively high, and the safety of water supply is seriously threatened.(3)
By comparing the results of water supply risk indices based on series with change point to those of previous research based on stationary stochastic simulation, the impact of change point on water supply risk is assessed. The time-based reliabilityrt, the volume-based reliability rv, the resilience res, and the sustainability index SUI of Panjiakou reservoir during 1956–2016 is lower than the corresponding indicators of 1956–1979 but higher than those of 1980–2016. However, the vulnerability v, the relative vulnerability vr and the drought risk index DRI of Panjiakou reservoir during 1956–2016 is higher than the corresponding indicators of 1956–1979 but lower than that of 1980–2016. Change point seriously affects the water supply risk, which also proves the significance of this study. This study provides a way for water supply risk analysis with change point and provides valuable information for reservoir managers to formulate water supply strategies and schemes in the changing environment.



---

*Source: 9619254-2019-09-08.xml*