Application of Empirical Mode Decomposition Method to Synthesize Flow Data: A Case Study of Hushan Reservoir in Taiwan

: Although empirical mode decomposition (EMD) was developed to analyze nonlinear and non-stationary data in the beginning, the purpose of this study is to propose a new method—based on EMD—to synthesize and generate data which be interfered with the non-stationary problems. While using EMD to decompose ﬂow record, the intrinsic mode functions and residue of a given record can be re-arranged and re-combined to generate synthetic time series with the same period. Next, the new synthetic and historical ﬂow data will be used to simulate the water supply system of Hushan reservoir, and explore the di ﬀ erence between the newly synthetic and historical ﬂow data for each goal in the water supply system of Hushan reservoir. Compared the historical ﬂow with the synthetic data generated by EMD, the synthetic data is similar to the historical ﬂow distribution overall. The ﬂow during dry season changes in signiﬁcantly ( ± 0.78 m 3 / s); however, the ﬂow distribution during wet season varies signiﬁcantly ( ± 0.63 m 3 / s). There are two analytic scenarios for demand. For Scenario I, without supporting industrial demand, the simulation results of the generation data of Method I and II show that both are more severe than the current condition, the shortage index of each method is between 0.67–1.96 but are acceptable. For Scenario II, no matter in which way the synthesis ﬂow is simulated, supporting industrial demand will seriously a ﬀ ect the equity of domestic demand, the shortage index of each method is between 1.203 and 2.12.


Introduction
Water resources are necessary for human survival, and are an important part of maintaining socio-economic development and the sustainability of the ecological environment as well. Taiwan is located at the border of tropical and subtropical areas. The main source of water resources is rainfall. The annual average rainfall in Taiwan is about 2500 mm, which is 2.5 times of the world's average. In recent years, statistics show that rainfall characteristics have become more extreme in Taiwan. Not only the raining days gradually decreasing, but the intensity gradually is increasing. On the whole, although Taiwan seems to be a country with abundant rainfall, the overall water resources environment is not easy to manageme because of climate change.
The impact of climate change tends to cause non-stationary conditions in hydrological data. Traditionally, hydrological analysis and data generation only deal with stationary data. Empirical mode decomposition (EMD) can be applied to any type of signal decomposition. Unlike singular spectrum analysis, Fourier transform, and Wavelet transform, EMD does not require any pre-determined basis functions and can extract intrinsic mode function (IMF) components from the original signal in a self-adaptive way [1]. Through the characteristics of EMD, non-stationary and nonlinear signals can be processed, which effectively overcomes the limitations of traditional methods.

Hushan Reservoir
Hushan reservoir is an off-channel reservoir. It is located between Douliu City and Gukeng Township, Yunlin County, 10 km southeast away from the Douliu City center (as shown in Figure 2). The reservoir was set up in 2016 and first reached its full water level in June 2019. The main purpose of the Hushan Reservoir is to reduce pumping groundwater and to slow down the problem of land subsidence in Yunlin county. The elevation-area-capacity relationship of Hushan Reservoir is shown in Table 1. The elevation of the remaining water level is 165 m, the full water level is 211.5 m, and the effective storage capacity of the reservoir is 53.47 × 106 m 3 . The elevations of the water outlets are 165 m and 180 m respectively. The designed flow of the domestic water channel and the permanent river outlet are 12.27 m 3 /s and 1.00 m 3 /s, respectively. The watershed of Hushan reservoir is only 6.58 km 2 and leads to limited water resources. Therefore, the Tongtou weir, the intake of Hushan reservoir, on the Qingshui river needs to be guided through the diversion tunnel.

Hushan Reservoir
Hushan reservoir is an off-channel reservoir. It is located between Douliu City and Gukeng Township, Yunlin County, 10 km southeast away from the Douliu City center (as shown in Figure 2). The reservoir was set up in 2016 and first reached its full water level in June 2019. The main purpose of the Hushan Reservoir is to reduce pumping groundwater and to slow down the problem of land subsidence in Yunlin county. The elevation-area-capacity relationship of Hushan Reservoir is shown in Table 1. The elevation of the remaining water level is 165 m, the full water level is 211.5 m, and the effective storage capacity of the reservoir is 53.47 × 106 m 3 . The elevations of the water outlets are 165 m and 180 m respectively. The designed flow of the domestic water channel and the permanent river outlet are 12.27 m 3 /s and 1.00 m 3 /s, respectively. The watershed of Hushan reservoir is only 6.58 km 2 and leads to limited water resources. Therefore, the Tongtou weir, the intake of Hushan reservoir, on the Qingshui river needs to be guided through the diversion tunnel.

Hushan Reservoir
Hushan reservoir is an off-channel reservoir. It is located between Douliu City and Gukeng Township, Yunlin County, 10 km southeast away from the Douliu City center (as shown in Figure 2). The reservoir was set up in 2016 and first reached its full water level in June 2019. The main purpose of the Hushan Reservoir is to reduce pumping groundwater and to slow down the problem of land subsidence in Yunlin county. The elevation-area-capacity relationship of Hushan Reservoir is shown in Table 1. The elevation of the remaining water level is 165 m, the full water level is 211.5 m, and the effective storage capacity of the reservoir is 53.47 × 106 m 3 . The elevations of the water outlets are 165 m and 180 m respectively. The designed flow of the domestic water channel and the permanent river outlet are 12.27 m 3 /s and 1.00 m 3 /s, respectively. The watershed of Hushan reservoir is only 6.58 km 2 and leads to limited water resources. Therefore, the Tongtou weir, the intake of Hushan reservoir, on the Qingshui river needs to be guided through the diversion tunnel.

