Daily maximum runoff frequency analysis under non-stationary conditions due to climate change in the future period: Case study Ghareh Sou basin

‘Floods are considered as most destructive among all natural hazards which impose lots of damages on human societies’ (Daliran Firouz et al. 2016). Therefore, predicting the flood time and its intensity rate in the future is one of the most important issues which can increase the accuracy of constructing the hydraulic structures and decrease the amount of damage. ‘The design of hydraulic structures and flood risk management is often based on instantaneous peak flow’ (Jimeno-Saez et al. 2017). One of the methods of estimating the maximum instantaneous flood is to analyze the statistics of the hydrometric station data in the basin, while the time series of the flow are very scarce, and they usually have limited time intervals. Therefore, calculating the maximum instantaneous flood using the experimental methods is an efficient solution.

To deal with the dangers of floods, it is necessary to know this phenomenon, the factors affecting it, and the appropriate methods for estimating it. In general, the calculation of floods in different return periods is provided by analyzing the hydrometric station data. However, in the absence of these stations or the incomplete or short-term data, the indirect methods are utilized to estimate floods, one of which is the use of experimental formulas (Dastorani et al. 2011). Some of these empirical methods are Fuller, Sangal, Fill Steiner, and one of the statistical methods is artificial intelligence. Using these methods, attempts have been made to connect the maximum instantaneous flood in an area with some physical characteristics of the basin.

For example, Dastorani et al. (2011) used four methods of Fuller, Fill Steiner, Sangal, and artificial intelligence on 12 stations in Iran and calculated the maximum instantaneous flood from the data of the average daily discharge. Their studies demonstrate that artificial intelligence has better performance than other methods in all stations. Shabani (2013) utilized empirical and artificial intelligence methods to calculate the maximum instantaneous flood in the Kharestan basin in Iran. In this research, comparing the four methods of Fuller, Sangal, Fill Steiner, and artificial intelligence indicated that the artificial intelligence method has better performance in contrast to other methods. Jimeno-Saez et al. (2017) used three methods of Fuller, Sangal, and Fill Steiner and also two methods of the artificial intelligence and Anfis (adaptive network-based fuzzy inference system) in the peninsular basin of Spain in order to estimate the instantaneous peak flow and compared them with each other. The results represented that the Fuller equation has fewer errors and more accurate performance compared with the other two equations.

‘The definition of a stationary time series is one where all finite dimensional distributions are invariant for changes in time’ (Madsen 2007). ‘Traditional approaches assume that hydrological processes evolve in an environment where the hydrological cycle is stationary over time. However, in recent years, it has become increasingly evident that in many areas of the world the foregoing assumption may no longer apply, due to the effect of anthropogenic and climatic induced stressors that cause non-stationary conditions’ (Salas et al. 2018). According to recent studies, the main reasons for non-stationarity can be climate change, anthropogenic effects on the hydrological cycle, and low-frequency climatic variability (Milly et al. 2008; Bayazit 2015; Cancelliere 2017; Salas et al. 2018). Because of these reasons, the hydrological regime of rivers is nonstationary. Thus, it is necessary to review the flood frequency analysis, which was previously provided based on the stationary assumption and estimate based on the assumption of non-stationarity. Different studies have been provided in the flood frequency analysis field, considering the non-stationary assumption. Rajagopalan (2010) considered two approaches for the non-stationary frequency analysis. The first approach is to use the local polynomial probability method. The second approach is to use the General Extreme Value (GEV) with the auxiliary variables, in which the scale parameter is variable over time, and depends on the auxiliary variables. Also, Cheng & Aghakouchak (2014) introduced a framework for the stationary and non-stationary estimations using the Bayesian inference. This framework was performed in Non-Stationary Extreme Value Analysis (NEVA) software. The results demonstrate that NEVA can confidently describe the extreme values and their return levels. Also, Condon et al. (2015), in a case study, analyzed the future flood frequency using the non-stationary extreme value models and design-life methodology. The historical flood was simulated in two locations utilizing the variable infiltration capacity (VIC) model and the GEV model. This paper presents an end-to-end analysis using the GEV methods along with the current downscaled climate predictions. Limitations can be omitted using the GEV model because the model parameters like the mean and the variance can be considered as a function of time or auxiliary variables.

