Peak flows are the most important flood parameter which relatively reflects the highest level and potential destructive power of a flood. Understanding peak flow changes can effectively capture a flood characteristic and is essential for developing flood control strategies. This study aims to reveal how regional peak flows evolved in recent decades, mainly from a non-linear perspective. The Beijiang River Basin (BRB) was chosen for the analysis, and hydrological data from four hydrologic stations were used. Methods including ensemble empirical mode decomposition and rescaled range analysis were applied to advance the research. Results indicate a non-significant uptrend and a multiple periodicity of peak flows in BRB. However, short periods were more distinct than long ones. In the future, peak flows may continue to increase over time. Such changes in peak flows are possibly due to local reservoir operations and the changing South Asian Summer Monsoon (SASM). The research suggests an increasing flood risk and recommends more regional flood adaptations to avoid flood losses for BRB. Synchronously, it provides a reference for studies regarding periodicity and the future trend of peak flows in other regions.

## INTRODUCTION

Flood hazards are well known to have a strong impact on society and the environment and, thus, are currently drawing the attention of the hydrological research community. Previous investigations of flood events usually included studies on their mechanisms, impacts and characteristics (Leviandier *et al.* 2000; Marchi *et al.* 2010; Haddad *et al.* 2012; Munoz *et al.* 2012). Among these studies, the analysis of flood parameters, such as duration, frequency and flow, has enabled a more direct understanding of flood mechanisms. As such, it is crucial to analyze variations of these flood parameters in terms of flood assessment and flood forecasts.

Peak flow, the maximum flow during a flood process, is one of the most important flood parameters. On the one hand, peak flows often reveal the aspects of the intensity and even the total flow of a flood. By extension, it corresponds to the highest flood level and relatively reflects the potential destructive power of a flood. On the other hand, estimation of the design floods of hydraulic structures, such as dikes, spillways and storm water evacuation canals, requires the determination of peak flows based on the stream flow fluctuations (Fill & Steiner 2003; Munoz *et al.* 2012). Consequently, flood control and disaster mitigation planning usually depend a great deal on peak flow examination.

It is an accepted fact that abnormal climate and human activities have triggered extreme precipitation, causing flood events (IPCC 2012). Flood peak flows are expected to become more variable and unpredictable, due to the impacts of global change on climate, storm-weather systems and river discharge conditions. Recent studies have attempted to discuss how flood peak flows will respond to changes in precipitation and land use. They have proven the close relationship between flood peak flows and precipitation (Gottschalk & Weingartner 1998; Calenda *et al.* 2005; Daniels 2007). Furthermore, these studies support the view that human activities, such as urbanization, deforestation and afforestation, affect flood peak flows (Beschta *et al.* 2000; Guillemette *et al.* 2005; El Alfy 2016; Mei *et al.* 2016). In addition to these issues, several researchers have focused on the variability of flood peak flows. For instance, previous studies have revealed that peak flow variability would decrease when average flows are expected (Pearce *et al.* 1980; Hewlett 1982). The variation coefficient of annual flood peak flows tends to decrease as the catchment area increases (Kuzuha *et al.* 2009). An example of applying a distribution model to flood peak flows is seen in the study of the River Tiber in Rome; the Gumbel distribution best fits the high return period quantiles of peak flows in the River Tiber Basin (Calenda *et al.* 2009).

Many of the studies referenced above focus on peak flow trends and influential factors, rather than periodicities and fractal properties of annual peak flows. A thorough understanding of the periodicity of peak flows can provide a reference for regional flood prediction and help to further explore the influential factors, whilst studying the fractal properties of peak flows can provide a reference for relevant researches of peak flows from a non-linear perspective. In addition, it is particularly necessary for small regions with prosperous economies and large populations, where increased flood peak flows are expected to threaten local residents (Kuzuha *et al.* 2009). However, up to now few studies have focused on changes in peak flows in the Beijiang River Basin (BRB), a prosperous area in South China where major floods often occur in the rainy season (Zhang *et al.* 2007; Wang *et al.* 2012, 2015).

