**1. 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.

EMD has been proposed [2] as an adaptive time-frequency data analysis method. In this study, the EMD method and the decomposed IMFs were used for data analysis and synthesis. Many studies use EMD to analyze different time series and decompose the implied fluctuation periods [3–6]. In recent years, there are many carefully studies of nonlinear and non-stationary for the hydrological field: Lee and Ouarda [7,8] proposed model reproduces Nonstationary oscillation (NSO) processes by utilizing EMD and nonparametric simulation techniques; Molla et al. [9,10] show that the EMD successively extracts the IMFs with the highest local frequencies in a recursive way, which yields effectively a set low-pass filters based entirely on the properties exhibited by the data. In addition, there are many studies based on the EMD method to predict and estimate streamflow data: Huang et al. [11] and Meng et al. [12] developed a modified EMD-based support vector machine for non-stationary streamflow prediction; Zhang et al. [13] improved the efficiency of a prediction approach which can be used as a general technique for non-stationary time series forecasting.

However, the mode mixing problem may occur during the EMD. To avoid mode mixing, Wu and Huang [14] proposed a new noise-assisted data analysis method, the Ensemble EMD (EEMD), 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. Zhao and Chen [15] and Zhang et al. [16] used the EEMD and other methods to build a hybrid model for annual runoff time series forecasting. Many of EEMD applications in hydrology focus on the diagnostic analysis of historical data, such as features of temperature variation trends [17], forward prediction of runoff data in data-scarce basins [18], and forecast monthly reservoir inflow [19].

The analysis methods of water resources allocation system are mainly divided into two categories: simulation method and optimization method. Yeh [20] and Wurbs [21] introduced the theory development of the simulation model and optimization model in reservoir system analysis. Chang et al. [22] present the usefulness of the neural network in deriving general operating policies for a multi-reservoir system. An optimization-based approach is typically required for the optimal operation of a reservoir system, in order to obtain the optimal solutions to support the decision-making process [23]. In this study, we use the simulation method for the analysis of water resources systems. Lee and Huang [24] evaluated the impact of climate change on the water supply of the reservoir in central Taiwan.

The Hushan reservoir is a new water supply resource system in central Taiwan (Figure 1). Since some problems exist in the water supply system of the Zhuoshui river, such as its high concentration of sediment, and low efficiency for the Hushan reservoir while it considers a joint operation with Jiji weir [25]. Therefore, Chu and Huang [26] recommends the Hushan reservoir should operate independently. Under the assumption that simulation results show that the Hushan reservoir is capable of supplying more than 100,000 m3 per day (up to 380,000 m3/day) and the reservoir utilization rate can reach 2.41 from 1.18 while Shortage Index (SI) = 1.

This study proposes a novel procedure, which is based on EEMD and its derivation (IMF and residue) as an implement for synthesizing non-stationary data, which will be applied in reservoir inflow and impact assessment for water supply. This study will continue the research results of Huang et al. [27]. Huang et al. [27] proposed a surface temperature and rainfall synthesis method base on EMD. The method utilizes the recombination of the IMF of the segmented data, as well as the characteristics of the residuals to generate the data. There is a plan to use the IMFs and residues obtained, after EMD decomposition to permutation and combination, to generate new flow data. There is also a plan to try to apply the newly synthetic flow data to simulate the water supply system of the Hushan reservoir (as shown in Figure 1). Finally, several evaluation indicators are used to evaluate water supply system simulation results.

**Figure 1.** Research flowchart. EMD, empirical mode decomposition components.

#### **2. Materials and Method**

#### *2.1. Study Area*

#### 2.1.1. 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 <sup>×</sup> 106 m3. 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 m3/s and 1.00 m3/s, respectively. The watershed of Hushan reservoir is only 6.58 km2 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.

**Figure 2.** Location of study area.


**Table 1.** Elevation-area-capacity relationship of Hushan reservoir.

#### 2.1.2. 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 m3/s, with a gravity type.

#### *2.2. Research Method*

#### 2.2.1. 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:


*Water* **2020**, *12*, 927

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]:


Finally, original time series x(t) can be denoted as sum of IMFs ci(t) and residue r(t).

$$\mathbf{x}(\mathbf{t}) = \sum\_{i=1}^{n} \mathbf{c}\_{i}(\mathbf{t}) + \mathbf{r}(\mathbf{t}) \tag{1}$$

where n, ci(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].

**Figure 3.** Results of EMD (empirical mode decomposition) method decomposition.

#### 2.2.2. Ensemble Empirical Mode Decomposition

EEMD, an improvement of EMD, was proposed by Wu and Huang [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:

$$\mathbf{X(t)} = \mathbf{x(t)} + \mathbf{w(t)}\tag{2}$$


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]

$$
\varepsilon\_n = \frac{\varepsilon}{\sqrt{\mathcal{N}}} \tag{3}
$$