The purpose of the present study is to compare the results of maximum instantaneous flood analysis under the stationary and non-stationary assumption; it also answers the question: ‘whether there is a difference between estimating the return period of maximum instantaneous flood in stationary and non-stationary conditions’.

The study area of this research is the Ghareh Sou sub-basin, which is located in the northwest of the Karkheh basin in Kermanshah province, in the west of Iran. Table 1 shows the summarized information of the study area. Figure 1 shows the location of the study area in Iran.

The Ghareh Sou River is one of the most important permanent rivers in the Karkheh basin. It has more water in winter and spring and decreases in summer. In a 32-year period, the average flow of Ghareh Sou was 289 million cubic meters per year. Figure 2 shows Ghareh Sou sub-basin, the Ghareh Sou River's position and the hydrometric station in the study area.

According to Figure 2, the study area was divided into 29 cells of 0.125 × 0.125° (13.75 × 13.75 km). The daily data of the minimum temperature, maximum temperature, rainfall and wind speed of the European Centre for Medium-Range Weather Forecasts (ECMWF) database at scales of 0.125 × 0.125° from 2010 to 2013 and the daily runoff data of the Doab Mark hydrometric station have been used to calibrate and validate the VIC model, which will be explained in the following.

‘Ghareh Sou basin has experienced severe droughts, and floods, in recent years’. (Ashrafi et al. 2020). Kermanshah province has about 550 villages that suffered from a water shortage crisis. The source of income of most villagers in this area depends on water resources (agriculture and animal husbandry). Due to the decrease in rainfall in recent years, it has caused economic problems for the people of the region (Shafiei et al. 2019). On the other hand, the occurrence of floods in the region has caused financial and human losses (the last destructive flood related to March 2019 and the effects of it can be seen to this day). Therefore, predicting the maximum instantaneous flood and the amount of rainfall reduction in the future period can help in better management and planning.

Figure 3 illustrates the flowchart of the stages of the research.

The present study includes a main purpose and four side purposes. The main purpose is investigation of the effect of non-stationarity on maximum instantaneous flood and the side purposes include:

• 1.

  evaluation of VIC model;

• 2.

  evaluation of methods for estimating the maximum instantaneous flood from maximum daily discharge;

• 3.

  trend analysis on maximum instantaneous flood data;

• 4.

  comparison of daily temperature, precipitation and runoff in the past and future periods.

In order to achieve the purposes mentioned above, the following steps were performed:

• (1)

  In this study, the VIC model was used to generate daily runoff. In order to calibrate the model, ECMWF daily data (maximum and minimum temperature, wind speed and rainfall) and daily runoff data (Doab Marak station) for the period 2010–2012 were used. For validation of the VIC model, ECMWF daily data and daily runoff data for the period 2012–2013 were also used.

• (2)

  Hadley Centre Coupled Model, version 3 (HadCM3) (under two scenarios of representative concentration pathway (RCP) 2.6 and 8.5) was used to generate daily runoff by the VIC model in the future period. This was downscaled by the Statistical DownScaling Model (SDSM). The SDSM model requires two series of data (called predictable and predictor). The daily data of ECMWF in the period of 1971–2001 were used as predictable and the daily data of NCEP in the same period of time were used as predictor.

• (3)

  Using HadCM3 downscaled data and VIC model, daily runoff for the future period (2020–2049) was generated.

• (4)

  Four methods of Fuller, Senegal, Fill Steiner and artificial neural network (ANN) were compared with each other to generate maximum instantaneous flood from daily discharge data, and finally the ANN method was selected to generate maximum instantaneous flood for the future period.

• (5)

  Trend analysis and autocorrelation analysis were performed on the data of maximum instantaneous flood in the past (1984–2013) and future periods.

• (6)

  Finally, using the NEVA software package, the amount of maximum instantaneous flood in the stationary and non-stationary conditions in the past and future period under different return periods were calculated.

‘Presently, the most reliable tool for producing climate scenarios is the coupled three-dimensional model of Atmosphere-Ocean General Circulation (AOGCM)’ (Mansouri et al. 2015).