Taking BRB as a study case, the main objectives of this study are to: (1) detect whether basin-scale peak flows showed significant trends in recent decades; (2) reveal how basin-scale peak flows evolved from a non-linear perspective; (3) discuss whether climate change and human activity have impacts on regional peak flow changes and, if so, how. The study would be useful for local regulators to manage and can provide a reference for relevant researches on periodicity and future trends of peak flows across other regions.

The paper is organized as follows. The study area and data sources are described first; methodologies such as Ensemble Empirical Mode Decomposition (EEMD) and R/S are presented next; and results are listed afterwards. Finally, the discussion and conclusion sections are presented.

## STUDY AREA AND DATA SOURCES

### Study area

*et al.*2011). It has 13 tributaries and its watershed area is more than 1,000 km

^{2}. The annual mean temperature in the BRB is approximately 21 °C. Precipitation in the BRB is usually concentrated in the warm season from April to September, which accounts for nearly 70% of the annual total (Luo

*et al.*2008). During spring and summer, warm air masses from the sea are blocked by the Nanling Mountains, causing frequent rainstorms and leading to subsequent floods in the BRB. The extraordinary historic flood that occurred in June 1994 had a peak flow of up to 19,700 m

^{3}/s at Shijiao station and caused serious damage to the basin, including the loss of 371 lives (see www.gdbjdd.com. cn/). The recent historic flood that occurred in July 2006 also brought enormous damage to the local society.

### Data sources

The Hydrology Bureau of Guangdong Province (HBGP) collected the flow data used in the study from four hydrological stations located in the main stream of Beijiang River (Figure 1 and Table 1). The quality of the flow data is strictly controlled by HBGP, which is officially responsible for the regional hydrological data. For each station, the annual peak flow series was obtained using the annual maximum value method (Mitosek *et al.* 2006). Therefore, there is only one maximum instant peak flow for each year for each station.

Station | Fenshi | Pingshi | Lishi | Shijiao |
---|---|---|---|---|

Measurement period | 1970–2007 | 1964–2007 | 1955–2007 | 1954–2007 |

Control area (km^{2}) | 880 | 3,567 | 6,976 | 38,363 |

Station | Fenshi | Pingshi | Lishi | Shijiao |
---|---|---|---|---|

Measurement period | 1970–2007 | 1964–2007 | 1955–2007 | 1954–2007 |

Control area (km^{2}) | 880 | 3,567 | 6,976 | 38,363 |

For investigating the possible climatic causes of changes in the annual peak flows in BRB, this study considers five climatic indices: the Multivariate El Niño-Southern Oscillation index (MEI), the East Asian Summer Monsoon Index (EASMI), the South Asian Summer Monsoon Index (SASMI), the South China Sea Summer Monsoon Index (SCSSMI), and the Dipole Mode Index of India Ocean (DMI). EASMI/SASMI/SCSSMI is defined as an area-averaged seasonally dynamical normalised seasonality at 850/850/925 hPa within the East Asian monsoon domain (10–40 °N, 110–140 °E)/South Asian domain (5–22.5 °N, 35–97.5 °E)/South China Sea monsoon domain (0–25 °N, 100–125 °E) (Li & Zeng 2005). These monsoon indices were developed by Li & Zeng (2005). The MEI data set was obtained from the US National Oceanic and Atmospheric Administration Earth System Research Laboratory (www.cdc.noaa.gov/people/klaus.wolter/MEI/). The EASMI, SASMI and SCSSMI data sets were obtained from the United States NCEP/NCAR monthly mean data sets (http://ljp.lasg.ac.cn/dct/page/65544), whilst the DMI data set was obtained from the Japan Marine Science and Technology Center (www.jamstec.go.jp/frcgc/research/d1/iod/HTML/Dipole%20Mode%20Index.html).

## METHODOLOGY

### Linear regression

*et al.*2013). The linear regression equation is described as: where