. Tongtou Weir
Tongtou weir is the intake of Hushan reservoir. In order to ensure the river base flow and the downstream water requirements-such as agriculture and domestic-of Tongtou weir will be satisfied. The weir diversion operation will give priority to the base flow and water right; next, the intake of Hushan reservoir from Tongtou weir on Qingshui river be considered. The plan of the Tongtou weir diversion water accounts for about 10% of the Qingshui river during the wet season (May to October), and about 13% during the dry season. The Tongtou weir is located 70 m downstream of the Tongtou suspension bridge in the Qingshui river. There is a sedimentation basin before the water be taken from Tongtou weir to Hushan reservoir. The intake channel is 2.8 km long and designed intake capacity is 20 m 3 /s, with a gravity type.

Empirical Mode Decomposition
EMD is proposed for nonlinear and nonstationary time series analysis [2]. For the hydrological time series, EMD divides the sequence into two parts: IMFs and residue. The IMFs represent the implied period sequence of the data, and the residue represents the non-period sequence of the data. Although the methods used in traditional signal analysis fields such as fast Fourier transform, harmonic analysis, and wavelet analysis can also express the original sequence as a superposition of multiple period functions; these Fourier basic function-based analyses usually only deal with stationary time series. Contrary to the previous methods, EMD is intuitive, direct, a posteriori, and adaptive, with the basis of the decomposition based on, and derived from, the data [2]. EMD is more suitable for analyzing non-stationary hydrological time series.
EMD decomposes the original signal into different time scale of oscillation components called IMF. Unlike singular spectrum analysis, Fourier transform, and Wavelet transform, EMD does not require any pre-determined basis functions and can extract IMF components from the original signal in a self-adaptive way [1]. Each IMF component should satisfy the following two conditions:

1.
In the whole data series, the number of extrema and the number of zero-crossings must either equal or differ at most by 1; 2.
At any point, the mean value of the envelope defined by the local maxima and the envelope defined by the local minima is zero. An IMF represents a simple oscillatory mode as a counterpart to the simple harmonic function, but it is much more general. For time series data x(t) (t = 1, 2, . . . , n), the procedure of EMD can be described as follows [16]: Identify all the local maxima and minima of the original time series x(t); 2.
Using the three-spline interpolation function to create the upper envelopes e up (t) and the lower envelopes e low (t) of the time series; 3.
Calculate mean value m(t) of the upper and lower envelopes m(t) = e up (t) + e low (t) /2 ;

4.
Calculate the difference value d(t) between time series x(t) and mean value m(t), Check the difference value d(t): (a) if d(t) satisfies the two IMF conditions, then d(t) is defined as the ith IMF, the residue is not an IMF, then d(t) replace the x(t). 6.
Repeat (1)-(5) until the residue item r(t) becomes a monotone function or the number of extrema is less than or equal to 1, so that the IMF component cannot be decomposed again.
Finally, original time series x(t) can be denoted as sum of IMFs c i (t) and residue r(t).
where n, c i (t) and r(t) represent the number of IMF, the ith IMF and the residue, respectively. The residue r(t) also represents the overall trend or the mean value of the original time series data. Figure 3 shows the results of the EMD method decomposition. The procedure is illustrated in Huang et al. [2] and Zhang et al. [16].
Water 2020, 12, 927 6 of 20    [14]. Although EMD has obvious advantages in signal analysis, there are unavoidable defects such as the boundary-effect and mode-mixing. Mode mixing is a problem often encountered when using EMD to decompose sequences. Mode mixing means that a certain IMF contains waves with very different periodicity, often caused by intermittent signal interference. In particular, the mode-mixing will not only cause the mixing of various scale vibration modes but can even lose the physical meaning of individual IMF. To overcome this problem of the EMD method, a new noise-assisted data analysis method is proposed-the EEMD [14]-which defines the true IMF components as the mean of an ensemble of trials, each consisting of the signal plus a white noise of finite amplitude.
The main step of EEMD is described as follows [16]: 1.
Add white noise w(t) to the original time series x(t). The new time series can be defined as: 2. Decompose the new time series into IMFs using EMD method; 3.
Obtain the mean of the ensemble corresponding IMFs of the decompositions as the final result.
Differently from EMD, the EEMD method first adds white noise to the original data before decomposition and uses the statistics feature of white noise, the expectation is equal to zero, to perform the screening process. After adding white noise many times, the effect of adding white noise will be offset by the mean of the ensemble IMF. Therefore, the decomposition using EEMD not only keeps the inherent feature of the original signal, but also overcomes the mode-mixing.
The effect of the added white noise should decrease using the equation [14] ε where N is the number of ensemble members, ε is the amplitude of the added noise and ε n is the final standard deviation of error, which is defined as the difference between the input signal and the corresponding IMF(s). According to Equation (3), when the average number of times increases, the error will reduce, and when n reaches a certain number of times, its error value can be ignored.

Zero Up-Cross Method
Statistical analysis of a IMF, a wave set, requires establishment of the definition of individual wave heights and periods. The standard is the employment of the zero-crossing definition. The zero line or the line of the mean wave level is set for a record of wave registration, and each wave is defined with the reference point where the instantaneous water level crosses the zero line. When the upward crossing point is employed as the beginning of an individual wave, the method is called the zero up-crossing method [28]. The definition of upward crossing point is [29] where η i is the wave height in given time.
The wave period, T, of an individual wave will be defined as the time between two successive upward crossing point, as Figure 4 showing. For any individual wave, the peak and valley are expressed as In this study, we use the IMF decomposed by the EEMD to perform a zero up cross analysis method, analyze the periods of each IMF, and discuss it.
After finding the period in the fluctuation using the zero up cross analysis method, the sum of squares of each IMFs is regarded as its energy, and its value can correspond to the period relationship in the fluctuation. The total energy value of each group is added up to form the IMF percentage of each group compared with the energy of each group. The relational expressions of the energy percentages are where x is the square value at each time point of each group of IMF; E is the total energy; W is weight of energy.