In this research, the HadCM3 model that was described by Gordon et al. (2000) and Pope et al. (2000) under two emission scenarios of RCP2.6 and RCP8.5 has been used to predict the future time intervals (2020–2049) on the daily scale. The scenario RCP2.6 represents the most optimistic state (the descending trend of the greenhouse gas emissions until 2100), and the scenario RCP8.5 represents the most pessimistic state (the ascending trend with the highest slope of the greenhouse gas emissions until 2100). Since greenhouse gas emissions in the future are uncertain, in this study both the most optimistic and the most pessimistic states have been utilized to predict the future.

One of the main problems in using the output of the AOGCM models is the large temporal and spatial scale of the computational cell, relative to the study area.

In this research, the SDSM model has been used to downscale the rainfall, maximum temperature, minimum temperature, and wind speed data of the HadCM3 model under two emission scenarios of RCP2.6 and RCP8.5 for the future periods (2020–2049). The SDSM software provides a relation with the maximum correlation coefficient value between the time series of the daily rainfall and the large-scale regional variables. Then, it predicts and generates the time series of the daily rainfall by utilizing the parameters achieved from the above relation and the output of the HadCM3 model under two scenarios of the RCP2.6 and RCP8.5 for the future period. The SDSM model inputs consist of the daily observational data of the station (predictable), large-scale variables of national centers for environmental prediction (NCEP) (predictor) in a similar observational period, and the output of the large-scale variables of GCM under the different emission scenarios for the future study period. The ECMWF database with the spatial resolution of 0.125° for the period 1971–2001 was used to provide the daily observation data. Also, the large-scale data of NCEP for a similar period were downloaded. The SDSM model output is utilized to estimate future runoff. Therefore, these outputs are considered as the VIC model input. As will be mentioned in the following sections, the VIC model needs the networked input data for generating future runoff. Thus, the study area was divided into 29 cells of 0.125 × 0.125°, and the downscaling operations were performed separately for each cell. Therefore, in this study, the SDSM model was performed once for each cell and each parameter (rainfall, maximum temperature, minimum temperature, and wind speed) and one output was obtained for each of them.

The large-scale hydrological models which simulate the hydrological processes of the land surface in the large-scale basins have been considerably developed in recent decades. These models are very close to the models of the land surface utilized in the GCMs, but their main purpose is to model the runoff, river flow, and basin water balance. The large-scale models act in the distributive form because of the physical basis and by considering the major effective processes in the basin energy and water balance, they have the ability to simulate some issues like the spatial distribution of runoff in the basin, spatial distribution of evaporation and transpiration, spatial distribution of soil moisture, soil temperature, etc. Among different models existing in this field, the VIC model is more popular and useable among researchers around the world (Das et al. 2018; Liu et al. 2018; Wang et al. 2019; Gou et al. 2020).

Figure 4 represents an overview of the VIC model, along with three soil layers. The surface of each cell has been covered with N + 1 different vegetation covers and n = N + 1 represents bare soil (Liang et al. 1994). In this model, the evaporation and transpiration are estimated using the Penman–Monteith equation. The evaporation and transpiration is a function of the net radiation and vapor pressure deficit.

The VIC model is one of the large-scale land surface models that simulate the water and energy balance in some cells with the dimensions of 0.125–2° and in the time steps from hours to 1 day. In this research, the water balance has been used.

Although most of the VIC model parameters are based on remote sensing or field observations, some of them cannot be measured by direct measurement. These parameters are conceptual and sometimes have non-physical values. Therefore, they need to be calibrated and validated.

The routing model estimates the time of concentration required for the runoff to reach the output of a cell and also the time of concentration of the channel direction flow (in the channel network). In this model, it is assumed that the main horizontal flow inside the cell enters the available channel network inside the desired cell before reaching the side cell. According to the D8 algorithm, eight different directions can be considered for each cell's output flow, but it should be considered that there is only one output direction for the flow in a cell. When water flows into the channel, it is assumed that there is no possibility of returning water or outflow from the channel; consequently, once it reaches the channel, it is no longer part of the water budget. The daily surface runoff generated by the VIC model in each cell is first transferred to each cell's output using a triangular unit hydrograph and then routed in the channel network to reach the basin output. Figure 5 shows how water transferred in the channel, which is provided using a simple linear routing model that follows the linear form of the Saint-Vénant equation. In this model, it is assumed that all runoff exits a cell in a single flow direction.

