Assessing future changes in flood frequencies under CMIP6 climate projections using SWAT modeling: a case study of Bitlis Creek, Turkey

Floods are dynamic environmental threats that can lead to serious economic, societal, and ecological repercussions, including the loss of lives (Klijn et al. 2015; Bai et al. 2019). Many of these consequences of flooding stem from the destruction of civil engineering infrastructures. The assessment of flood risk for these structures commonly relies on the magnitude and frequency of flood occurrences (Mani et al. 2014). For instance, when hydraulic structures, dams, or similar constructions are to be built on a river, flood hydrographs for anticipated floods with frequencies of 5, 10, and 25 years are used in the design of diversion facilities that control streamflow and mitigate flooding during the construction period. The magnitude of rarer events, such as 50- or 100-year floods, is crucial in the design of structures like bridges and culverts located at the outlet of small and medium watersheds. In situations with a significant risk of loss of life, it is appropriate to design a facility to withstand 500-year or even more extreme flood events. Civil engineers typically rely on historical streamflow data when designing water-related infrastructures (Quintero et al. 2018). However, it is anticipated that future climate will increase the likelihood of occurrence and strength of extreme rainfall events, leading to more frequent river flooding and heightened flood risks (Fleming et al. 2014; Shrestha & Lohpaisankrit 2017). Previous studies have demonstrated that flood intensity is highly sensitive to changes in climate (Kure & Tebakari 2012; Prudhomme et al. 2013; Quintero et al. 2018; Meresa et al. 2022; Alaminie et al. 2023; Kiran et al. 2023). In the coming decades, global climate change will undoubtedly become a predominant factor altering hydrological regimes and flood risk conditions. The nature of those alterations varies locally, underlining the necessity for regional evaluations to guide flood risk management and planning (Zhang et al. 2023).

Developing mitigation and adaptation strategies in response to climate change relies on the use of robust simulations of climate change scenarios (Eyring et al. 2019). While global circulation models (GCMs) serve as state-of-the-art tools for climate change analyses through predefined scenarios of societal development, their direct use in predicting local-scale future flood magnitudes and frequencies is hindered by uncertainties stemming from multiple sources. Climate projections of GCMs include various systematic errors, commonly referred to as bias (Wang & Chen 2014). Researchers employ diverse techniques to reduce model biases (Li et al. 2010; Teutschbein & Seibert 2012; Pierce et al. 2015; Yuan et al. 2016). Huang et al. (2011) proposed utilizing a minimum of 30 years of time series data for bias correction to establish a statistically relevant and consistent connection between past and future scenarios. However, due to the use of unique simplifications, parameterizations, and numeric approximations related to climate systems in different GCMs, there exists model uncertainty in long-term climate simulations (Murphy et al. 2004; Knutti et al. 2010, 2019). Several approaches have been suggested in the literature to quantify and minimize the uncertainty arising from GCM selection (Shiogama et al. 2016; Hosseinzadehtalaei et al. 2017; Mateus & Tullos 2017; Kaczmarska et al. 2018; Lehner et al. 2019). While most of these methods rely on comparing individual GCMs to their ensemble mean, it is important to note that the return period of extreme events heavily depends on climate variability rather than climatic mean (Salman et al. 2020). Therefore, uncertainty in climate projections can significantly under- or overestimate the return periods of extremes (Harris et al. 2013; Hewitt et al. 2016). Wang et al. (2020) emphasized that GCMs are the leading sources of uncertainty in evaluating changes in annual and peak streamflows, and they recommended using at least 10 GCMs when assessing the impacts of climate change on hydrology to thoroughly quantify uncertainty.