Data Synthesis
The impact of climate change tends to cause non-stationary conditions in hydrological data. Therefore, we use the EEMD [14] in this study, and proposes a new way to synthesize and generate data that can solve the non-stationary data problems. For a given flow time series, the IMFs and residue, the EEMD decomposed result, would be replace by corresponding component in another period; next, the new combination of IMFs and residue are summed up according to the completeness of EMD [2]; then, the synthetic new flow time series will be produced. Thus, if there are many EEMD decomposed results with the same length but different period, we could synthesize a flow data set through permutation and combination of the corresponding IMFs and residues.
In this study, the historical 60-year (1956-2015) flow data of Tongtou streamflow gauge are segmented in groups of 10 years. Through the EEMD method divides the sequence into two parts: IMFs represent the inherent period sequence of the data, and the residue represents the non-period sequence of the data. We can go through the permutations and combinations of the IMFs to get n sequences. Where n is the number of groups of data grouping; i is the number of IMFs by the group.
After completing the permutations and combinations of IMFs as shown in Figure 5, we can get the new synthesis data of equal length streamflow. In this study, two permutations and combinations methods are used to synthesize new 10-year flow series data. The method is described as follows: 1. Regression analysis is performed on the different residue of each group to fit a new residue to represent the residue value. By adding the results of the IMFs permutation and combination with the representative values of the residue, we can get 10-year new synthetic flow data for the n groups.
2. Permutation and combination of the different remainders of each group with IMFs, we can get 10-year new synthetic flow data for the n groups. The average value of η max the highest point (peak) and the lowest point (valley) of η min is used as the amplitude of each wave, as shown in Equation (7) amplitude : (η max + η min )/2 (7) In this study, we use the IMF decomposed by the EEMD to perform a zero up cross analysis method, analyze the periods of each IMF, and discuss it.
After finding the period in the fluctuation using the zero up cross analysis method, the sum of squares of each IMFs is regarded as its energy, and its value can correspond to the period relationship in the fluctuation. The total energy value of each group is added up to form the IMF percentage of each group compared with the energy of each group. The relational expressions of the energy percentages are where x j 2 is the square value at each time point of each group of IMF; E i is the total energy; W i is weight of energy.

Data Synthesis
The impact of climate change tends to cause non-stationary conditions in hydrological data. Therefore, we use the EEMD [14] in this study, and proposes a new way to synthesize and generate data that can solve the non-stationary data problems. For a given flow time series, the IMFs and residue, the EEMD decomposed result, would be replace by corresponding component in another period; next, the new combination of IMFs and residue are summed up according to the completeness of EMD [2]; then, the synthetic new flow time series will be produced. Thus, if there are many EEMD decomposed results with the same length but different period, we could synthesize a flow data set through permutation and combination of the corresponding IMFs and residues.
In this study, the historical 60-year (1956-2015) flow data of Tongtou streamflow gauge are segmented in groups of 10 years. Through the EEMD method divides the sequence into two parts: IMFs represent the inherent period sequence of the data, and the residue represents the non-period sequence of the data. We can go through the permutations and combinations of the IMFs to get n i sequences. Where n is the number of groups of data grouping; i is the number of IMFs by the group.
After completing the permutations and combinations of IMFs as shown in Figure 5, we can get the new synthesis data of equal length streamflow. In this study, two permutations and combinations methods are used to synthesize new 10-year flow series data. The method is described as follows:

1.
Regression analysis is performed on the different residue of each group to fit a new residue to represent the residue value. By adding the results of the IMFs permutation and combination with the representative values of the residue, we can get 10-year new synthetic flow data for the n i groups.

2.
Permutation and combination of the different remainders of each group with IMFs, we can get 10-year new synthetic flow data for the n i+1 groups.

Water Supply System Simulation of Hushan Reservoir
The simulation method is commonly used in the analysis of water resources systems. It conceptualizes the actual operation of the system to explore the use of water resources. In this study, we use the simulation method to simulate the water supply system and explore the difference between the newly synthesized flow data for each goal in the water supply system of Hushan reservoir.
First of all, the historical (1956-2015) flow data of Qingshui river were used to simulate accordance with the Hushan reservoir operation rules. Because the Hushan reservoir takes water from the Qingshui river through a water diversion tunnel, and there is some agricultural water right downstream of the Tongtou weir (as shown in Figure 6). Therefore, before taking water to the Hushan reservoir, Tongtou weir has to consider and satisfy the agricultural water downstream. The remaining flow is the amount of water that can be taken to the Hushan reservoir. The main thing of simulation is that the water taken into the Hushan reservoir should not affect the downstream agricultural water right.
This study will assume different scenarios to simulate the water supply system of Hushan reservoir. The main difference of these two scenarios is the object of water supply (as shown in Table  2). The local water supply system is shown in Figure 6, and the simulation program calculation process is shown in Figure 7.
In the past, the Linnei water purifying plant only supplied 1.2 × 10 5 m 3 water per day, and the deficit of demand was filled by pumping groundwater. After the completion of the Hushan reservoir, the situation of over pumping will be improved. The Hushan reservoir was originally intended to be used in conjunction with the Jiji weir to supply water for domestic demand in the Yunlin area. However, the conjunction operation of Hushan reservoir and Jiji weir will result in a low reservoir utilization rate of only 1.18 [25]. Therefore, we assume the Hushan reservoir operates independently in this study. The domestic demand in the Yunlin area is 2.8 × 10 5 /day, and the support demand for industrial is 3.0 × 10 5 /day. Based on the previous study, the operational regulations of reservoir were added [26], and the upper and lower rule curve is set at the water level of 190 m and 185 m for Hushan reservoir, respectively. In order to explore the situation of water resources utilization in the Yunlin area of the Hushan reservoir under different demand scenarios, the operation of reservoir should follow these principles: 1. Supply of the Hushan reservoir must not be less than the projected demand if the water surface elevation exceeds the upper curve. 2. Supply of the Hushan reservoir should satisfy 100% of the projected domestic demand, 90% of the projected industrial demand, and 50% of the projected irrigation demand if the water surface