The VIC and routing models run in series, which means that the output of the VIC model (fluxes) is one of the fundamental inputs of the routing model.

Since the VIC model's output is the daily runoff data for the future period, it is necessary to use the appropriate method to convert this into the maximum instantaneous flood. In this study, an attempt was made to use both empirical and statistical methods, so three well-known empirical methods (Fuller, Sangal, and Fill Steiner) and the ANN's statistical method were used.

Fuller (1914) studied flood data of 24 watersheds in the United States with basin areas between 3.06 and 151,592 km². Its formula is the most important and widely accepted due to its simplicity.

Sangal (1983) tested his method using data from 387 gauging stations in Canada for basin areas measuring less than one to more than 100,000 km². He obtained good estimations in the majority of basins, although in small basins the peak flow could be underestimated.

Fill & Steiner (2003) created a study based on Sangal's formula. They used data from 14 stations of basins with drainage areas between 84 and 687 km² in Brazil and developed a simple formula, suitable for drainage areas from 50 to 700 km².

In this study, first all four methods were utilized for the past data (1984–2013). Afterwards, the output values of the methods were compared with the actual values and the best method was selected. Finally, the selected method was used to generate the maximum instantaneous flood from the daily runoff for the future period. Table 2 describes the equations of empirical methods.

‘ANN is one of the popular computational intelligence-based approaches in flood prediction providing high accuracy, long lead-time and low computational cost’ (Fotovatikhah et al. 2018). The ANN is the simulation of the natural nervous system, and it consists of a set of neural units called neurons, which is connected by some connections called axons. In ANNs, it is attempted to design a structure similar to the biological structure of the human brain and the nervous system to have the learning, generalizing, and decision-making abilities the same as the biological structure. In these networks, the purpose is to train the model by introducing the performance history of a dynamical system and how the system works is stored in the memory to be used in issues that the model has not previously encountered. The general structure of the ANN is made up of three layers with different tasks. The input layers have the data distribution role in the network. The middle or hidden layer has the role of processing information. The output layer, in addition to processing per input vector, also controls its output network.

One of the most common neural network architectures is called the Multilayer Perceptron (MLP) (Nazzal et al. 2008). These networks are made up of different layers. Some neurons are considered in each layer, which are connected by some connections to the neurons in adjacent layers. The neurons of the first layer take the input information and transfer them to the hidden layer's neurons using the related connections. In these networks, each neuron's effective input is the multiplication of the outputs of the previous layer neurons in the weights among those neurons.

In neural networks, data is generally divided into three parts: training, validation, and testing. Normally, in the training part, the model trains the network with the training set to adjust the weights. To avoid overfitting the network and also fine-tune models, it is essential to input the validation set to the network and check if the error is within some range. This set is not being used directly to adjust the weights but is used to give the optimal number of hidden units or determine a stopping point for the back-propagation algorithm. Finally, the accuracy of the model on the test data gives a realistic estimate of the performance of the model on completely unseen data and in order to confirm the actual predictive power of the network.

In the current research, NeuroSolutions (version 5) software with MLP architectures (two layers), 10 neurons and the feed-forward back-propagation function has been utilized. In designing this network, the LM training function (this function updates weight and bias values according to Levenberg–Marquardt optimization) has been used due to it being the fastest medium-sized backpropagation algorithm. Also, the data of maximum instantaneous flood and daily runoff in the period 1984–2013 were considered as input data. Seventy per cent of the data were selected for training, 15% for validation and 15% for testing.

In this research, to evaluate the performance of models, the statistics shown in Table 3 were used.

In this relation, O_i is the observational values, S_i is the simulated values, the average of the observational values, is the average of the simulated values and n is the number of data.