Streamflow projections and flood frequency analysis in the context of a changing climate are crucial for strategic planning and management of water resources, especially in highly susceptible snow-fed watersheds (Adam et al. 2009). The Euphrates and Tigris are the two major snow-fed rivers of the Middle East region, defined as a prominent climate change hotspot warming almost two times more rapidly than the global mean (Giorgi 2006; Zittis et al. 2022). Bozkurt & Sen (2013) proposed that the snow-dominated highlands of the Euphrates-Tigris Basin in Turkey are expected to face more pronounced direct consequences of climate change compared to the other riparian countries. According to Sen et al. (2011), Bozkurt et al. (2015), and Yucel et al. (2015), projections indicate shifts toward earlier times of year in snowmelt runoff timings for the northern highlands of the watershed. Despite the anticipated changes in the future streamflow regime of the basin, as indicated by these and several other studies (e.g., Nohara et al. 2006; Kitoh et al. 2008; Özdoğan 2011; Şen 2019; Peker & Sorman 2021; Şensoy et al. 2023) based on GCM projections from the third and fifth phases of the coupled model intercomparison project (CMIP), the situation could be worse than expected according to the latest release of CMIP, the sixth phase (CMIP6). Bağçaci et al. (2021) compared temperature and precipitation projections for Turkey using ensembles from the best-performing GCMs belonging to the CMIP5 and CMIP6 databases. Their findings indicate that CMIP6 projections foresee a slightly warmer climate than CMIP5, and CMIP6 exhibits more noticeable extreme precipitation changes in both maximum decrease and increase, addressing an improved representation of future climate variability. To the best of the knowledge, there is no specific study in the literature quantifying the potential impacts of climate change on return levels of flows for the most threatened highlands region of the basin under the CMIP6 climate scenarios. Such quantification enables a better understanding of any potential under- or overestimation of flood discharges when historical streamflow records are traditionally used in the design of water resources projects.

Bitlis Creek is one of the main branches of the Tigris River within the borders of Turkey, along with Garzan, Botan, and Batman Creeks. It originates from the southern slopes of the Southeast Taurus Mountain Range and gains the character of a river at elevations of around 1,500 m during its southward course toward the Tigris River. Along this flow route, it conjoins with Kezer Creek, its largest branch, at about 500 m elevation, and later merges with Botan Creek at an elevation of approximately 450 m before flowing into the lake of the Ilisu Dam situated on the Tigris River (Figure 1). Although the Bitlis Creek basin experiences a continental climate overall, the upper part of the watershed falls within a distinctive micro-climate zone distinguished by abundant snowfalls, leading to snow depths of up to 5 m (Aydin & Işhik 2015). Five cascade hydroelectric power plant (HEPP) projects (i.e., Deliktas Weir and HEPP, Bitlis Weir and HEPP, Baykan Dam and HEPP, Baykan II Weir and HEPP, and Basoren Dam and HEPP, from upstream to downstream) have been planned on the mainstream of Bitlis Creek to utilize the hydropower potential of this snow-fed mountainous watershed, draining an area of about 1,076 km² from its source to the conjunction with Kezer Creek (DSI 2023). Bitlis Creek is selected as the case for this study due to its long-term historical flow records, which accurately reflect natural flow conditions in the basin, given that the basin's dam projects are not yet in operation.

There are two stations measuring the runoffs of Bitlis Creek, namely Karinca (Station ID: E26A016) and Baykan (Station ID: E26A010) (Figure 1). While the Karinca SGS, draining an area of 346.4 km², has 5-year daily flow data from 1965 to 1969, the Baykan SGS has a drainage area of 636.5 km², and its 56-year available daily flow data cover the period from 1955 to 2010 (DSI 2022). The streamflow measurements of the Baykan SGS are the fundamental input for the reconnaissance and feasibility assessments of the hydropower projects on Bitlis Creek in optimizing power plant and reservoir storage capacities, as well as in determining spillway design discharges (EIE 1990; Yolsu 2009).

The synoptic meteorological stations (MSs) with daily precipitation, maximum and minimum air temperature, wind speed, relative humidity, and solar radiation records inside or near the Baykan SGS basin are Bitlis (Station ID: 17207) and Siirt (Station ID: 17210) located at altitudes of 1,573 and 895 m, respectively (Figure 1). It is a standard procedure to use all available historical streamflow data in flood frequency analysis (Quintero et al. 2018). The precipitation and temperature records of the Bitlis and Siirt MSs allow the historical analysis period of the study to be from 1955 to 2010. On the other hand, the hydrological modeling timeline of the study is constrained to the period from 1988 to 2009 due to the availability of wind speed, relative humidity, and solar radiation records from these stations (MGM 2022b).

