Calculation of Comprehensive Ecological Flow with Weighted Multiple Methods Considering Hydrological Alteration

Instream ecological flow is an essential determinant of river health. Intra- and interannual distribution characteristics of runoff have been altered to different degrees by dam construction. Historical runoff series with alterations, as basic data for ecological flow calculation, provide minimal instream hydrological process information, which affects the credibility of calculation results. Considering the influence of the alterations in runoff series on ecological flow calculation, the Gini coefficient (GI) is introduced to study the evenness degrees of the intra-annual runoff distribution of four hydrological stations located in the Naolihe basin of the Sanjiang Plain. The hydrological alteration diagnosis system is used to examine the alteration points in the GI series of each hydrological station for selecting reasonable subsequences. Based on the selected subsequences, the ecological flow of each station is calculated using three hydrological methods, and the comprehensive ecological flow is calculated using weighted calculation results from the three hydrological methods. The study results show that ecological flow and natural flow have similar processes with two peaks occurring in the process in May and August, respectively. Also, dams decrease the ecological water requirement damage frequency in dry seasons, but overuse of water resources increases the ecological water requirement damage frequency in flood seasons.


Introduction
Riverine ecosystems play a significant role in the survival and development of human society, so protecting riverine ecosystems is vital. Ecological flow is the runoff that maintains the riverine ecosystem health and biodiversity [1][2][3]. Studies on ecological flow have developed rapidly [4]. Many methods have been presented to calculate ecological flow [5,6], of which hydrological methods are widely used due to the simplicity of the data requirements and low calculation costs. A common hypothesis states that the ecological regime becomes well-adjusted to the hydrological regime, and the ecological requirements of river-living organisms have been satisfied from a hydrological standpoint [7]. The hydrological method calculates ecological flow by analyzing the frequency distribution characteristics of historical runoff series. Some representative hydrological methods, such as the Tennant method [8], NGPRP (Northern Great Plains Resource Program) method [9], minimum monthly mean flow method [10], month-by-month minimum ecological flow method [11], flow duration curve method [12], and the 7Q10 (seven-day low flow with a 10-year recurrence interval) Four hydrological stations are present in the Naolihe basin: Baoqing, Baoan, Caizuizi, and Hongqiling, as shown in Figure 1. The Baoqing hydrological station is located in the town of Baoqing, with a drainage area of 3689 km 2 . The Caizuizi hydrological station is located at the trunk stream of Naolihe River, in the town of Sanlitun, and has a control area of 20,896 km 2 . The Baoan hydrological station is located in the town of Youyi, having a drainage area of 1344 km 2 . The Hongqiling hydrological station is located on the Hongqiling farm. Its drainage area is 1147 km 2 . Two leading reservoirs are found in the Naolihe basin: the Longtouqiao Reservoir and the Hamatong Reservoir, shown in Figure 1. The first reclamation in the Naolihe basin began in 1956, and no agricultural activities were performed in the study area prior to 1956. Therefore, 1956 was chosen as the start year for this study. The historical daily flow data, except for December to March in the next year due to the freezing period, is from 1956 to 2012 for the four hydrological stations in the Naolihe Basin for this study. Location of the study area is shown in Figure 1. Four hydrological stations are present in the Naolihe basin: Baoqing, Baoan, Caizuizi, and Hongqiling, as shown in Figure 1. The Baoqing hydrological station is located in the town of Baoqing, with a drainage area of 3689 km 2 . The Caizuizi hydrological station is located at the trunk stream of Naolihe River, in the town of Sanlitun, and has a control area of 20,896 km 2 . The Baoan hydrological station is located in the town of Youyi, having a drainage area of 1344 km 2 . The Hongqiling hydrological station is located on the Hongqiling farm. Its drainage area is 1147 km 2 . Two leading reservoirs are found in the Naolihe basin: the Longtouqiao Reservoir and the Hamatong Reservoir, shown in Figure 1. The first reclamation in the Naolihe basin began in 1956, and no agricultural activities were performed in the study area prior to 1956. Therefore, 1956 was chosen as the start year for this study. The historical daily flow data, except for December to March in the next year due to the freezing period, is from 1956 to 2012 for the four hydrological stations in the Naolihe Basin for this study. Location of the study area is shown in Figure 1.

Methodology
The comprehensive ecological flow was calculated using the weighting mean method. The flow chart of the study process is shown in Figure 2.