In trend analysis it is crucial to assess whether or not the time series data are independent. If a serial correlation is present in the time series, it can significantly affect the trend analysis outcome. Positive autocorrelation can artificially induce a trend in a time series, while negative autocorrelation can weaken the trend if it is present. In this research, before analyzing the trend in the maximum instantaneous floods, these values should be tested for serial correlation. For this purpose, the Lagrange multiplier (LM) test is used. The Breusch–Godfrey test is a test for autocorrelation in the errors in a regression model. It makes use of the residuals from the model being considered in a regression analysis, and a test statistic is derived from these. The null hypothesis is that there is no serial correlation of any order up to lag order p, where p is a pre-specified integer. Because the test is based on the idea of LM testing, it is sometimes referred to as an LM test for serial correlation.

In order to analyze the trend in the maximum instantaneous flood data, trend tests are utilized. These tests consist of the nonparametric tests including Mann–Kendall and Spearman and the parametric test including the linear regression test. In this study, Trend version 1.0.2 software was used. The input data used in this software include the maximum instantaneous flood data for the past periods (1984–2013) and future periods (2020–2049).

The positive values of S demonstrate the ascending trend, and its negative values demonstrate the descending trend. Also, the critical values of the statistic z are available in the normal probability tables for the different significance levels.

For large samples, the value of has nearly normal distribution with the average of zero and the variance is one. The test statistic's critical values are available in the normal probability tables for the different significance levels.

The test statistic S follows a student's t-distribution with a degree of freedom of n–2. The test statistic's critical values are available in the statistical tables of student's t-distribution for the different significance levels. The linear regression test assumes that the data has the normal distribution and independent errors, and they follow a normal distribution with zero average.

The frequency analysis of the maximum instantaneous flood is mostly based on the stationary assumption and designing the hydrological structures is done according to the stationary conditions. However, since natural phenomena are in non-stationary conditions in practice, it is necessary to handle the issues in a non-stationary way. Also, it should be noted that the effects of climate change are increasing and, therefore, the non-stationary conditions will be increased as well. Consequently, new methods should also be discovered to analyze the frequency of the maximum instantaneous flood using the non-stationary data.

The NEVA software package has been used to achieve the goal of this study, which is to study the non-stationary conditions of the maximum instantaneous flood in the future (2049–2020) and the past period (1984–2013).

The software package of NEVA has been improved to make the extreme value analysis under stationary and non-stationary assumptions easy. In a Bayesian approach, NEVA estimates the extreme value parameters with a Differential Evolution Markov Chain (DE-MC) approach for global optimization over the parameter space. In NEVA, the time-variant parameter (μ_t) can be derived using different quantiles from the DE-MC. In this paper μ_t is computed as the median of μ_t and the 95th percentile of the DE-MC sampled μ_t values (low risk).

NEVA includes posterior probability intervals (uncertainty bounds) of estimated return levels through Bayesian inference, with its inherent advantages in uncertainty quantification. The software presents the results of non-stationary extreme value analysis using various exceedance probability methods. Generally, NEVA has two sections:

• 1.

  GEV distribution for analysis of annual maxima (block maxima).

• 2.

  Generalized Pareto distribution (GPD) for analysis of extremes above a certain threshold (i.e. peak-over-threshold (POT) method).

The GEV and GPD approaches can be utilized to analyze the stationary (independent of time) and non-stationary extreme values. The GEV method has been utilized in this study.

In this relation, μ is the location parameter which determines the distribution center, ɛ is the shape parameter which has effects on the tails of the GEV distribution, and σ is the scale parameter which determines the amount of deviation around the location parameter.

Now, if ɛ is close to zero then it is the Gumbel distribution, if ɛ > 0 it is the Frisch distribution and if ɛ < 0 it is the Weibull distribution. In this study, the extreme events do not change with time in the stationary condition, and in the non-stationary process, the GEV distribution parameters are variable with time.

NEVA identifies the trend and non-stationary existences in the extreme values of input data using the Mann–Kendall test and the significance level is determined by the user. By default, this significance level is equal to 0.05, which is hugely useful in hydrologic studies. If the null hypothesis is rejected, NEVA provides the extreme value analysis under the stationary assumption. If the trend is recognized in the level of 5%, the GEV distribution parameters are estimated under the non-stationary assumption (μ₁, μ₀, σ, and ɛ). NEVA utilizes a Bayesian technique to comprehend the GEV distribution parameters under the stationary and non-stationary conditions.