The daily precipitation, maximum air temperature, and minimum air temperature simulations of the 21 CMIP6 GCMs for the coordinates of the Bitlis and Siirt MSs are utilized in forcing the hydrologic model to project daily streamflow time series at the Baykan SGS location for the periods 1955–2010 and 2025–2099. The grid-based GCM datasets for the CMIP6 historical experiment and future scenarios of SSP245 and SSP585 are retrieved from the Earth System Grid Federation (ESGF) website (ESGF 2022). The model IDs, institutional information, and horizontal resolutions of the selected GCMs from the CMIP6 database are listed in Table 1. The model selection is based on the availability of daily outputs under the first ensemble member for the considered climate scenarios. Several ensemble members with different variant labels, constructed using the indices of realization (r), initialization (i), physics (p), and forcing (f), are provided by the CMIP6 models, and a few models have more than one ensemble member (Sun et al. 2022). To ensure equal weighting for different CMIP6 GCMs, this study considers only the simulation results of the first ensemble member (i.e., r1i1p1f1) for each GCM. Additionally, the climate datasets with different grid sizes of the selected GCMs are uniformly interpolated to a spatial resolution of 0.5° × 0.5° under the first-order conservative remapping method (Jones 1999) before extracting the raw simulation data for the coordinates of the Bitlis and Siirt MSs.

The hydrologic model of the Baykan SGS basin is constructed with SWAT, developed by the United States Department of Agriculture – Agricultural Research Service (USDA-ARS) (Neitsch et al. 2011; Arnold et al. 2013). The SWAT model setup is carried out within ArcSWAT 2012 (revision 664), using the void-filled 1 arc-second resolution digital elevation model (DEM) obtained from the Shuttle Radar Topography Mission (SRTM) database (USGS 2014), the 1 km spatial resolution land cover map extracted from the Global Land Cover 2000 (GLC2000) version 2.0 dataset (EC-JRC 2006), and the soil map at a scale of 1:5 million derived from the Digital Soil Map of the World (DSMW) version 3.6 dataset (FAO 2007) (Figure 1). In the SWAT modeling, the delineated watershed is firstly divided into 81 subbasins based on the inputted DEM data. Then, these subbasins are further divided into hydrological response units (HRUs) based on the inputted grid-based land cover and soil characteristics of the basin, along with the applied topographic slope discretization considering five slope classes (i.e., 0–5%, 5–15%, 15–25%, 25–50%, and >50%). To maximize computational efficiency by eliminating minor HRU formations, the threshold percentages of 5% for land cover, 5% for soil, and 15% for slope are utilized, resulting in a total of 586 HRUs.

SWAT utilizes station-based daily weather data for precipitation, maximum and minimum air temperatures, wind speed, relative humidity, and solar radiation to conduct runoff simulations at the HRU scale, which are then aggregated at the subbasin scale and routed downstream. Considering the observation periods of the Bitlis and Siirt MSs, the 1984–2009 timeframe is identified as the common period to calibrate and validate the SWAT model against the streamflow measurements of the Baykan SGS, including 4 years for the model warm-up (i.e., 1984–1987). Hence, the daily records from the Bitlis and Siirt MSs for the years 1984 to 2009 are inputted into the model, along with the geographic position data and long-term mean monthly weather statistics of the stations (MGM 2022a, 2022b, 2022c). Accordingly, SWAT assigns the climate data of the Bitlis and Siirt MSs to 71 and 10 subbasins, accounting for 87 and 13% of the overall catchment area, respectively. Moreover, five elevation bands are defined for all subbasins, except those with minimal maximum reliefs, to account for the impact of the mountainous landscape on precipitation and temperature when simulating snow-related processes (Abbaspour 2015).

The developed SWAT model is calibrated and validated against the daily historical streamflow time series from the Baykan SGS (DSI 2022) for the period from 1988 to 2009 within the SWAT Calibration and Uncertainty Procedures (SWAT-CUP) software package (Abbaspour 2015). The first 14 years of the records (i.e., 1988–2001) are utilized to calibrate the model for daily streamflow simulations, and the records of the subsequent 8 years (i.e., 2002–2009) are used to verify the accuracy of the calibrated model. For the parameter sensitivity, calibration, validation, and uncertainty analyses, the Sequential Uncertainty Fitting Version 2 (SUFI-2) algorithm (Abbaspour et al. 2004, 2007) is employed as the optimization procedure using the bR², defined as the coefficient of determination (R²) multiplied by the slope of the zero-intercept linear regression line between observed and simulated datasets (b), as the objective function.