Water Supply System Simulation of Hushan Reservoir
The simulation method is commonly used in the analysis of water resources systems. It conceptualizes the actual operation of the system to explore the use of water resources. In this study, we use the simulation method to simulate the water supply system and explore the difference between the newly synthesized flow data for each goal in the water supply system of Hushan reservoir.
First of all, the historical (1956-2015) flow data of Qingshui river were used to simulate accordance with the Hushan reservoir operation rules. Because the Hushan reservoir takes water from the Qingshui river through a water diversion tunnel, and there is some agricultural water right downstream of the Tongtou weir (as shown in Figure 6). Therefore, before taking water to the Hushan reservoir, Tongtou weir has to consider and satisfy the agricultural water downstream. The remaining flow is the amount of water that can be taken to the Hushan reservoir. The main thing of simulation is that the water taken into the Hushan reservoir should not affect the downstream agricultural water right.
Water 2020, 12, 927 9 of 20   This study will assume different scenarios to simulate the water supply system of Hushan reservoir. The main difference of these two scenarios is the object of water supply (as shown in Table 2). The local water supply system is shown in Figure 6, and the simulation program calculation process is shown in Figure 7.

Indices for Impact and Risk Assessment
To assess the impact and the risk under different scenarios, the shortage index (SI), satisfaction and reliability of the water supply were determined. Equations (10)-(16) define these indices.

Satisfaction:
2. Reliability: 3. Shortage Index: In the past, the Linnei water purifying plant only supplied 1.2 × 10 5 m 3 water per day, and the deficit of demand was filled by pumping groundwater. After the completion of the Hushan reservoir, the situation of over pumping will be improved. The Hushan reservoir was originally intended to be used in conjunction with the Jiji weir to supply water for domestic demand in the Yunlin area. However, the conjunction operation of Hushan reservoir and Jiji weir will result in a low reservoir utilization rate of only 1.18 [25]. Therefore, we assume the Hushan reservoir operates independently in this study. The domestic demand in the Yunlin area is 2.8 × 10 5 /day, and the support demand for industrial is 3.0 × 10 5 /day. Based on the previous study, the operational regulations of reservoir were added [26], and the upper and lower rule curve is set at the water level of 190 m and 185 m for Hushan reservoir, respectively. In order to explore the situation of water resources utilization in the Yunlin area of the Hushan reservoir under different demand scenarios, the operation of reservoir should follow these principles:

1.
Supply of the Hushan reservoir must not be less than the projected demand if the water surface elevation exceeds the upper curve.

2.
Supply of the Hushan reservoir should satisfy 100% of the projected domestic demand, 90% of the projected industrial demand, and 50% of the projected irrigation demand if the water surface elevation is between the upper and lower curve.

3.
Supply of the Hushan reservoir should satisfy 80% of the planning domestic demand and 0% of the projected irrigation and industrial demand if the water surface elevation is under the lower curve.

Indices for Impact and Risk Assessment
To assess the impact and the risk under different scenarios, the shortage index (SI), satisfaction and reliability of the water supply were determined. Equations (10)-(16) define these indices.

1.
Satisfaction: 2. Reliability: 3. Shortage Index: where Q d and Q sup are the quantities of water demand and supply, respectively; N is the given time scale of the simulation, and N sat is the number of days for which the demand has been satisfied. 4.
The average duration of water shortage events: where M is the number of occurrences of water shortage events, l j is the duration of each water shortage event, l n is the average duration of water shortage events. 6.
The average water shortage in the water shortage event (10 6 m 3 ): where d j is the total water shortage in each event, and d n is the average water shortage in the water shortage event. 7.
The daily average water shortage in the water shortage event (10 6 m 3 /day):

60-Year Flow Data Analysis with EEMD
We analyze the historical flow data from 1956 to 2015 in this study. Because the skewness of flow data will affect the decomposition process, the historical flow data is taken nature log transformation before decomposing by the EEMD. The 60-year historical data is decomposed by the EEMD method, and the decomposed IMFs are shown in Figure 8.
As can be seen from Figure 9, the energy value is mostly concentrated in IMF4, reaching 39% and the energy of IMF4 is higher than that of other IMFs. Corresponding to IMF4 in Figure 10, it can be found that the corresponding period is about 36 ten-day, which represents that the flow data has a strong annual (36 ten-day) period characteristic. In addition, the IMF1 and IMF2 still have a certain proportion of the energy percentage, which are 29% and 22% respectively, and their periods correspond to about 3 and 6 ten-day, respectively. Therefore, it can be seen that seasonal characteristics exists in the flow data.    Figure 8 shows that there are 10 IMFs and 1 residue. The termination condition of the EEMD is according to following literatures: Wu and Huang [14] show that the number of IMFs approaching log 2 n can be obtained by EEMD decomposition of the original data; in addition, Huang et al. [27] also explained that the IMF obtained after the data is decomposed will be between log 2 n − 1 and log 2 n. Therefore, the result that deals with 60-year daily flow by EEMD will introduce 10 IMFs and 1 residue. Among them, the oscillation of IMF 1 -IMF 10 gradually decrease from ± 3 to ± 0.008, but the range of the residue gradually increased from 1.73 to 2.09 and then decreased to 1.75, which is a monotonic function. From the decomposition, we can find that the residue part plays an important role, which often affects the trend of whole time series.
Historical flow data is regarded as a signal for a certain period of time. After decomposing the flow data through EEMD, the physical meaning of IMFs will be diagnosed. The period of each IMF is estimated according to zero up-cross analysis. However, the true energy value cannot be calculated, the square of the IMF value only be regarded as direct proportion with the true energy. In other words, the energy representation, the summation of square of each IMF, should be regarded as the 'weight' of oscillation for IMFs.
The values of the IMFs decomposed by the EEMD are squared and the sum is taken as the representative value of energy, as shown in Figure 9 and Table 3. After the zero up cross analysis method is used to average the IMF period of each group, the average period represented by the IMF can be obtained. IMF 10 is NAN because the period of the sequence cannot get after the zero up cross analysis method.