Because the output of the AOGCM models is large scale, it should be downscaled. The SDSM model has been used to downscale the rainfall, maximum temperature, minimum temperature, and wind speed data of the HadCM3 model under two emission scenarios of RCP2.6 and RCP8.5 for future periods (2020–2049). For calibration of the model, ECMWF daily data for the period of 1971–2001 were used as predictable and NCEP daily data for the period of 1971–2001 were used as predictor. The model performance in calibration for a sample cell (the first cell from the bottom left in Figure 2) is shown in Table 4.

As is obvious in Table 4, the minimum temperature performed better compared to the other two parameters with R² of 0.85 and RE of 3.05%.

After calibrating the SDSM model for all the study area cells, the values of temperature, rainfall, and wind speed for the period of 2020–2049 were simulated for all the cells. The simulation results of minimum temperature and rainfall for the sample cell are represented in Figures 6 and 7.

According to Figure 6, the minimum temperature increase is generally observed in the future compared to the past period for all months. Also, the difference between past and future minimum temperatures in summer is greater than in other months. In both the past and the future, July is the warmest month. The increase reached its maximum level in July, which demonstrates an increase of about 5 °C in the future compared to the past. Also, the lowest difference is in November, which is close to 1 °C.

As is obvious in Figure 7, the rainfalls in the future will be decreased in all months. Also, the difference between past and future rainfalls in spring is greater than in other months. The maximum reduction is related to March and April, which is more than 50 mm in the future compared to the past. The minimum amount of rainfall reduction is in January, which is close to 2 mm.

It can be concluded that the rate of decrease in precipitation will be higher than the rate of increase in temperature in the future period. Also, in all cells and most of the months, the temperature increase and rainfall reduction are predictable in the future compared to the past.

The model calibration and validation were performed for the period between 2010–2012 and 2012–2013, respectively. Also, the Doab Mark station's daily runoff was utilized to compare the model's output data in calibration and validation. Table 5 shows the model assessment.

According to Table 5, the results demonstrate the proper performance of the model and its proper reliability in the runoff estimation. Figures 8 and 9 represent the simulated and observed flow hydrograph during the calibration and the validation, respectively.

Figures 8 and 9 demonstrate that the model mostly has an acceptable ability to simulate the flows, and the model weakness is in simulating the peak flow.

Table 6 represents the values of the soil parameters achieved from the calibration stage of the VIC model and Table 7 represents the description of each parameter.

As is obvious in Table 6, the Ds value is small, which represents a low base discharge in the lower layer. The smaller value of bi shows an increase in the infiltration and a decrease in the surface runoff.

Figure 10 represents the long-term average monthly runoff at the Doab Mark hydrometric station in the Ghareh Sou basin for the past period (1984–2013) and the future period (2020–2049) under two emission scenarios of RCP2.6 and RCP8.5 (obtained from the VIC model).

According to Figure 10, the decrease in runoff is seen in all months, in the future period under both scenarios. The runoff reduction in March and April reaches its maximum level (more than 30% decrease in scenario RCP2.6 and more than 50% decrease in scenario RCP8.5) and in September and August reaches its minimum level under both scenarios (more than 1% decrease in scenario RCP2.6 and more than 5% decrease in scenario RCP8.5). It seems that the significant decrease in runoff in spring is mainly due to the decrease in rainfall in spring. Generally, it can result that increasing the temperature and decreasing the rainfall in the future compared to the past will reduce the runoff in the future.

Additionally, the runoff under scenario RCP8.5 is less than the runoff under scenario RCP2.6 in all months. The highest difference occurs in March, and the lowest difference in June. In other words, a further increase in greenhouse gas emissions will have a greater impact on decreasing the amount of runoff in the future in the study area.

To estimate the maximum instantaneous flood, the daily discharge data of four methods of Fuller, Sangal, Fill Steiner, and the ANN were utilized. Table 8 represents the evaluation of these four methods for the past period (1984–2013).

According to Table 8, all four methods have acceptable results. However, due to the fact that the ANN method generated the most appropriate result, this method has been selected for further research.