In the SUFI-2 algorithm, uncertainties originating from the conceptual model, input datasets, and modeling parameters are iteratively distributed over the ranges of sensitive parameters. The model output uncertainty is quantified by the 95% prediction uncertainty (95PPU) band, which is produced at the 2.5 and 97.5% levels of the cumulative distribution of an output variable obtained through Latin hypercube sampling, in terms of the statistics of P-factor and R-factor. The P-factor indicates the proportion of observed data bracketed by the 95PPU band, and the R-factor represents the ratio of the average width of the 95PPU band to the standard deviation of observed data. Since there is a trade-off between these two indices, achieving a higher P-factor necessitates accepting a larger R-factor. To strike a balance, starting with wide parameter ranges, the SUFI-2 algorithm is iterated a few times by gradually narrowing the assigned parameter ranges at each iteration until an optimal set of parameter ranges is achieved. Accordingly, the final iteration has the calibrated parameter ranges, and the best simulation of the final iteration with the highest objective function value (i.e., bR²) has the best-performing parameter values (Abbaspour et al. 2015).

In this study, the SWAT model is calibrated in three steps against the daily streamflow records of the Baykan SGS in the 1988–2001 period, using the weather datasets from 1984 to 2001, to avoid parameter interaction and identifiability problems in simulating snowpack and snowmelt processes: (1) calibration of the precipitation lapse rate (PLAPS) and temperature lapse rate (TLAPS) parameters in the subbasin (.sub) input files and fixing them to their best simulation values, (2) calibration of the snow-related sensitive parameters in the basin (.bsn) input file and fixing them to their best simulation values, and (3) calibration of the other streamflow-related sensitive modeling parameters in the groundwater (.gw), HRU (.hru), main channel (.rte), soil (.sol), and management (.mgt) input files. At each step of the applied calibration procedure, the sensitivity of each parameter is evaluated individually by performing a one-at-a-time analysis through a single iteration with 50 simulations. After determining the sensitive parameters and their initial ranges in these analyses, combined iterations with 500 simulations are repeated for the sensitive parameters of that step by narrowing the initially assigned parameter ranges each time until acceptable P-factor and R-factor values are obtained (Abbaspour et al. 2004; Yalcin 2019).

The calibrated parameter ranges and best-performing parameter values are verified against the streamflow records in the 2002–2009 period by performing a single combined iteration with 500 simulations and a single simulation, respectively, using the weather datasets from 1998 to 2009. To evaluate the streamflow forecasting performance of the SWAT model under the use of the best-performing parameter values, the statistics of R² (with b), Nash–Sutcliffe efficiency (NSE) (Nash & Sutcliffe 1970), percent bias (PBIAS), and the ratio of the root mean square error to the standard deviation of observed data (RSR) are utilized in addition to the daily time series graphics.

GCM simulations require bias adjustment to rectify potential systematic errors and align them with specific observed climate characteristics, such as distribution, sequencing, and magnitude (Tan et al. 2020). Hence, the extracted GCM datasets for the coordinates of the Bitlis and Siirt MSs are bias-corrected against the daily records of the stations before utilization as input weather data for the calibrated SWAT model. In this process, the linear scaling (LS) method proposed by Lenderink et al. (2007) is applied within the Climate Model Data for Hydrologic Modelling (CMhyd) tool (Rathjens et al. 2016) to reduce possible systematic biases at the point scale for both historical and future periods. The LS method is based on adjusting the GCM simulations by imposing the mean monthly values to match entirely the long-term observed means. While the daily temperature simulations are adjusted by adding the difference between the long-term monthly means of the observed and simulated datasets, the daily precipitation simulations are adjusted through a monthly multiplication factor, keeping the same mean as the observed data.

To estimate the daily streamflow rates at the Baykan SGS location under the climate conditions of the CMIP6 historical experiment in the period from 1955 to 2010, the bias-corrected historical simulations of each GCM for the coordinates of the Bitlis and Siirt MSs in the 1951–2010 period are inputted into the calibrated SWAT model, and a total of 21 model runs are performed using the first 4-year climate data for the model warm-up. Then, the calibrated model runs are repeated 42 times over the 2025–2099 period for the future scenarios of SSP245 and SSP585, using the bias-corrected climate datasets of each GCM for the 2021–2099 period. In all these model runs, the representative daily values of wind speed, solar radiation, and relative humidity, simulated using the pre-inputted monthly weather statistics by the built-in weather generator within SWAT (Neitsch et al. 2011), are utilized to facilitate fair comparisons between the historical and future periods. Before performing flood frequency analysis, the median of the 21 historical streamflow projection time series is analyzed against the station records on monthly, seasonal, and annual time scales to evaluate its appropriateness for use as the baseline scenario in assessing changing climate impacts on the streamflow rates under the SSP245 and SSP585 scenarios.