*X*is the independent variable,

*Y*is the dependent variable,

*a*is the regression coefficient that reflects the rate of changes of the time series, and

*b*is the regression constant. In this study, the parametric

*t*-test was employed and the

*P*values for the

*t*-test were computed to identify whether the trends in the peak flow series are statistically significant (Yue & Pilon 2004).

### Multiple linear regression

*Y*is the response variable and

*X*,

_{1}*X*,…,

_{2}*X*are the predictors, and , , … , are the coefficients corresponding to each of the predictors. is a constant. For a predictor,

_{n}*X*, the contribution of

_{i}*X*to

_{i}*Y*is determined by the following equation:

### Mann–Kendall analysis

The Mann–Kendall (M-K) method is widely used for identifying abrupt changes in a time series. The advantage of M-K is that the series does not need to comply with a certain distribution of samples, avoiding interference from abnormal values. A detailed algorithm has been described elsewhere (Mosmann *et al.* 2004; Tabari & Hosseinzadeh Talaee 2013; Wu *et al.* 2015).

Then define another statistic *UB _{k}*. Similarly, the values of

*UB*are computed backward, starting from the end of the time series.

_{k}By drawing the *UF _{k}* and

*UB*curves on the same graph, the point of abrupt change in the series can be identified. If the

_{k}*UF*and

_{k}*UB*curves intersect, and then diverge, acquiring a specific threshold value of 1.96 at a 0.05 significance level, there is a statistically significant trend. The point where they intersect shows the approximate year at which the abrupt change occurred.

_{k}### Ensemble empirical mode decomposition

This specific method is derived from the method, Empirical Mode Decomposition (EMD). Before looking at the EEMD, a brief description of the original EMD is proposed. Presented in 1998, EMD is a method for the time series analysis of non-stationary and non-linear signals (Huang *et al.* 1998). It decomposes non-stationary data into a collection of intrinsic mode function (IMF) components with a residue component (Res.) using different time scales (Bi *et al.* 2010). The original series is decomposed into different time scale fluctuations, from which the periodicities of the original series can be detected (Ding *et al.* 2010). The procedure for extracting IMFs from a signal is as follows (Kaluzynski 2014; Kim *et al.* 2014):

Identify all local extrema over the entire time period of the signal,

*x*(*t*).Create an envelope of the local maxima,

*e*(_{up}*t*), and minima,*e*(_{low}*t*).Check whether

*d*_{1}(*t*) satisfies the following criteria of an IMF:the number of zeros and extrema are equal or differ by no more than 1;

the sum of the envelopes of maxima and of minima of an IMF is zero.

*d*

_{1}(

*t*) does not satisfy the criteria of an IMF, set

*d*

_{1}(

*t*) to

*x*(

*t*), then repeat steps (1)–(5) until

*d*

_{1}(

*t*) satisfies the criteria. The first signal to satisfy the criteria is termed as the first IMF,

*h*

_{1}(

*t*).

Reiterate steps (1)–(6) until all IMFs are identified. The process is completed when the residual is constant, monotonic or has a single extremum.

Suggested by Wu & Huang (2009), mode mixing can be the significant drawback of EMD (Wu & Huang 2009; Wang *et al.* 2012). To overcome the shortcoming of EMD, the EEMD method is proposed. It is a new noise assisted analysis method to obtain the actual time–frequency distribution of the seismic signal. With the procedure of extracting IMFs by EMD in mind, the EEMD is developed as follows (Wu & Huang 2009):

Add a white noise series to the targeted data.

Decompose the data with added white noise series into IMFs.

Repeat steps (1) and (2), but with different white noise series each time.

Calculate the means of corresponding IMFs of the decompositions as the final results.