where N is the number of ensemble members, ε is the amplitude of the added noise and ε<sup>n</sup> 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.

#### 2.2.3. 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]

$$
\mathfrak{n}\_{\text{li}} \times \mathfrak{n}\_{\text{li}+1} < 0 ; \mathfrak{n}\_{\text{li}+1} > 0 \tag{4}
$$

where η<sup>i</sup> 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

$$
\eta\_{\text{li}-1} < \eta\_{\text{li}} \cap \eta\_{\text{li}} > \eta\_{\text{li}+1} \text{ (highest point)}\tag{5}
$$

$$
\mathfrak{l}\,\,\eta\_{\mathrm{i}-1} > \eta\_{\mathrm{i}} \cap \eta\_{\mathrm{i}} < \eta\_{\mathrm{i}+1} \,\,\big(\mathrm{lowest}\,\,\mathrm{point}\,\big)\tag{6}
$$

**Figure 4.** Show of the zero up-cross method [30]. H, The wave high; nmax: The average value of the highest point (peak).

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)

$$\text{amplitude} : (\mathfrak{n}\_{\text{max}} + \mathfrak{n}\_{\text{min}})/2 \tag{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

$$\mathbf{E}\_{\mathbf{i}} = \sum\_{\mathbf{i}=1}^{n} \mathbf{x}\_{\mathbf{i}} \,^2 \tag{8}$$

$$\text{W}\_{\text{i}} = \frac{\text{E}\_{\text{i}}}{\sum\_{i=1}^{n} \text{E}\_{\text{i}}} \times 100\% \tag{9}$$

where xj <sup>2</sup> is the square value at each time point of each group of IMF; Ei is the total energy; Wi is weight of energy.

#### 2.2.4. 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 ni 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:


**Figure 5.** Graphs of the data synthesis process. h: segment.

#### 2.2.5. 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.

**Figure 6.** Water supply system of Hushan reservoir.

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.


**Table 2.** Scenario description.

**Figure 7.** Program simulation flowchart.

In the past, the Linnei water purifying plant only supplied 1.2 <sup>×</sup> 105 m<sup>3</sup> 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 <sup>×</sup> 105/day, and the support demand for industrial is 3.0 <sup>×</sup> <sup>10</sup>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:


#### 2.2.6. 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.

*Water* **2020**, *12*, 927

1. Satisfaction:

$$\text{Satisfaction} = \frac{\text{Q}\_{\text{sup}}}{\text{Q}\_{\text{d}}} \times 100\% \tag{10}$$

2. Reliability:

$$\text{Reliability} = \frac{\text{N}\_{\text{sat}}}{\text{N}} \times 100\% \tag{11}$$

3. Shortage Index:

$$\text{SI} = \frac{100}{\text{N}} \sum\_{t=1}^{\text{N}} \left( \frac{|\text{Qd} - \text{Q}\_{\text{sup}}|}{\text{Q}\_{\text{d}}} \right)^{2} \tag{12}$$

where Qd and Qsup are the quantities of water demand and supply, respectively; N is the given time scale of the simulation, and Nsat is the number of days for which the demand has been satisfied.

4. Reservoir Efficiency:

$$\text{RE} = \frac{\text{Annual actual water supply}}{\text{reservoir effective capacity}} \tag{13}$$

5. The average duration of water shortage events:

$$\mathbf{\tilde{l}}\_{\rm n} = \frac{1}{\mathbf{M}} \sum\_{j=1}^{\mathbf{M}} \mathbf{l}\_{j} \tag{14}$$

where M is the number of occurrences of water shortage events, lj is the duration of each water shortage event, ln is the average duration of water shortage events.

6. The average water shortage in the water shortage event (106 m3):

$$\overline{\mathbf{d}}\_{\mathbf{n}} = \frac{1}{\mathbf{M}} \sum\_{\mathbf{j}=\mathbf{l}}^{\mathbf{M}} \mathbf{d}\_{\mathbf{j}} \tag{15}$$

where dj is the total water shortage in each event, and dn is the average water shortage in the water shortage event.

7. The daily average water shortage in the water shortage event (10<sup>6</sup> m3/day):

$$\frac{\overline{\mathbf{d}\_n}}{\overline{\mathbf{l}}\_n} \tag{16}$$

#### **3. Results and Discussion**

#### *3.1. 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.

**Figure 8.** Decomposition of 60-year flow data by EMD.

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 log2 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 log2 n − 1 and log2 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 IMF1–IMF10 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. IMF10 is NAN because the period of the sequence cannot get after the zero up cross analysis method.

**Figure 9.** Weight of each IMF for a 60-year flow.


**Table 3.** Results of 60-year data by EMD with energy and period.

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 10.** Average period of each IMF for a 60-year flow.

#### *3.2. 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.


**Table 4.** Results of decadal data by EMD with energy and period.

**Figure 11.** Weight of each IMF for a 10-year flow.

**Figure 12.** Average period of each IMF for a 10-year flow.

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.

#### *3.3. Data Synthesis*

This section explains the synthesis of flow data. The historical flow data (1956–2015) 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 ni (6<sup>7</sup> = 279,936) new IMF (IMF1 + IMF2 + ... + IMF7) 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<sup>i</sup> <sup>+</sup> <sup>1</sup> (67 <sup>+</sup> <sup>1</sup> = 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.

**Figure 13.** Flow data residues and their representative equation. CMS, m3/s; r\*(t): representative values of the residue.

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 m3/s) and the lowest flow occurs in January (2.05 m3/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.


**Table 5.** Comparison of historical and synthesis flow data.

**Figure 14.** Comparisons of 10-year flow data.

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.

**Figure 15.** Comparisons of 10-year flow data with Method I. P(X > x): probability of exceedance.

**Figure 16.** Comparisons of 10-year flow data with Method II.

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 m3/s, which is more than the annual average of Method I, 13.71 m3/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.

### *3.4. 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.



Note: CMD, m3/day; ACDD, average consecutive dry days; ADWS, average daily water shortage.

**Figure 17.** Ten-day reliability and satisfaction of the water supply in present each scenario. CMD, m3/day.

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 <sup>×</sup> 10<sup>6</sup> m3 and 0.27 <sup>×</sup> <sup>10</sup><sup>6</sup> m3, 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 <sup>×</sup> <sup>10</sup><sup>6</sup> <sup>m</sup>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 m3, the annual water supply will

reach 165.14 <sup>×</sup> <sup>10</sup><sup>6</sup> <sup>m</sup>3, and the annual water shortage will also increase significantly to 9.74 <sup>×</sup> 106 <sup>m</sup>3. Although the SI only increases to 0.309, and the ADWS only increases to 0.06 <sup>×</sup> <sup>10</sup><sup>6</sup> <sup>m</sup><sup>3</sup> 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 m3/day of industrial water, and the annual water shortage has decreased to 2.18 <sup>×</sup> 106 m3. 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 m3 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 <sup>×</sup> 106 m3, the annual water shortage is only 0.27 <sup>×</sup> 10<sup>6</sup> m3, and the SI is only 0.002. While the water shortage occurring, the ADWS is only 0.03 <sup>×</sup> 106 <sup>m</sup>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 <sup>×</sup> 106 m3, the annual water shortage increased to 3.59 <sup>×</sup> 10<sup>6</sup> m3, and the SI increased significantly to 0.668. While the water shortage occurring, the ADWS is increased slightly to 0.05 <sup>×</sup> 106 <sup>m</sup>3. However, Figure <sup>18</sup> 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 <sup>×</sup> 106 m3, the annual water shortage has increased significantly to 6.68 <sup>×</sup> 106 m3. The SI increased significantly to 1.96 and the ADWS increased to 0.07 <sup>×</sup> <sup>10</sup><sup>6</sup> m3 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.


**Table 7.** Comparison of historical and synthesis flow data. SI, Shortage Index.

**Figure 18.** Monthly (**a**) reliability and (**b**) satisfaction of the water supply for Scenario I.

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 m3 of industrial demand daily, the current annual water supply is 128.30 <sup>×</sup> 106 m3, the annual water shortage is 2.18 <sup>×</sup> 106 m3 and the SI is 0.1. While water shortage occurs, the average daily water shortage is 0.05 <sup>×</sup> 106 m3.

According to the simulation result of water supply in Method I, the annual water supply decreased slightly to 116.82 <sup>×</sup> 10<sup>6</sup> m3, the annual water shortage increased to 11 <sup>×</sup> 106 m3, and the SI increased significantly to 1.203. While the water shortage occurring, the ADWS is still maintained at 0.05 <sup>×</sup> <sup>10</sup><sup>6</sup> m3. 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 <sup>×</sup> 106 <sup>m</sup>3, the annual water shortage increased sharply to 13.46 <sup>×</sup> <sup>10</sup><sup>6</sup> m3, and the SI increased significantly to 2.12. While occurs the water shortage event, the ADWS also increased to 0.07 <sup>×</sup> <sup>10</sup><sup>6</sup> <sup>m</sup>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.


**Table 8.** Comparison of historical and synthesis flow data.

**Figure 19.** Monthly (**a**) reliability and (**b**) satisfaction of the water supply for Scenario II.

#### **4. 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 IMF1 and IMF4, which IMF4 has the largest weight (39%), and the corresponding period is about 36 ten-day periods (1 year). However, IMF1 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 m3/s), and the flow difference during the dry season is insignificant (±0.78 m3/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.

**Author Contributions:** T.-Y.C. analyzed the data and drafted the manuscript; W.-C.H. revised the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** This research was sponsored by the Ministry of Science and Technology in Taiwan (MOST 106-2625-M-019-001).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