Methodology
The comprehensive ecological flow was calculated using the weighting mean method. The flow chart of the study process is shown in Figure 2.

Evaluation of the Intra-Annual Runoff Distribution Degrees of Evenness
The Gini coefficient is an index used to quantify and evaluate the degrees of distribution evenness. In this paper, it is used to evaluate intra-annual runoff distribution evenness degrees. The GI of runoff series is calculated as follows: (1) Group and sort historical monthly flow data. Divide the 12 monthly flows in the ith year into a group, and sort the monthly flows into a group in ascending order. (2) Accumulate the ascending runoff data in each group: where k is the ranking order.

Evaluation of the Intra-Annual Runoff Distribution Degrees of Evenness
The Gini coefficient is an index used to quantify and evaluate the degrees of distribution evenness. In this paper, it is used to evaluate intra-annual runoff distribution evenness degrees. The GI of runoff series is calculated as follows: (1) Group and sort historical monthly flow data. Divide the 12 monthly flows in the ith year into a group, and sort the monthly flows into a group in ascending order. The ascending runoff series of the ith year is {r i1 , r i2 , . . . , r i12 }, where i represents the ith year, r i1 < r i2 < . . . < r i12 . (2) Accumulate the ascending runoff data in each group: where k is the ranking order. (3) Draw a Lorenz curve of the annual distribution for monthly flow. Take k/12 as the abscissa, R ik /R i12 as the ordinate, and draw a Lorenz curve of the annual distribution for monthly flow. The Lorenz curve is shown in Figure 3. (4) The ith year runoff Gini coefficient GI i is computed as: where S Ai and S Bi are the area acreages of A i and B i , respectively.
The calculation details of GI series in each station are shown in Appendix A.
2018, 10, x Draw a Lorenz curve of the annual distribution for monthly flow. Take k/12 as the abs 12 / ik i R R as the ordinate, and draw a Lorenz curve of the annual distribution for monthly The Lorenz curve is shown in Figure 3.
The ith year runoff Gini coefficient GIi is computed as: where SAi and SBi are the area acreages of Ai and Bi, respectively. The calculation details of GI series in each station are shown in Appendix A. The smaller the obtained value, the more even the runoff distribution in a year. Conversel r the obtained value, the more uneven the runoff distribution in a year.