Table 9 shows the results of the LM test.

Based on the results of Table 9, the probability of F statistic is more than 5% and it can be concluded that the null hypothesis is accepted and there is no serial correlation in the maximum instantaneous floods data (both in the past and future data).

In order to analyze the trend in the maximum instantaneous flood data for the past and the future periods under two emission scenarios of RCP2.6 and RCP8.5, three tests of Mann–Kendall, Spearman, and linear regression were provided. Table 10 demonstrates the results of these three tests on the maximum instantaneous flood.

According to the results of Table 10, the values of the maximum instantaneous flood discharge in the past period have a significant nonlinear trend at the medium significance levels (0.05) and it has a significant linear trend at the lower significance levels (0.1). There is a trend in the maximum instantaneous flood of the region in the past, and therefore they are in the non-stationary form. The results also demonstrate that the maximum instantaneous flood values in the period of 2020–2049 under the effect of climate change have a significant trend at high levels (0.01). Comparing the results of the maximum instantaneous flood trend under the RCP2.6 scenario with RCP8.5 shows that the trend rate in the RCP8.5 scenario was more than the RCP2.6 scenario, In other words, further increase in the greenhouse gas emissions has a greater impact on the trend rate of the maximum instantaneous flood discharge.

Another point obtained from Table 10 is that the instantaneous flood discharge in both the past and future has a negative trend, and the negative trend will increase in the future. In other words, the region will observe a reduction in the maximum instantaneous flood intensity in the future.

The purpose of this research is to study the non-stationary conditions of the maximum instantaneous flood in the past and future period under the two scenarios of RCP2.6 and RCP8.5. For this reason, the NEVA software package was utilized. This software package was run in two stationary and non-stationary forms.

In this section, the outputs of the NEVA software package estimating the maximum instantaneous flood relative to the different return periods in the stationary and non-stationary form under two scenarios (RCP2.6 and RCP8.5) are discussed.

Figure 11 shows the maximum instantaneous flood in different return periods under the stationary assumption, while Figure 12 exhibits the maximum instantaneous flood under the non-stationary assumption for the observation period (here, the past period is 1984–2013 and the future period is 2020–2049). Figure 13 displays the maximum instantaneous flood under non-stationary conditions for 100 years beyond observations using median and 95th percentile of sampled location parameters, respectively. Also, the uncertainty bounds of the computed maximum instantaneous flood can be derived based on 5 and 95% posterior probability intervals of the ensemble.

In hydrology studies, design return periods vary typically from 10 to 100 years. Table 11 represents the maximum instantaneous flood for different return periods (10–100 years) under the stationary and non-stationary assumptions in the past and future periods for both scenarios of RCP2.6 and RCP8.5 (this table is a summary of Figures 11–13).

Table 11 is designed to answer whether there is a difference between estimating the return period of the maximum instantaneous flood in stationary and non-stationary conditions.

According to Table 11, the maximum instantaneous floods under the stationary assumption are mostly smaller than the maximum instantaneous floods under the non-stationary assumption in both the past and future periods. In most cases, the difference increases by expanding the return period. The results also show that the difference between the maximum instantaneous flood under the stationary and non-stationary assumptions in the future is more than in the past period under both scenarios. The maximum difference in the past period is 9.3%, and in the future, under scenarios RCP2.6 and RCP8.5, are 30.5 and 36.8%, respectively. Moreover, the maximum instantaneous floods in all return periods and scenarios in the future are smaller than the maximum instantaneous floods in the past. Also, the maximum instantaneous floods under the RCP8.5 scenario in all return periods are smaller than under the RCP2.6 scenario.

As the results of the trend test on the maximum instantaneous flood data in the past and future periods show, more trends will occur in the future period under both scenarios compared to the past period in the study area which indicates that the difference of the maximum instantaneous flood between the stationary and non-stationary assumptions is perhaps for this reason.