The added white noise series would cancel each other in the final mean of the corresponding IMFs. The mean IMFs stay within the natural dyadic filter windows, and therefore reduce significantly the chance of mode mixing and preserve the dyadic property. In this study, the ratio between the standard deviation of the white noise series and that of the targeted data is 0.2, whilst the repeat times of steps (1) and (2) is 20.

### Rescaled range analysis

The R/S is used to estimate the fractal properties of a time series. The main idea of the R/S is that one looks at the scaling behavior of the rescaled cumulative derivations either from the mean or from the distance the system travels as a function of time (Karakasidis & Liakopoulos 2004).

Equation (7) presents a fitted straight line in a log-log plot of R/S as a function of *τ*, and *H* is the slope of the straight line. In the case of only short-range correlations (uncorrelated), the walk proﬁle displays the properties of a standard random walk with *H* = 0.5. If fluctuations in subsequent values are positively correlated, then *H* > 0.5; if fluctuations in subsequent values are negatively correlated, *H* < 0.5.

## RESULTS

### Trends

^{3}/s·a, respectively. The maximum, minimum and mean values of peak flows in Table 2 indicate that the annual peak flows tended to increase from the Fenshi to the Shijiao station. The

*P*values of the series are respectively 0.375, 0.098, 0.105 and 0.324 (Table 2), all of which are greater than 0.05. In other words, none of the upward trends are statistically significant. The maximum eigenvalues of the time series at the Fenshi, Pingshi, Lishi and Shijiao stations are 1,990, 4,830, 8,800 and 17,400 m

^{3}/s, respectively. With the exception of the maximum at the Fenshi station that occurred in 1994, all the stations' maxima occurred in 2006. Furthermore, the maxima of the series increased from the Fenshi to the Shijiao station. The coefficients of deviation (

*C*) of the four series, on the other hand, decrease in value from the Fenshi to the Shijiao station. This signifies that the annual peak flow series of the Fenshi station was more discrete, though with a lower magnitude, compared with the other three stations.

_{v}Station | Fenshi | Pingshi | Lishi | Shijiao |
---|---|---|---|---|

Maximum (m^{3}/s) | 1,990 | 4,830 | 8,800 | 17,400 |

Minimum (m^{3}/s) | 184 | 636 | 630 | 2,740 |

Mean (m^{3}/s) | 629 | 1,591 | 2,538 | 9,876 |

P value | 0.375 | 0.098 | 0.105 | 0.324 |

Coefficient of deviation (C) _{v} | 0.569 | 0.495 | 0.521 | 0.319 |

Coefficient of skew (C) _{s} | 1.810 | 1.992 | 2.402 | 0.234 |

Station | Fenshi | Pingshi | Lishi | Shijiao |
---|---|---|---|---|

Maximum (m^{3}/s) | 1,990 | 4,830 | 8,800 | 17,400 |

Minimum (m^{3}/s) | 184 | 636 | 630 | 2,740 |

Mean (m^{3}/s) | 629 | 1,591 | 2,538 | 9,876 |

P value | 0.375 | 0.098 | 0.105 | 0.324 |

Coefficient of deviation (C) _{v} | 0.569 | 0.495 | 0.521 | 0.319 |

Coefficient of skew (C) _{s} | 1.810 | 1.992 | 2.402 | 0.234 |

*P* ≤ 0.05 denotes the trend in the series is beyond the 0.05 significance level.

### Abrupt changes

*UF*and

_{k}*UF*curves for the four series intersect at more than one time point, though none of the

_{B}*UF*curves exceed the critical value of ±1.96. Therefore, none of the abrupt change points meet the 0.05 significance level, i.e. the abrupt changes of the annual peak flows in BRB were non-significant for their own measurement periods.

_{k}### Stationary characteristics and periodicities