. Hydrological Alteration Diagnosis System
A single hydrological alteration diagnosis method occasionally produces unreliable result iple hydrological alteration diagnosis methods produce different results. The hydrol ation diagnosis system synthetically considers the diagnosis results from multiple diag ods and examines the alteration form and alteration point holistically. The steps we foll lteration diagnosis are as follows. Firstly, we used the Hurst coefficient [25] method an ing average method to form a primary diagnosis and judge whether or not the series conta ation. If so, then various examination methods were used to conduct a detailed diag ding three trend diagnosis methods (the correlation coefficient method, Spearman lation method, and Kendall rank correlation method) and 11 jump diagnosis methods (the inan method, rank test, slide F test, R/S method, Mann-Kendall method, Bayesian method diagnosis results were also classified into two types, trend results and jump results, an  The smaller the obtained value, the more even the runoff distribution in a year. Conversely, the larger the obtained value, the more uneven the runoff distribution in a year.

Hydrological Alteration Diagnosis System
A single hydrological alteration diagnosis method occasionally produces unreliable results, and multiple hydrological alteration diagnosis methods produce different results. The hydrological alteration diagnosis system synthetically considers the diagnosis results from multiple diagnosis methods and examines the alteration form and alteration point holistically. The steps we followed for alteration diagnosis are as follows. Firstly, we used the Hurst coefficient [25] method and the moving average method to form a primary diagnosis and judge whether or not the series contains an alteration. If so, then various examination methods were used to conduct a detailed diagnosis, including three trend diagnosis methods (the correlation coefficient method, Spearman rank correlation method, and Kendall rank correlation method) and 11 jump diagnosis methods (the Lee-Heghinan method, rank test, slide F test, R/S method, Mann-Kendall method, Bayesian method, etc.). The diagnosis results were also classified into two types, trend results and jump results, and the results were synthesized. Finally, efficiency coefficients were calculated for identifying the alteration form, and the alteration form was judged to determine if one coefficient was bigger than the other one.

1.
In the range of the variability approach (RVA) method, ecological flow calculations should consider the hydrological regime. Parameters of the indicators of the hydrologic alteration (IHA) method are all closely related to runoff, and monthly flow influences river-living organisms, soil, etc. Therefore, the RVA analyses the monthly flow frequency distribution of each month and selects the flow corresponding to 25% and 75% in the frequency distribution as the upper and lower limits of the monthly flow, respectively [26]. Ecological flow is calculated as: where Q ejRVA is the ecological flow of the jth month, Q j is monthly mean flow of the jth month, Q j25% is the upper limit of RVA, and Q j75% is the lower limit of RVA.

2.
The steps for ecological flow calculation using the monthly frequency method (MFM) are as follows: (1) Calculate the monthly flow distribution empirical frequency of the jth month. (2) Select the probability distribution function and draw the monthly flow distribution theoretical frequency curve of the jth month. The Pearson III distribution curve and generalized extreme value (GEV) distribution curve are commonly used. The GEV distribution curve is better fitting with the runoff data [27].
Select the flow corresponding to the maximum frequency in the GEV distribution curve as the ecological flow of the jth month Q ejMFM [28].

3.
The annual daily mean flow frequency method (ADMFFM) assumes that daily flow may appear in a month with a certain probability. The frequency of daily flow in the jth month of the historical runoff series is calculated and the daily flow corresponding to 60% is selected as the ecological flow [24] of the jth month, Q ejADMFFM . The ecological flow calculated using this method meets the ecological water requirements in different cases.

4.
In the comprehensive ecological flow calculation method, ecological flow maintains riverine ecological function. Too much or too little will affect the riverine ecosystem. Therefore, the weights of each ecological flow calculation method are used to comprehensively calculate ecological flow. The steps for this method are as follows: (1) Evaluate ecological flows from RVA, MFM, and ADMFFM, and the comprehensive evaluation values of RVA, MFM, and ADMFFM of the jth month are W jRVA , W jMFM and W jADMFFM , respectively.
(2) Calculate the comprehensive ecological flow of the jth month as follows:

Weight Calculation of Comprehensive Ecological Flow
The three evaluation indexes listed in Section 2.2.4 were used to evaluate the calculation results following the above ecological flow calculation methods. Three single evaluation results were synthesized into a comprehensive evaluation result. The weights of each single ecological flow calculation result were obtained. The detailed calculation steps are as follows.

1.
We first determined the deviation rate of monthly flow. The median of the natural daily mean flow series in the jth month of N years was called the standard value. This index is a ratio of ecological flow in the jth month to the standard value. It reflects the deviation degree between the ecological flow and the natural flow in the same period [28,29]. It is computed as follows: (1) Calculate the natural daily mean flow in the jth month of N years Q ij , where i is the year, i = 1, 2, . . . , N; N is equal to the length of the runoff series; and j is the month, j = 1, 2, . . . , 12; (2) Sort Q ij in the ascending order Q 1j , Q 2j , . . . , Q nj , (Q 1j < Q 2j < . . . < Q nj , n is the serial order), where the standard value of the jth month is Q mid,j ; (3) The calculation formula of the deviation rate of the monthly flow is: where C j is the deviation rate of the monthly flow of the jth month and Q ej is the ecological flow of the jth month, in m 3 /s. If the value is close to 1, the calculated flow is approaching the natural flow and the riverine ecosystem is healthy.

2.
Satisfaction degree of the monthly ecological flow. This a ratio of the days where the natural flow is equal to or greater than the ecological flow to the total number of days in the same month. The formula to calculate the satisfaction of the monthly ecological water requirement is: where P j is the satisfaction degree of the monthly ecological flow of the jth month, D ej is the number of days where the natural flow was equal to or greater than the ecological flow in the jth month, D j is the total number of days in the jth month in N years, K num,j is the number of days in the jth month, Q ijk is the natural daily flow on the kth day of the jth month in the ith year in m 3 /s, and Q ej is the ecological flow of the jth month in m 3 /s. The greater the satisfaction degree of the monthly ecological flow, the healthier the riverine ecosystem.

3.
Suitability degree of the monthly ecological flow. The monthly ecological flow discrete coefficient is the sum of the discrete degree between the median and the characteristic extreme value of the ecological flow and natural flow. It reflects the suitability degree between the ecological flow and natural flow in a month.
where F j is the discrete coefficient of the monthly ecological flow in the jth month, Q emid,j is the median of the jth month's ecological flow of over N years, Q mid,j is the median of the jth month's natural flow of N years, Q emax,j is the maximum of the jth month's ecological flow of N years, and Q max,j is the maximum of the jth month's natural flow of N years. In the case of F j being greater than 10, it is considered as completely discrete and the value of F j is taken as 10. The suitability degree of the monthly ecological flow is computed as: When F j equals 1, a complete suitableness exists between the ecological flow and natural flow in the jth month; when it equals 0, absolutely no suitableness exists between the ecological flow and natural flow in the jth month. This index reflects the degree of suitability between the total ecological flow requirement and the natural flow in a month. When F j is 1, the total ecological flow requirement and the natural flow completely match. When F j is 0, the flow is completely unsuitable. A discrete coefficient greater than 10 can be considered completely discrete and the value of the discrete coefficient is equal to 10.

4.
Ecological flow comprehensive evaluation. The above indexes compare the ecological flow and natural flow from different aspects. The deviation rate of the monthly flow evaluates the degree and magnitude of deviation between the ecological flow and natural flow. The satisfaction degree of the monthly ecological flow temporally analyses the degree of satisfaction between the ecological flow and natural flow. The suitability degree of the monthly ecological flow evaluates the degree of suitability between the ecological flow and natural flow in discrete degrees. The above three indexes were synthesized into a comprehensive index α j to evaluate ecological flow in the jth month. The formula for calculating α j is: where α j ∈ [0, 1] and a value of α j close to 1 means the calculated ecological flow is approaching the natural flow and the riverine ecosystem is healthy.

5.
Weight calculation of the comprehensive ecological flow. The geometric mean method was used to calculate the weights for each ecological flow calculation method result.

Hydrological Alteration Diagnosis System
The results of the change in trend in the GI series obtained using the moving average method are shown in Figure 4; the moving average curve for the Hongqiling and Caizuizi stations showed a rising trend. In addition, slight fluctuations were observed in the moving average curves for the Baoqing and Baoan stations. Also, a downward trend exists in the mean annual flow curves shown in Figure 5. Additionally, GI is the index used to quantify and evaluate the evenness degrees of mean annual flow distribution, so there were different trends in GI curves and mean annual flow curves for the same station. Water 2018, 10, x 9 of 17 Based on the alteration analysis above, runoff data from 1966 to 2012 were selected as the subsequence to calculate ecological flow at the Baoqing station. Runoff data for the periods 1956-1990, 1956-1987, and 1956-1990 were selected as the subsequence to calculate ecological flow at the Baoan, Hongqiling, and Caizuizi stations, respectively. Results of the hydrological alteration diagnosis system are shown in Table 1.   Based on the alteration analysis above, runoff data from 1966 to 2012 were selected as the subsequence to calculate ecological flow at the Baoqing station. Runoff data for the periods 1956-1990, 1956-1987, and 1956-1990 were selected as the subsequence to calculate ecological flow at the Baoan, Hongqiling, and Caizuizi stations, respectively. Results of the hydrological alteration diagnosis system are shown in Table 1.  1956-1990, 1956-1987, and 1956-1990 were selected as the subsequence to calculate ecological flow at the Baoan, Hongqiling, and Caizuizi stations, respectively. Results of the hydrological alteration diagnosis system are shown in Table 1. Note: + in the table represents a strong alteration, and − represents a weak alteration. The runoff series was altered when the Hurst coefficient was not equal to 0.5 [25].

Results of Ecological Flow Calculation
Based on the selected subsequence above, the RVA method, monthly frequency method, and annual daily mean flow frequency method were used to calculate the ecological flows of the four hydrological stations. Comprehensive evaluation results for each method were determined by comprehensive evaluation in Section 2.2.4. The weights of the ecological flow calculated using each method were calculated from comprehensive evaluation results, and the comprehensive ecological flow based on the above methods was obtained. The ecological flows determined using the above methods and their weights are exhibited in Tables 2 and 3, respectively. The comprehensive ecological flow calculation results are shown in Figure 6.  (1) is the ecological flow calculated by subsequence series selected by the hydrological alteration diagnosis system. (2) is the ecological flow calculated by the whole series. RVA is the range of variability approach method, MFM is monthly frequency method and ADMFFM is annual daily mean flow frequency method. As shown in Figure 6, comparing the ecological flow values determined by different calculation methods, the ecological flow values from the RVA method were generally low and the ecological flow values from the monthly frequency method were relatively high. However, the comprehensive ecological flow values were all within the extreme values of the three calculation results. The comprehensive ecological flow and the ecological flow obtained from the three methods were all less than the monthly average flow, which is consistent with the characteristics of the ecological flow. In addition, the ecological flow process exactly coincides with the natural flow process, which is shown by the identical timing of extreme ecological flow and extreme natural flow.
The monthly flow meets the ecological water requirements when the monthly flow is greater than the ecological flow, otherwise the ecological water requirement is not satisfied, which means the riverine ecosystem may be damaged. The ecological water requirement damage frequency is a ratio of the number of times the monthly flow was less than the ecological flow to the series length. The As shown in Figure 6, comparing the ecological flow values determined by different calculation methods, the ecological flow values from the RVA method were generally low and the ecological flow values from the monthly frequency method were relatively high. However, the comprehensive ecological flow values were all within the extreme values of the three calculation results. The comprehensive ecological flow and the ecological flow obtained from the three methods were all less than the monthly average flow, which is consistent with the characteristics of the ecological flow. In addition, the ecological flow process exactly coincides with the natural flow process, which is shown by the identical timing of extreme ecological flow and extreme natural flow.
The monthly flow meets the ecological water requirements when the monthly flow is greater than the ecological flow, otherwise the ecological water requirement is not satisfied, which means the riverine ecosystem may be damaged. The ecological water requirement damage frequency is a ratio of the number of times the monthly flow was less than the ecological flow to the series length. The calculation results of the ecological water requirement damage frequency of the subsequences and the whole series are shown in Table 4.

Discussion
The ecological flow determined using the hydrological method, in which there is no consideration of the effect of some special situations such as overbank flow on the riverine ecosystem, was calculated from historical runoff data. Therefore, the calculated ecological flow regime is in concordance with the natural runoff regime. As the Naolihe basin is located in the cold region of China, precipitation from December to March in the next year is mostly dominated by sleet and solid precipitation involving the formation of snow. The sleet and solid precipitation flow into rivers in the form of snowmelt runoff in April and May in the following year. Flow in April and May accounts for a large proportion of annual flow. Therefore, rivers in the Naolihe basin generally have two flood processes: the spring freshet in April and May and the summer freshet in July and August. As shown in Figure 5, the ecological flow process, which corresponds to the natural flow process, also has two peaks. The maximum peak appears in May, which is consistent with organism reproduction needs in the alpine region.
In Figure 6, the comprehensive ecological flow values were between the values calculated by the three methods. Compared to the bigger values of the three methods, it is can be considered that more water in the river could be utilized for irrigation, drinking, power generation, and so on based on values of the comprehensive ecological flow. At the same time, the water utilization, according to the comprehensive ecological flow, is enough to guarantee the water requirement for the river ecosystem, because they exceed the smaller values calculated by the three methods.
In the Baoqing and Baoan stations, located in the upstream region of the Naolihe river, human activities have relatively little impact on the riverine ecosystem. Therefore, differences in the ecological water requirement damage frequency for the subsequence series and the whole series were not significant. However, the ecological water requirement damage frequency increased obviously in July, August, and September. This is because the study area is located in an agricultural area, and crops need more water in this period. Human use of water resources for irrigation tends to enhance the ecological water requirement damage frequency. However, ecological water requirement damage frequency averages for the subsequence series and the whole series of the Hongqiling station and Caizuizi station were basically the same. This is because the Hongqiling station and Caizuizi station are located in the middle and downstream portions of the Naolihe river, respectively, and the construction of the Longtouqiao and Hamatong Reservoirs effectively solves the problem of instream flow decreasing during the non-flood season. Similarly, during the flood season, due to the impoundment process of reservoirs, part of the natural flow is stored. The instream flow decreases and the ecological water requirement damage frequency increases. In addition, the average ecological water requirement damage frequencies in the Hongqiling station and Caizuizi station were more than in the Baoqing station and Baoan station. This finding is also because of human development. First, sort the monthly flows into a group in ascending order; the result is listed in column (3) in Table A2.
Second, accumulate the ascending monthly flow series by Equation (1); the result is listed in column (4).
Fourth, draw a Lorenz curve with the series in column (5) and (6) in Table A2.  (2), and the result was 0.5752.
The calculation results for each year at the four stations are shown in Table A3.  Then, calculate the areas of A i and B i ; the area of A i is 0.2876 and the sum of the areas of A i and B i is 0.5. So, the GI value was calculated by Equation (2), and the result was 0.5752.
The calculation results for each year at the four stations are shown in Table A3.