Annual maximum daily streamflow projections at the Baykan SGS location under the CMIP6 historical experiment and future scenarios of SSP245 and SSP585 are extracted from the 21 daily streamflow projection time series for the 1955–2010 period and 42 daily streamflow projection time series for the 2025–2099 period, respectively. All these sets of annual daily streamflow peaks obtained for each GCM's bias-corrected historical and future climate simulations are utilized in the flood frequency analysis, along with the observed annual maximum daily streamflow rates. Analyzing the frequency of extreme events is crucial for comprehending the likelihood of floods occurring again at specific return periods. Nonetheless, high flow frequency estimates mainly rely on the type of probability distribution model and the number of distribution parameters (Meresa et al. 2022). Since there is currently no universally acceptable probability distribution model, depending on the environmental and climatic background, different distribution functions are deployed to estimate the frequencies of high flow (Ukumo et al. 2023). In this study, the normal, log-normal (with two parameters), log-normal (with three parameters), Pearson type 3, log-Pearson type 3, and Gumbel distributions, mostly conforming to the high flows in the upper Tigris River basin (Yalcin & Tigrek 2016), are utilized in assessing flood discharges for the 5-, 10-, 25-, 50-, 100-, and 500-year return periods.

The first stage of the analysis involves conducting the nonparametric Kolmogorov–Smirnov (K-S) test (McCuen 1993) to assess the appropriateness of the considered distribution functions at a significance level of 5% for each set of annual daily streamflow peaks. Then, the medians of the flood discharges attained under the historical climate data of each GCM are compared with the flood discharges based on the observed annual maximum daily streamflow rates for each distribution function to quantify the overall uncertainty resulting from the GCM simulations and hydrologic modeling. After testing the use of the medians of the large ensemble against the flood discharges based on the observed runoff data in the 1955–2010 period for all distribution functions, the future responses of the flood discharges to changing climate are assessed over the 2025–2074 and 2025–2099 periods under the medium- and high-forcing scenarios of SSP245 and SSP585, considering the median flood discharges attained for the 1955–2010 period as the baseline scenario. Following the distribution function-based comparisons of the median flood discharges, to estimate flood discharges that can be used for design purposes in applications, the same procedure is applied under the use of the best-fitted distribution function according to the K-S test results for each set of annual daily streamflow peaks. For the historical period, the median flood discharges attained under the use of the best-fitted distribution functions are checked against the flood discharges of the most appropriate distribution function for the station records, and the projected changes are quantified again considering the median flood discharges of the historical period as the baseline scenario.

The last combined iteration of the calibration process, performed using the calibrated parameter ranges, provides a 95PPU band with the P-factor and R-factor values of 0.89 and 0.95, respectively. Theoretically, the P-factor can vary between 0 and 1, while the R-factor can range from 0 to infinity. A P-factor value of 1 and an R-factor value close to 0 suggest a precise match between observed and simulated datasets. According to Abbaspour et al. (2015), a P-factor value greater than 0.7 and an R-factor value less than 1.5 (around 1) are considered adequate to establish a balance between these two indices for streamflow simulations. Therefore, the resultant 95PPU band indicates that the parameter uncertainties are within the desired ranges in terms of the achieved P-factor and R-factor values for the study watershed.

Regarding the best simulation of the last combined iteration, the R² (with b), NSE, PBIAS, and RSR statistics are attained to be 0.76 (b = 1.00), 0.68, 4.6%, and 0.56, respectively, as detailed in Table S2 (Supplementary Material). According to the statistical performance evaluation criteria presented by Moriasi et al. (2015), the performance of watershed-scale models can be evaluated as ‘satisfactory’ if an R² value higher than 0.6, an NSE value higher than 0.5, and a PBIAS value less than ±15% are achieved for daily streamflow simulations. Additionally, Moriasi et al. (2007) suggest an RSR value of less than 0.7 for a ‘satisfactory’ streamflow simulation. Accordingly, the statistical performance of the SWAT model under the use of the best-performing parameter values reflects an acceptable agreement between the daily observed and simulated streamflow rates for the calibration period.

The daily simulation results attained for the validation period of 2002–2009 are presented in Figure 3(b), along with the station records of the corresponding dates. Similar to the calibration period, the calibrated parameter ranges provide a reasonable 95PPU band for the validation period, with the P-factor and R-factor values of 0.87 and 1.10, respectively (Abbaspour et al. 2015). In the case of the use of the best-performing parameter values, an acceptable simulation performance is achieved in terms of the R² (with b), NSE, PBIAS, and RSR statistics, similar to the calibration results, as detailed in Table S2 (Moriasi et al. 2007, 2015). The statistics and graphical representation of the simulation results over the validation period demonstrate that the calibrated SWAT model reasonably captures both the average and extreme daily flow rates at the Baykan SGS location, making it suitable for use in climate change analyses.