Analysis of 6 Groups of 10-Year Flow Data
The historical flow data (1956-2015) is divided into 6 decades, and the EEMD screening is implement to get the IMFs and the residue. The relationship between the period, energy, and flow data of different segment IMFs will be discussed. The period and energy percentage distribution of IMFs for each decade are shown in Table 4 and Figures 11 and 12.
The corresponding relationship between the energy percentage of each IMFs and the average period is shown in Figures 11 and 12. We can find the energy distribution concentrates in IMF1-IMF4, which correspond to an average period from 2.96 to 37.03 ten-days. In addition, it can be found that the energy percentage of IMF4 for the 6 groups has the most highest weight, with an average of about 20.5-44.5% compared with total weight of all IMFs, and it corresponds to an average period around 36 ten-days. Moreover, the minor IMFs, IMF1 and IMF2, are have around 28-36 % and 18-27 % energy of the total respectively, and their corresponding periods are about 3 and 6 ten-days. These results are similar to the 60-year flow analysis, as shown in Table 4.   As can be seen from Figure 9, the energy value is mostly concentrated in IMF 4 , reaching 39% and the energy of IMF 4 is higher than that of other IMFs. Corresponding to IMF 4 in Figure 10, it can be found that the corresponding period is about 36 ten-day, which represents that the flow data has a strong annual (36 ten-day) period characteristic. In addition, the IMF 1 and IMF 2 still have a certain proportion of the energy percentage, which are 29% and 22% respectively, and their periods correspond to about 3 and 6 ten-day, respectively. Therefore, it can be seen that seasonal characteristics exists in the flow data.

Analysis of 6 Groups of 10-Year Flow Data
The historical flow data (1956-2015) is divided into 6 decades, and the EEMD screening is implement to get the IMFs and the residue. The relationship between the period, energy, and flow data of different segment IMFs will be discussed. The period and energy percentage distribution of IMFs for each decade are shown in Table 4 and Figures 11 and 12.
The corresponding relationship between the energy percentage of each IMFs and the average period is shown in Figures 11 and 12. We can find the energy distribution concentrates in IMF1-IMF4, which correspond to an average period from 2.96 to 37.03 ten-days. In addition, it can be found that the energy percentage of IMF4 for the 6 groups has the most highest weight, with an average of about 20.5-44.5% compared with total weight of all IMFs, and it corresponds to an average period around 36 ten-days. Moreover, the minor IMFs, IMF1 and IMF2, are have around 28-36 % and 18-27 % energy of the total respectively, and their corresponding periods are about 3 and 6 ten-days. These results are similar to the 60-year flow analysis, as shown in Table 4.

Analysis of 6 Groups of 10-Year Flow Data
The historical flow data (1956-2015) is divided into 6 decades, and the EEMD screening is implement to get the IMFs and the residue. The relationship between the period, energy, and flow data of different segment IMFs will be discussed. The period and energy percentage distribution of IMFs for each decade are shown in Table 4 and Figures 11 and 12.

Analysis of 6 Groups of 10-Year Flow Data
The historical flow data (1956-2015) is divided into 6 decades, and the EEMD screening is implement to get the IMFs and the residue. The relationship between the period, energy, and flow data of different segment IMFs will be discussed. The period and energy percentage distribution of IMFs for each decade are shown in Table 4 and Figures 11 and 12.
The corresponding relationship between the energy percentage of each IMFs and the average period is shown in Figures 11 and 12. We can find the energy distribution concentrates in IMF1-IMF4, which correspond to an average period from 2.96 to 37.03 ten-days. In addition, it can be found that the energy percentage of IMF4 for the 6 groups has the most highest weight, with an average of about 20.5-44.5% compared with total weight of all IMFs, and it corresponds to an average period around 36 ten-days. Moreover, the minor IMFs, IMF1 and IMF2, are have around 28-36 % and 18-27 % energy of the total respectively, and their corresponding periods are about 3 and 6 ten-days. These results are similar to the 60-year flow analysis, as shown in Table 4.    The corresponding relationship between the energy percentage of each IMFs and the average period is shown in Figures 11 and 12. We can find the energy distribution concentrates in IMF 1 -IMF 4 , which correspond to an average period from 2.96 to 37.03 ten-days. In addition, it can be found that the energy percentage of IMF 4 for the 6 groups has the most highest weight, with an average of about 20.5-44.5% compared with total weight of all IMFs, and it corresponds to an average period around Water 2020, 12, 927 14 of 21 36 ten-days. Moreover, the minor IMFs, IMF 1 and IMF 2 , are have around 28-36% and 18-27% energy of the total respectively, and their corresponding periods are about 3 and 6 ten-days. These results are similar to the 60-year flow analysis, as shown in Table 4.