The current research was an attempt to study the non-stationary condition on the maximum instantaneous flood in the future period in the Ghareh Sou basin. Thus, the VIC model was used to generate the runoff of the future periods (2020–2040) under two emission scenarios of RCP2.6 and RCP8.5. The results obtained from the calibration (Nash coefficient 0.82) and validation (Nash coefficient 0.66) of the model indicate its proper performance. This has also been proven in the research of Aziziyan et al. (2017). They calculated the Nash coefficient of the VIC model in the calibration stage as 0.84 and the validation stage as 0.74. The main weakness of the model was in simulating the peak flow. This study calibrated the most sensitive parameters of the VIC model while the other parameters were not calibrated. For example, the LAI and canopy fraction have a strong influence on results (Bennett et al. 2018). The source of these parameters is the remote sensing data (satellites) and for calibration of them, it needed to network to ground-based observation data which is not available.

The ANN method was utilized to achieve the maximum instantaneous flood from the daily runoff. Finally, the stationary and non-stationary analysis were provided on data using the NEVA software package and using the GEV method. In the non-stationary analysis, it was assumed that the location parameter (μ) has a linear relationship with time. In future studies, the relationship between the location parameter and time can be assumed to be nonlinear and also the impacts of this assumption on the results can be examined. Furthermore, the variability of the scale parameter with time (as linear and nonlinear) can also be investigated. The effect of simultaneously considering the assumption of non-stationarity in the scale and shape parameters on the results can be examined.

Considering the non-stationary assumption in the analysis of maximum instantaneous flood is the most important advantage of the present study compared to other studies in this field (Ganamala & Sundar Kumar 2017; Lima et al. 2017; McCollum et al. 2019). On the other hand, the lack of access to long-term data to investigate the non-stationarity of the shape parameter can be considered as one of the limitations of the present study.

According to the results of downscaling the HadCM3 data, the temperature increases while rainfall decreases in the future (2020–2049), in contrast to the past period (1971–2001). Also, the VIC model results show that the daily runoff in the future period of 2020–2049 will reduce compared to the past period of 1984–2013. The reduction is also obvious in the peak values and the maximum instantaneous flood. In this regard, many studies have been carried out in Karkheh basin. For example:

• 1.

  Malmir et al. (2016) used the HadCM3 model under the A2 scenario. The results indicate an increase in temperature and a decrease in precipitation in the period 2011–2040 compared to the period 1971–2000 and a 32.62% decrease in discharge in the region.

• 2.

  Adib et al. (2020) used the HadCM3 model under the A2 and B2 scenarios. The results of this research indicate an increase in temperature and a decrease in precipitation and discharge in the future period (2030–2073).

• 3.

  Ahmadi et al. (2020) used the CORDEX model under two scenarios (RCP4.5 and RCP8.5) for projection of precipitation and temperature data for 2020–2090. The results showed that under RCP8.5 and RCP4.5 the amount of annual precipitation will decrease by 88.2 and 56.1 mm, respectively. Also, according to RCP8.5 and RCP4.5, the minimum temperatures will rise by 4.5 and 2.5 °C, respectively.

The results of these studies are similar to the results of the present study.

Based on the results of Table 11, it can be concluded that considering the non-stationary assumption can have an impact on the maximum instantaneous flood values under different return periods. Numerous studies have been performed around the world to prove the effect of non-stationary on extreme values. For example, Faulkner et al. (2019), using the Ballysadare River annual maximum flow data from 1945 to 2015, examined the 100-year return period under stationary and non-stationary conditions. In situations where the location parameter was variable with time, the estimated design flow increases linearly over the period of record. Jake & Richard (2017), estimated 100-year streamflow quantile for Aberjona River using the stationary and nonstationary approaches. They estimated the streamflow of 1,570 m³ for stationary conditions and 2,170 m³ for non-stationary conditions.

Considering the non-stationary assumption for maximum instantaneous flood analysis is a new issue in the study area that has not been examined before. Due to the existence of water dams, urban areas, rural areas, and arable lands in the study area, predicting the return period values of the maximum instantaneous flood in the future period is important for the construction of water structures and water resources management. According to the results of NEVA, the lack of consideration of the non-stationary condition and ignoring the trend in the maximum instantaneous flood data causes errors in estimating the maximum instantaneous flood in different return periods (underestimating values) which increases with an increasing return period. This creates errors in designing hydraulic structures and causes financial and human losses.

All relevant data are included in the paper or its Supplementary Information.