The bias-corrected GCM simulations are used as input in the calibrated SWAT model to project streamflow rates at the Baykan SGS location for the periods 1955–2010 and 2025–2099. Achieving an appropriate bias-adjusting performance guarantees reliable streamflow projections (Girvetz et al. 2013; Sun et al. 2020). In this study, the LS method is preferred to account for potential systematic errors in the raw GCM simulations due to its simplicity and accuracy (Salman et al. 2020). The study employs four error metrics, namely the modified index of agreement (md) (Legates & McCabe 1999), normalized root mean square error (nRMSE) (Almeida et al. 2015), Kling–Gupta efficiency metric (KGE) (Gupta et al. 2009), and fractions skill score (FSS) (Roberts & Lean 2008), to analyze the performance of both the raw and bias-corrected historical simulations of the 21 CMIP6 GCMs. Perfect agreement is characterized by an md value equal to 1, an nRMSE value equal to 0, a KGE value equal to 1, and an FSS value equal to 1.

Considering the ensemble medians attained under the CMIP6 historical experiment as the baseline scenario, the future changes in the climate of the Baykan SGS basin are analyzed on a seasonal and annual basis over the periods of 2025–2049 (near-future), 2050–2074 (mid-future), and 2075–2099 (long-future), using the median time series of the areal-averaged daily climate projections of the GCMs under the future scenarios of SSP245 and SSP585, as summarized in Table 2. Under the SSP245 scenario, while no substantial changes are projected for the mean annual daily precipitation amount in the future periods, the mean daily precipitation rate during the summer season decreases by 12.6 and 12.0% in the mid- and long-future periods, respectively. The projected gradual increases in the mean annual maximum and minimum air temperature rates throughout the future periods reach 4.3 and 4.0 °C in the long-future period, respectively. Under the SSP585 scenario, there are no significant differences in the mean annual daily precipitation amounts between the historical period and the near- and mid-future periods. However, the mean annual daily precipitation amount in the long-future period is 11.2% less than the precipitation rate attained under the CMIP6 historical experiment. Moreover, the projected decreases in summer precipitation during the mid- and long-future periods are more pronounced than under the SSP245 scenario. For the summer seasons of these periods, the percentages of precipitation decrease are determined to be 17.6 and 27.4%, respectively. Furthermore, the most substantial seasonal precipitation decreases are projected over the long-future period, with rates of 9.7, 8.5, 13.8, and 27.4% for the autumn, winter, spring, and summer seasons, respectively. Regarding the temperature anomalies, the gradual increases in the mean annual maximum and minimum air temperature rates reach 7.5 and 6.4 °C, respectively, in the long-future period.

Before conducting flood frequency analysis, the daily streamflow simulations of the calibrated SWAT model, attained by using the bias-corrected climate projections of the 21 CMIP6 GCMs under the CMIP6 historical experiment, are assessed against the Baykan SGS records in the 1955–2010 period, both seasonally and annually. Accordingly, while the mean annual observed streamflow rate is 18.58 m³/s, the ensemble median from the 21 SWAT model simulations is determined to be 17.17 m³/s. Regarding the seasonal streamflow averages, the Baykan SGS measurements show mean streamflow rates of 5.38, 15.19, 46.20, and 7.33 m³/s for the autumn, winter, spring, and summer seasons, respectively. Correspondingly, the ensemble median flow rates for these seasons are determined to be 5.96, 13.72, 42.52, and 6.28 m³/s. It is important to note that the calibration and validation of the developed SWAT model are conducted over a shorter timeframe (i.e., 1988–2009) compared to the historical analysis period (i.e., 1955–2010). Moreover, the calibrated SWAT model is executed with representative daily wind speed, relative humidity, and solar radiation values produced by the weather generator within SWAT based on the inputted long-term weather statistics. Despite these limitations, the agreement between the mean annual and seasonal average flow rates from the station records and ensemble medians emphasizes the robustness of the model calibration and the efficacy of employing a large ensemble to mitigate the uncertainty contribution of GCMs (Wang et al. 2020).