Data Synthesis
This section explains the synthesis of flow data. The historical flow data  is divided into six groups of 10-year segments in order. After decomposing through EEMD, multiple IMFs and residue corresponding to them can be obtained. The decadal residues are shown in Figure 13. Regression analysis of different residue in each group is fitted to a new residue representing the value.
This study used two methods to synthesize new decadal flow data as the basis for simulation research. Method (I): first permutations and combinations all IMFs to get n i (6 7 = 279,936) new IMF (IMF 1 + IMF 2 + . . . + IMF 7 ) set. Where n is the number of groups and i is the number of IMFs each group. After permutations and combinations of the IMFs, add a representative value of the residue, which is the regression of six residues. The representative value of the residue is shown in Figure 13, the bold dashed line. Regression analysis fits the representative value of the residue, and then adds this residue to the IMF of the 279,936 group to obtain the synthesis flow data of the 279,936 group for 10 years. Method (II): all the IMFs and the residue in each group are permutations and combinations to obtain a total of n i + 1 (6 7 + 1 = 1,679,616) new synthetic flow data. In the following, the new flow data synthesized from the two methods and historical flow data are compared with each other, and the differences are explored. In addition, this study will apply new flow data synthesized from two different methods to simulate Hushan reservoir water supply systems, and sequentially discuss the simulation results of different synthetic flow data in different water supply scenarios. Result of comparing the synthetic data with historical data divided by 10 years is shown in Table  5 and Figure 14. For the historical data, the flow in Qingshui river is significantly different between the wet and dry season; the flow from wet season, May to October, is significantly higher than November to April. The ratio of wet/dry is about 9:1. The highest flow occurs in August (62.09 m 3 /s) and the lowest flow occurs in January (2.05 m 3 /s). Compared the historical with the synthetic flow, the monthly distribution is similar between both Method I and Method II. The synthetic flow is concentrated in May-October, and the highest flow always occurs in August, but the lowest flow occurs in the dry season in January and February. In terms of the overall flow during the wet season, the average synthesized by Method I is significantly less than the current situation; in contrast, the average synthesized by Method II is much higher. In the dry season, the average synthesized by Method I and Method II is not much different from the current.
The distribution of the probability of exceedance is shown in the Figure 15 and 16. For the 10year synthesis, the overall time distribution is consistent with the historical one; furthermore, the synthesis keeps the distribution of wet and dry seasons as well. Through the probability of exceedance of 0.05, the extreme flow of Method II is greater than Method I. In contrast, while the probability of exceedance is 0.95, the flow data of Method II is much less than Method I. This will also lead to the opposite result of the average flow in the simulation of the water supply system.
Although the preceding says that, there are not much variance of the distribution trend between Result of comparing the synthetic data with historical data divided by 10 years is shown in Table 5 and Figure 14. For the historical data, the flow in Qingshui river is significantly different between the wet and dry season; the flow from wet season, May to October, is significantly higher than November to April. The ratio of wet/dry is about 9:1. The highest flow occurs in August (62.09 m 3 /s) and the lowest flow occurs in January (2.05 m 3 /s). Compared the historical with the synthetic flow, the monthly distribution is similar between both Method I and Method II. The synthetic flow is concentrated in May-October, and the highest flow always occurs in August, but the lowest flow occurs in the dry season in January and February. In terms of the overall flow during the wet season, the average synthesized by Method I is significantly less than the current situation; in contrast, the average synthesized by Method II is much higher. In the dry season, the average synthesized by Method I and Method II is not much different from the current.

Application of Synthesis Data
Next, we apply the synthesized flow data to simulate the water supply system according to the given two demand scenarios. Using the historical flow data, Method I and the Method II synthetic flow data are used for simulation, and the simulation results are shown in the Table 6 and Figure 17.
First, we apply the historical flow to simulate the water supply system. For the case of no supporting industrial water in Scenario I, the annual water supply and shortage are 101.93 × 10 6 m 3 and 0.27 × 10 6 m 3 , respectively. The water shortage index (SI) is only 0.002. While the water shortage occurring, the average daily water shortage (ADWS) is only 0.03 × 10 6 m 3 . The preceding indices show that if the Hushan reservoir only supplies water for domestic demand, the water supply situation The distribution of the probability of exceedance is shown in the Figures 15 and 16. For the 10-year synthesis, the overall time distribution is consistent with the historical one; furthermore, the synthesis keeps the distribution of wet and dry seasons as well. Through the probability of exceedance of 0.05, the extreme flow of Method II is greater than Method I. In contrast, while the probability of exceedance is 0.95, the flow data of Method II is much less than Method I. This will also lead to the opposite result of the average flow in the simulation of the water supply system.

Application of Synthesis Data
Next, we apply the synthesized flow data to simulate the water supply system according to the given two demand scenarios. Using the historical flow data, Method I and the Method II synthetic

Application of Synthesis Data
Next, we apply the synthesized flow data to simulate the water supply system according to the given two demand scenarios. Using the historical flow data, Method I and the Method II synthetic flow data are used for simulation, and the simulation results are shown in the Table 6 and Figure 17. Although the preceding says that, there are not much variance of the distribution trend between the 10-year flow data synthesized by these two different methods. However, the annual average between two synthesis methods is more significant. The flow data synthesized by Method II has an annual average flow of 24.34 m 3 /s, which is more than the annual average of Method I, 13.71 m 3 /s. The water supply system simulation will be based on these three flow sequences, historical data, Method I and Method II synthetic data, to simulate the water supply system of Hushan reservoir.