Although the trends and abrupt changes in peak flows at the four stations were non-significant, the stationarity characteristics should be further examined. We used the Augmented Dickey–Fuller (ADF) unit root test to analyze the stationarities of peak flow series for the four stations. The ADF method is a classic method which can examine the stationarity of a time series (Zhao & Chen 2015). Table 3 lists the *P* value statistics from ADF. Note that the *P* value statistics in the Fenshi, Pingshi and Lishi series are greater than the critical value 0.05, implying the nonstationary nature of peak flows at these stations. On the contrary, the *P* value statistic in the Shijiao series is smaller than 0.05, indicating the stationarity of the series. As the Feilaixia reservoir operation may alter, changing the properties of peak flows in the Shijiao station, we omitted the records after 1998, when the Feilaixia reservoir was constructed, and reanalyzed the stationarity of the new series. Results show, as expected, a nonstationary characteristic of the peak flow series covering 1954–1997 at the Shijiao station (Table 3). These results on the one hand can explain the nonstationarity of peak flows in BRB, and on the other hand infer a possible influence of the Feilaixia reservoir operation on the changing property of peak flows at the Shijiao station.

Station | Dickey–Fuller | Lag order | P value |
---|---|---|---|

Fenshi | −3.03 | 3 | 0.17 |

Pingshi | −3.08 | 3 | 0.15 |

Lishi | −3.15 | 3 | 0.11 |

Shijiao (during 1954–2007) | −3.74 | 3 | 0.03 |

Shijiao (during 1954–1997) | −3.46 | 3 | 0.06 |

Station | Dickey–Fuller | Lag order | P value |
---|---|---|---|

Fenshi | −3.03 | 3 | 0.17 |

Pingshi | −3.08 | 3 | 0.15 |

Lishi | −3.15 | 3 | 0.11 |

Shijiao (during 1954–2007) | −3.74 | 3 | 0.03 |

Shijiao (during 1954–1997) | −3.46 | 3 | 0.06 |

*P* ≤ 0.05 and *P**>* 0.05 indicate the series are stationary and nonstationary, respectively.

*et al.*2010). Therefore, the IMFs reflect different periods of the annual peak flow series, i.e. its short, medium and long periods. Taking the Shijiao station as an example, we analyzed the periodicity of peak flow series based on the IMFs, Res. and power spectra.

As shown in Figure 4(d), the Shijiao series showed oscillations with complex frequencies over the entire time interval in the IMF1; however, we saw a 2.1-year period and a 4-year period at the 0.05 significance level in the IMF1 based on the corresponding power spectrum (Figure 5). In the IMF2, the series fluctuates chaotically at first glance, however, a significant 5.9-year period is detected for IMF2 from the power spectrum. For the IMF3, its power spectrum reveals multiple medium and long periodicities of the series, in which >8-year periods are significant with high spectrum power beyond the 0.05 significance level. The IMF4 remains almost unchanged before 1990, after which it increases gradually. For the IMF5, it decreases before 1960, and during 1970–1990, whereas it increases from 1960 to 1970 and after 1990. By computing the power spectra of IMF4 and 5, however, we could not find any significant periods. In conclusion, the overall IMFs show short, medium and long periodicities of the annual peak flow series at the Shijiao station.

The periodicities of the other three series of the remaining stations are listed in Table 4 based on their own IMFs and Res., and their power spectra. Concisely, the series from the Fenshi station had 3.2, 4.2 and 7-year periodicities, while that of the Pingshi station had 2.2 and 4.2-year periodicities. In addition, the Lishi station series had 2.2, 4 and >6-year periodicities. Specifically, we found that, by comparison, the peak flow series at the Pingshi station was mainly dominated by short periods, while that at the Fenshi station was dominated by short and medium periods. For the Fenshi and Pingshi stations, not only short and medium periods were remarkable in the peak flow series, but also long periods.

Station | IMF1 | IMF2 | IMF3 | IMF4 | IMF5 |
---|---|---|---|---|---|

Fenshi | 3.2 | 4.2 | 7 | / | / |

Pingshi | 2.2 | 4.2 | / | / | / |

Lishi | 2.2 and 4 | 6 | >6 | / | / |

Shijiao | 2.1 and 4 | 5.9 | >8 | / | / |