To comprehend the impacts of the projected climate changes across the study basin on the streamflow regime at the Baykan station location, the median time series of the daily streamflow rates simulated using the climate projections from the 21 GCMs under the future scenarios of SSP245 and SSP585 are scrutinized both seasonally and annually for the near-, mid-, and long-future periods of 2025–2049, 2050–2074, and 2075–2099, as outlined in Table 2. In this evaluation, the median streamflow rates achieved for the climate conditions of the CMIP6 historical experiment are considered as the baseline scenario for fair comparisons. Accordingly, under the SSP245 scenario, the mean annual streamflow rate is projected to decrease by 4.4, 8.0, and 5.0% in the near-, mid-, and long-future periods, respectively. Meanwhile, under the SSP585 scenario, the corresponding decrease percentages are obtained to be 4.9, 7.5, and 18.2%. Noteworthy alterations are not anticipated for the autumn seasons of the considered future periods in either climate scenario. Regarding the summer flows, the decrease percentages for the near-, mid-, and long-future periods are projected to be 14.6, 17.1, and 14.7% under the SSP245 scenario, and 13.5, 19.2, and 34.5% under the SSP585 scenario, respectively. Although the projected streamflow reductions for the autumn and summer seasons are relatively consistent with the projected changes in the precipitation amounts, it is discerned that the progressively increasing temperatures across the basin throughout the upcoming periods will shift the hydrological regime from a snow-toward a rain-dominated season. Under the SSP245 scenario, the mean streamflow rates of the winter seasons are projected to be 33.0, 32.0, and 61.7% higher than the historical winter average for the near-, mid-, and long-future periods, respectively. Conversely, reductions of 15.5, 21.2, and 25.2% are anticipated for the spring seasons in these respective periods. These changes are even more pronounced under the SSP585 scenario due to more significant temperature increases. Under the SSP585 scenario, the mean streamflow rates during the winter seasons for the near-, mid-, and long-future periods are projected, in turn, to be 33.8, 58.7, and 76.3% higher than the historical winter average, while reductions of 17.3, 27.2, and 47.4% are expected for the spring seasons of these periods.

The utilization of the best-fitted distribution functions' results for calculating median flood discharges is also examined for the historical analysis period of 1955–2010. The results of the K-S test indicate that among the considered distribution functions, the Pearson type 3 is the most appropriate one for the annual maximum daily streamflow records of the Baykan SGS. However, the most appropriate distribution function varies for each set of the annual daily streamflow peaks simulated under the historical GCM data, as depicted in Figure S1 (Supplementary Material). Although the results of different distribution functions are utilized in the calculation of the median flood discharges for the considered return periods, it is observed that the median flood discharges based on the best-fitted distribution functions are highly consistent with the flood discharges obtained for the station records using the Pearson type 3 function, as presented in Figure 5(g) and detailed in Table 3. While values of 294.5, 344.7, 401.9, 441.0, 477.7, and 556.3 m³/s are obtained for the 5-, 10-, 25-, 50-, 100-, and 500-year flood discharges using the observed streamflow peaks under the Pearson type 3 function, the median flood discharges attained under the use of the best-fitted distribution functions are determined to be 264.8, 320.9, 392.8, 441.9, 484.1, and 578.2 m³/s, respectively.

Considering the projected changes in precipitation, temperature, and streamflow patterns within the study basin, it can be inferred that the impacts of climate change will intensify over the long-future period (i.e., 2075–2099). Hence, flood frequency analysis for the upcoming decades is conducted by encompassing the periods 2025–2074 and 2025–2099 to account also for the influence of the selected timeframe when assessing potential responses of flood discharges to the changing climate. The flood discharges calculated using the six distribution functions for the annual daily peak flow rates simulated under the climate projections of each GCM for the SSP245 and SSP585 scenarios are included in Figure 5(a)–5(f), and the medians of these flood discharges are provided in Table 3. To enable equitable comparisons, the median flood discharges obtained under the climate conditions of the CMIP6 historical experiment are considered as the baseline scenario for quantifying changes in the flood discharges projected for the 2025–2074 and 2025–2099 periods. The results reveal an upward trend in all projected flood discharges, and the increases are even more evident in the flood discharges projected for the 2025–2099 period under the SSP585 scenario. However, the rise percentages for the median flood discharges differ significantly based on the distribution function used, especially for the discharge rates with a return period of 25 years or longer. For instance, while the increase percentages for the 10-year median flood discharge in the 2025–2099 period analysis are determined to be in the range of 20.0–22.3% under the SSP245 scenario and 23.9–26.8% under the SSP585 scenario, the ranges for the 100-year median flood discharge are found to be 14.2–23.2% and 23.5–31.5% for the SSP245 and SSP585 scenarios, respectively.