Application of Synthesis Data
Next, we apply the synthesized flow data to simulate the water supply system according to the given two demand scenarios. Using the historical flow data, Method I and the Method II synthetic flow data are used for simulation, and the simulation results are shown in the Table 6 and Figure 17. Although the SI only increases to 0.309, and the ADWS only increases to 0.06 × 10 6 m 3 while the water shortage occurred. The reliability of the water supply reaches 0 already at the 9th and 11th ten-day, as depicted in Figure 17. The result show that the water supply risk in demand Scenario II is relatively severe than Scenario I. Because the simulation of demand Scenario II indicates severe impact on domestic, we set the domestic SI = 0.1, and try to find the capability of industrial supporting for Hushan reservoir. We found that Hushan reservoir can support 100,000 m 3 /day of industrial water, and the annual water shortage has decreased to 2.18 × 10 6 m 3 . The number of average consecutive dry days (ACDD) for domestic has been greatly reduced to 45 days, and the reliability of the water supply is also raised to 0.4-0.5. In the following analyses, we recommend the daily supporting amount for industry from the Hushan reservoir should be 100,000 m 3 as basis.  We simulate the water supply system with the historical flow as applied to method I (279,936 groups) and method II (1,679,616 groups) for newly synthesized flow data. The results are compared based on the ensemble average (as shown in Table 7). Scenario I without considers supporting industrial water, the current annual water supply is 101.93 × 10 6 m 3 , the annual water shortage is only 0.27 × 10 6 m 3 , and the SI is only 0.002. While the water shortage occurring, the ADWS is only 0.03 × 10 6 m 3 .
According to the water supply simulation results of 279,936 sets of flow data in Method I, the annual water supply slightly decreased to 98.25 × 10 6 m 3 , the annual water shortage increased to 3.59 × 10 6 m 3 , and the SI increased significantly to 0.668. While the water shortage occurring, the ADWS is increased slightly to 0.05 × 10 6 m 3 . However, Figure 18 shows that whether the reliability is lower First, we apply the historical flow to simulate the water supply system. For the case of no supporting industrial water in Scenario I, the annual water supply and shortage are 101.93 × 10 6 m 3 and 0.27 × 10 6 m 3 , respectively. The water shortage index (SI) is only 0.002. While the water shortage occurring, the average daily water shortage (ADWS) is only 0.03 × 10 6 m 3 . The preceding indices show that if the Hushan reservoir only supplies water for domestic demand, the water supply situation will be very stable and shortage will not be easy to occur. In the case of supporting industrial water in Scenario II. If the daily support for industrial water is 300,000 m 3 , the annual water supply will reach 165.14 × 10 6 m 3 , and the annual water shortage will also increase significantly to 9.74 × 10 6 m 3 . Although the SI only increases to 0.309, and the ADWS only increases to 0.06 × 10 6 m 3 while the water shortage occurred. The reliability of the water supply reaches 0 already at the 9th and 11th ten-day, as depicted in Figure 17. The result show that the water supply risk in demand Scenario II is relatively severe than Scenario I.
Because the simulation of demand Scenario II indicates severe impact on domestic, we set the domestic SI = 0.1, and try to find the capability of industrial supporting for Hushan reservoir. We found that Hushan reservoir can support 100,000 m 3 /day of industrial water, and the annual water shortage has decreased to 2.18 × 10 6 m 3 . The number of average consecutive dry days (ACDD) for domestic has been greatly reduced to 45 days, and the reliability of the water supply is also raised to 0.4-0.5. In the following analyses, we recommend the daily supporting amount for industry from the Hushan reservoir should be 100,000 m 3 as basis.
We simulate the water supply system with the historical flow as applied to method I (279,936 groups) and method II (1,679,616 groups) for newly synthesized flow data. The results are compared based on the ensemble average (as shown in Table 7). Scenario I without considers supporting industrial water, the current annual water supply is 101.93 × 10 6 m 3 , the annual water shortage is only 0.27 × 10 6 m 3 , and the SI is only 0.002. While the water shortage occurring, the ADWS is only 0.03 × 10 6 m 3 .
According to the water supply simulation results of 279,936 sets of flow data in Method I, the annual water supply slightly decreased to 98.25 × 10 6 m 3 , the annual water shortage increased to 3.59 × 10 6 m 3 , and the SI increased significantly to 0.668. While the water shortage occurring, the ADWS is increased slightly to 0.05 × 10 6 m 3 . However, Figure 18 shows that whether the reliability is lower than the current simulation result, they are still above 0.5, the lowest satisfaction is 0.87, which is generally acceptable.
According to the water supply simulation results of 1,679,616 sets of flow data in Method II. The annual water supply decreases to 95.52 × 10 6 m 3 , the annual water shortage has increased significantly to 6.68 × 10 6 m 3 . The SI increased significantly to 1.96 and the ADWS increased to 0.07 × 10 6 m 3 as well. Figure 18 shows that the reliability is similar to the method I in dry season. However, Method II's reliability is lower than method I value in the wet season. The satisfaction is generally lower than the current situation and the method I, but the minimum is still 0.84. The water supply situation of Method II is still acceptable. There is something odd: why the average annual flow of Method II is significantly more than the current situation and the average flow of Method I, but the simulation results are worse than the other cases? The main reason is that the extreme flow of the Method II often occurs; excessive flow cannot be used efficiently and there are many low flows during wet season. Therefore, the simulation result of Method II is the worst.  Finally, we simulate the water supply system using the historical flow with Method I and Method II's synthesis flow data. The results are compared to the ensemble average (as shown in Table  8). Under the condition of supporting 100,000 m 3 of industrial demand daily, the current annual water supply is 128.30 × 10 6 m 3 , the annual water shortage is 2.18 × 10 6 m 3 and the SI is 0.1. While water shortage occurs, the average daily water shortage is 0.05 × 10 6 m 3 .
According to the simulation result of water supply in Method I, the annual water supply decreased slightly to 116.82 × 10 6 m 3 , the annual water shortage increased to 11 × 10 6 m 3 , and the SI increased significantly to 1.203. While the water shortage occurring, the ADWS is still maintained at 0.05 × 10 6 m 3 . Figure 19 shows that the reliability decrease more earlier than the current situation, the lowest monthly reliability reaches 0.24 in March, and recoveries in May. The lowest fulfillment of Finally, we simulate the water supply system using the historical flow with Method I and Method II's synthesis flow data. The results are compared to the ensemble average (as shown in Table 8). Under the condition of supporting 100,000 m 3 of industrial demand daily, the current annual water supply is 128.30 × 10 6 m 3 , the annual water shortage is 2.18 × 10 6 m 3 and the SI is 0.1. While water shortage occurs, the average daily water shortage is 0.05 × 10 6 m 3 .
According to the simulation result of water supply in Method I, the annual water supply decreased slightly to 116.82 × 10 6 m 3 , the annual water shortage increased to 11 × 10 6 m 3 , and the SI increased significantly to 1.203. While the water shortage occurring, the ADWS is still maintained at 0.05 × 10 6 m 3 . Figure 19 shows that the reliability decrease more earlier than the current situation, the lowest monthly reliability reaches 0.24 in March, and recoveries in May. The lowest fulfillment of demand still occurs in March (0.74). In terms of dry season, the risk of water supply shortage is much worse than the current situation.
According to the simulation result of water supply in Method II, the annual water supply decreased to 113.82 × 10 6 m 3 , the annual water shortage increased sharply to 13.46 × 10 6 m 3 , and the SI increased significantly to 2.12. While occurs the water shortage event, the ADWS also increased to 0.07 × 10 6 m 3 . Figure 19 depicts that the reliability is similar to but lower than Method I in the dry season. Compared with the current situation, the water shortage occurs earlier, but it is lower than the Method I in the wet season. The satisfaction is generally lower than the current situation and the Method I, but similar to Method I. The satisfaction rises rapidly in May, which is higher than the current simulation results. However, the satisfaction is decreasing from August to October in wet season, it seems worse than the others.
In summary, no matter what kind of water supply scenario is applied, the current results are better than Method I and Method II. Although the simulation result of Method I is similar to Method II, the ensemble result of Method I is slightly better than the Method II. However, Method II synthesis annual average flow is significantly more than the current flow and the Method I flow, but the simulation results are worse than other cases. The main reason is that the extreme flow of the method II often occurs. Excessive flow cannot be used efficiently and there are many low flows in wet season. Therefore, the simulation result of Method II is the worst. season, it seems worse than the others. In summary, no matter what kind of water supply scenario is applied, the current results are better than Method I and Method II. Although the simulation result of Method I is similar to Method II, the ensemble result of MethodⅠis slightly better than the Method II. However, Method II synthesis annual average flow is significantly more than the current flow and the Method I flow, but the simulation results are worse than other cases. The main reason is that the extreme flow of the method II often occurs. Excessive flow cannot be used efficiently and there are many low flows in wet season. Therefore, the simulation result of Method II is the worst.