Station | IMF1 | IMF2 | IMF3 | IMF4 | IMF5 |
---|---|---|---|---|---|

Fenshi | 3.2 | 4.2 | 7 | / | / |

Pingshi | 2.2 | 4.2 | / | / | / |

Lishi | 2.2 and 4 | 6 | >6 | / | / |

Shijiao | 2.1 and 4 | 5.9 | >8 | / | / |

‘/’ denotes no significant periodicity is detected.

### Fractal property analysis

*H*value of the Shijiao station data is 0.503, which is very close to 0.5, demonstrating a nearly pure random walk characteristic of the series. Thus, it is unsure whether the future trend of the series will follow that of the past or not. On the other hand, by comparing the

*H*values of the four stations, it can be concluded that correlations between past and future trends gradually weaken from the Fenshi to the Shijiao station.

Station | Fenshi | Pingshi | Lishi | Shijiao |
---|---|---|---|---|

Hurst exponent | 0.687 | 0.640 | 0.582 | 0.503 |

Station | Fenshi | Pingshi | Lishi | Shijiao |
---|---|---|---|---|

Hurst exponent | 0.687 | 0.640 | 0.582 | 0.503 |

## DISCUSSION

### Reviews and implications

The characterization of the annual peak flows across BRB reveals the temporal patterns of its variability, which can not only enable comparison between such changes and those reported in other regions, but can also be useful for local farmers and regulators (Pisaniello 2016). Previous research has revealed an increase in annual peak stream flow in BRB (Zhang *et al.* 2015), and our study supports this point of view. Yang *et al.* (2012) implied some shifts together with an approximate six-month periodicity in peak flows during 1957–1967 in the Hailiutu River, Northwest China, whereas our study suggests no abrupt changes but multiple periodicities of peak flows in BRB. This may be ascribed to different precipitation regimes between these two regions (Wu *et al.* 2015). There is an increased necessity for flood control strategies in BRB during the past decades, as well as in the future, in light of the detected upward trends and long-term correlation characteristic in the recorded data. The unique multi-scale periodic oscillation characteristics of the peak flows, on the other hand, may provide a reference for flood prediction.

### Possible climatic and anthropogenic causes

In the current study, distinct differences between the variabilities of annual peak flows at the four stations are observed. Specifically, from the Fenshi to the Shijiao station, both the coefficients of deviation (Table 2) and long-term correlation (Table 5) of peak flows decrease in value. Such geographical contrasts may, in part, be due to the geomorphology (Whitaker *et al.* 2002). Geomorphologic factors are known to affect flood occurrences through two main mechanisms: orographic effects that augment precipitation, and topographic reliefs that promote the rapid concentration of flow (Costa 1987). BRB is considered to be a small watershed but with topographic contrast (Figure 1). Topographic reliefs in this watershed diminish from the upstream to the downstream, causing peak flows to vary spatially. In addition to these physiographic factors, human activities are assumed to be responsible for such phenomena. The peak flow series at the Shijiao station, for instance, is stationary for the period 1954–2007, but is non-stationary for the period 1954–1998, i.e. before the Feilaixia reservoir was constructed. Zhang *et al.* (2015) have confirmed the impacts of reservoir regulations on peak flows changes in the Pearl River Basin (with BRB involved), and this finding supports our point of view of Feilaixia reservoir affecting peak flows at the Shijiao station.