Due to the high level of consistency between the median discharge rates achieved using the results of the best-fitted distribution functions and the flood discharges of the most appropriate distribution function for the station records, projecting changes in flood discharges for future periods using the results of the best-fitted distribution functions can be considered a viable approach for use in hydraulic structure design. Figure S1 presents the best-fitted distribution function for each future streamflow peak time series, while Figure 5(g) and Table 3 include the projected flood discharges and their medians, respectively. Accordingly, under the SSP245 scenario, the 5-, 10-, 25-, 50-, 100-, and 500-year median flood discharges for the 2025–2074 period are 11.9, 13.0, 10.9, 5.7, 4.6, and 12.9% higher, respectively, compared to the median flood rates based on the streamflow peaks simulated under the climate conditions of the CMIP6 historical experiment. The corresponding percentage increases for the simulated streamflow peaks of the 2025–2099 period are 19.0, 20.0, 16.0, 10.7, 11.4, and 17.9%. Notably, the SSP585 scenario projects even greater increases. Under the SSP585 scenario, the percentage increases for the 5-, 10-, 25-, 50-, 100-, and 500-year events are 16.3, 11.5, 17.1, 21.7, 23.3, and 26.9% for the 2025–2074 period, and 21.1, 25.2, 25.3, 29.3, 27.6, and 31.7% for the 2025–2099 period. While it is advisable to utilize the median flood discharges projected for the 2025–2099 period under the SSP585 scenario in the design of hydraulic structures, other factors such as project significance, economic feasibility, and project lifespan will also influence the decision regarding which climate scenario or future period to consider for determining design discharges in response to the changing climate conditions in the upcoming years.

This study proposes an integrated framework to quantify potential changes in flood frequencies under the influence of future climate, with a focus on Bitlis Creek in the most threatened highlands region of the Euphrates-Tigris Basin. Climate projections from 21 GCMs belonging to the CMIP6 database, under both historical and future scenarios, are utilized to evaluate the increasing probability of floods due to climate-induced changes in the streamflow regime of the creek at the Baykan SGS location. The SWAT model is employed to simulate annual daily streamflow peaks based on the climate projections of each GCM. Five probability distribution functions are utilized to calculate the 5-, 10-, 25-, 50-, 100-, and 500-year flood discharges for the historical period from 1955 to 2010 and the two future periods from 2025 to 2074 and 2025 to 2099. To mitigate the uncertainty arising from differing GCM projections, the quantification of climate change impacts on the flood discharges is based on the medians of the discharge rates obtained under the climate data of each GCM. Given the variability in the projected median flood discharges based on the distribution function used, employing the best-fitted distribution functions is considered preferable when assessing the median flood discharges intended for design purposes. In this approach, the results indicate that the highest increases in the discharge rates for the 5-, 10-, 25-, 50-, 100-, and 500-year return periods are observed for the 2025–2099 period under the SSP585 scenario, with percentage increases of 21.1, 25.2, 25.3, 29.3, 27.6, and 31.7%, respectively.

The proposed methodology is expected to aid in determining the design of flood discharges for hydraulic structures planned in basins where climate change is foreseen to pose a significant threat to water resources in the coming decades (e.g., Cong et al. 2017; Shahid et al. 2018, 2021). For future studies, it is advisable to explore different bias correction techniques for the GCM datasets and to utilize higher-quality geospatial data (i.e., topography, land cover, and soil) in constructing the hydrological model to improve the consistency of the streamflow projections with the station records (Özcan et al. 2016; Wang et al. 2020). It is highly recommended that future assessments of climate change impacts rely on simulations from CMIP6-based regional climate models, which will become increasingly available in the coming years. Moreover, the proposed framework can be enhanced by considering future changes in other climatic variables (i.e., wind speed, relative humidity, and solar radiation) when simulating streamflow peaks (Gorguner & Kavvas 2020). While the proposed methodology allows for the assessment of how flood discharges for different return periods will change at the Baykan SGS location throughout the 21st century, it is recommended to conduct site-specific assessments for planning and designing water-related infrastructures within the basin. Additionally, for the dam projects on Bitlis Creek, a more comprehensive analysis incorporating future changes in the snowmelt and baseflow regimes should be carried out to provide reliable estimates of the maximum probable flood discharges necessary for designing spillway structures (Usul 2009).

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

The authors declare there is no conflict.