Conclusions
The IMFs and residue of the Tongtou weir flow data can be obtained by EEMD decomposition. This method can effectively solve the stationary or nonstationary impacts of hydrological data due to climate change. In this study, two methods are used to synthesize new 10-year flow series data. Method I uses the IMFs segmented data to reorganize the periodic characteristics and sums them up with the residue representation; in contrast, Method II uses the permutation and combination of all

Conclusions
The IMFs and residue of the Tongtou weir flow data can be obtained by EEMD decomposition. This method can effectively solve the stationary or nonstationary impacts of hydrological data due to climate change. In this study, two methods are used to synthesize new 10-year flow series data. Method I uses the IMFs segmented data to reorganize the periodic characteristics and sums them up with the residue representation; in contrast, Method II uses the permutation and combination of all IMFs and residues. The new synthetic flow data will provide the input of simulation for the Hushan Reservoir water supply system.
In the energy distribution of IMFs, most of the energy is concentrated between IMF 1 and IMF 4 , which IMF 4 has the largest weight (39%), and the corresponding period is about 36 ten-day periods (1 year). However, IMF 1 still has a certain proportion of energy percentage, the corresponding period is about 3 ten-day periods (1 month). This reflects that the flow of Tongtou weir is not only affected by a 1-year periodic sequence, but influenced by the monthly rainfall event. Comparing the current situation with the new data generated by EEMD, the generated data is similar to the current flow distribution, but the flow distribution during the wet season is significantly different (±0.63 m 3 /s), and the flow difference during the dry season is insignificant (±0.78 m 3 /s).
For the Scenario I without considers supporting industrial water, the simulation results of the generation data of Method I show that compared with the current situation, the SI has increased significantly from 0.002 to 0.668, and the number of ACDD of water shortage has also increased significantly from 9 to 51.2 days. The simulation results of Method II show that the SI is more severe than the current situation and Method I, it reaches 1.96. The number of ACDD increased slightly to 59.28 days. In addition, it can be found from the reliability and satisfaction Figure that the simulation results of the synthesis flow in this study will occur water shortages earlier than the current situation.
Scenario II considers supporting industrial water usage, the simulation results of Method I show that the SI has increased significantly to 1.023, and the number of ACDD of water shortage has also increased significantly to 111 days. The simulation results of Method II show that the SI is more severe to 2.12. The number of ACDD increased slightly to 113.82 days. In addition, it can be found from the reliability and satisfaction figure that no matter in which way the synthesis flow is simulated, supporting industrial water will seriously affect the equity of domestic water.
In this study, the multiplication data obtained through permutation and combination after EEMD can show the flow distribution characteristics of the Tongtou weir. In the future, the hydrological factors, such as temperature and rainfall in the same area, can be integrated to explore the correlation between the overall hydrological factors.