Climate change may influence flood occurrences through its effects on the variability of storm characteristics and the seasonality of rainfall and evapotranspiration, which affect the antecedent watershed conditions for individual storm events (Marchi *et al.* 2010). Assessments of the influences of climatic factors on hydrologic variables, such as precipitation and flow, are usually undertaken to interpret changes in these variables. To explore possible climatic causes of the annual peak flow changes in BRB, five climatic indices were selected after considering global and regional climate characteristics, i.e. EASMI, SASMI, SCSSMI, MEI and DMI. We used MLR to investigate the relationships between annual peak flows in the Fenshi, Lishi, Pingshi stations and these indices, and the contributions of these indices to peak flows. We did not consider the Shijiao station series because this station should be influenced by the Feilaixia reservoir, where the effect of the subjective operation may be greater than the climate change impact. The results indicate that the five indices, except the SCSSMI, may positively influence peak flows (Table 6). Of all the indices, SASMI contributed the most to the peak flow changes; a strengthened/weakened South Asian Summer Monsoon (SASM) may increase/decrease peak flows at the Fenshi, Lishi and Pingshi Stations in BRB.

Station | Indices | EASMI | SASMI | SCSSMI | MEI | DMI |
---|---|---|---|---|---|---|

Fenshi | Coefficient | 0.024 | 0.334 | 0.152 | 0.147 | 0.089 |

Contribution | 0.033 | 0.447 | 0.204 | 0.197 | 0.119 | |

Pingshi | Coefficient | 0.203 | 0.347 | −0.101 | 0.198 | 0.021 |

Contribution | 0.234 | 0.399 | 0.116 | 0.227 | 0.024 | |

Lishi | Coefficient | 0.129 | 0.283 | −0.070 | 0.131 | 0.075 |

Contribution | 0.188 | 0.411 | 0.102 | 0.190 | 0.109 |

Station | Indices | EASMI | SASMI | SCSSMI | MEI | DMI |
---|---|---|---|---|---|---|

Fenshi | Coefficient | 0.024 | 0.334 | 0.152 | 0.147 | 0.089 |

Contribution | 0.033 | 0.447 | 0.204 | 0.197 | 0.119 | |

Pingshi | Coefficient | 0.203 | 0.347 | −0.101 | 0.198 | 0.021 |

Contribution | 0.234 | 0.399 | 0.116 | 0.227 | 0.024 | |

Lishi | Coefficient | 0.129 | 0.283 | −0.070 | 0.131 | 0.075 |

Contribution | 0.188 | 0.411 | 0.102 | 0.190 | 0.109 |

Coefficient means MLR coefficient. Magnitude of contribution is calculated according to Equation (3).

The overall analyses reveal possible impacts of climatic and anthropogenic factors on annual peak flows in BRB. Nevertheless, the underlying mechanisms of the causes are far more complex, and a further study on this issue should be conducted to substantiate the MLR results.

## CONCLUSIONS

This study applied several methods including EEMD and R/S to examining the trends, abrupt changes, periodicities and fractal properties of annual peak flows in BRB, based on the hydrological data from four hydrological stations. In addition, the climatic and anthropogenic causes of peak flow changes were also explored. Below summarizes what can be inferred from this research:

For the past decades annual peak flows have increased but the uptrend was non-significant in the study region, which suggests a potential increase in regional flood risk. Meanwhile, no abrupt change occurred in peak flow series.

From a non-linear perspective, peak flows displayed a multi-periodic nature, in which both short and long periods were seen but the short ones were more distinct. Moreover, the fractal property analysis reveals that peak flows may continue to increase in the future, causing a higher flood risk than the past. Therefore, flood protection might need to be upgraded to avoid unacceptable losses.

Peak flow variability in BRB may be attributed to both local reservoir operations and the changing SASM. A strengthened/weakened SASM possibly causes an increase/decrease in peak flows, while a local reservoir might alter the stationary characteristics of peak flows.

These findings would be useful for flood prediction and defenses in the study region. More importantly, the research can provide a reference for other regions regarding periodicity and future trend of peak flows, as well as the impacts of human activities and changes in climate.

## ACKNOWLEDGEMENTS

The research is financially supported by the National Natural Science Foundation of China (Grant No. 51209095, 51579105, 51210013, 51479216), the National Science and Technology Support Program (Grant No. 2012BAC21B0103), the Water Science and Technology Innovation Project of Guangdong Province, the Fundamental Research Funds for the Central Universities (2014ZZ0027).