Hydrological model optimization using multi-gauge calibration (MGC) in a mountainous region

Hydrologic models have been widely used in hydrological research and water resource management studies. Hydrological models are increasingly used to simulate changes in the basin management, not only for investigating the impacts of external influences, but also for discovering the impacts of future data series (observing climate change, drought and flood studies) (Zhang et al. 2016). The physical hydrological model, namely Soil and Water Assessment Tool (SWAT), is a river basin scale model requiring specific information about weather, soil properties, topography, vegetation, and land management practices occurring in a basin (Arnold et al. 1998; Neitsch et al. 2011; Winchell et al. 2013). It has gained international acceptance as a robust interdisciplinary basin modeling tool. However, calibration has been one of the major challenges in developing hydrological models in such a way that simulations will be in good agreement with observed flows. The SWAT model calibration has been improving due to continuing advancements in the calibration techniques (Van Griensven et al. 2006; Arnold et al. 2012a, 2012b; Ercan et al. 2014).

In past decades, many calibration methods have been developed for SWAT, including manual calibration and automated procedures using the shuffled complex evolution (SCE) method and others. Arnold et al. (2012a, 2012b) developed the SWAT-CUP program that facilitates the calibration and validation of SWAT model using the semi-automated approach (SUFI2) (Abbaspour 2015). The sensitivity analysis of SWAT-CUP helps check model parameters that have a high impact on the simulations and make the necessary changes accordingly and model parameters must only be calibrated in a reasonable range that could give meaning on the ground (Van Griensven et al. 2006). For calibration studies, it is important to spatially account for hydrological processes in order to improve the efficiency of model run time and incorporate the impact of uncertainty in the model (Klemes 1986). This task is one of the major challenges that enhances SWAT model calibration. Therefore, the aim of this study is to investigate the impact of incorporating hydrological processes spatially in order to improve model calibration. In addition, the impact of the model calibration techniques on the hydrological modeling of the basin and climate change impact assessment will be studied.

A number of attempts have been conducted by hydrologists to incorporate model parameters spatially to improve model accuracy. In 2005, a multi-variable and multi-site approach to the calibration and validation of the SWAT model for the Motueka catchment, New Zealand, was conducted and found that the method improved the model calibration and validation (Nash–Sutcliffe coefficient of 0.78 and 0.72 respectively) (Cao et al. 2005). In another study, spatially distributed calibration at sub-basin level and temporal validation at the stream gauges outlet points was conducted to incorporate the spatial and temporal hydrological patterns for two river basins in the Ohio, USA, region and found out that spatially distributed calibration and validation of the basins improve the model predictive capacity (Santhi et al. 2008). Similar studies conducted by Lu et al. (2015), Bai et al. (2017) and Zhang et al. (2017) in the Yingluoxia, Miyun and Baihe basins, China, comparing single- and multi-site calibration and validation of the SWAT model, showed that the parameters determined from multi-site calibration and validation were far better than those from the single site calibrations. Haas et al. (2016) studied a joint multi-metric calibration of flow simulated by SWAT in Treene basin, Germany. The results showed that the approach used has improved the calibrated predictions, where adequate model runs with good performance for different hydrological conditions for the flow were detected. Contrary to the above, Shresthaa et al. (2016) indicated that multi-site calibration did not improve simulations of flow and sediments compared to single-site calibration.

Changes have occurred in several aspects of the atmosphere and surface that alter the global energy budget of the Earth causing the climate to change (Solomon et al. 2007). The Intergovernmental Panel on Climate Change (IPCC), set up by the World Meteorological Organization (WMO) to provide information on climate change, has confirmed with high confidence that the gradual increase in the average temperature of the Earth's surface is due to the growing concentration of greenhouse gases (GHGs) in the atmosphere, attributed to burning of fossil fuels and changes in land use and land cover (IPCC 2001). Today, Regional Circulation Models (RCMs) are applied to determine the likely effects of these changes in the GHGs. It is to be mentioned here that a GCM/RCM is a Global/Regional Circulation Model that could depict different pictures for the future weather over a region, and consequently affect the hydrologic regime of basins (Elshamy et al. 2009). Increases have occurred in the number of heavy precipitation events. The availability of observational data restricts the types of extremes that can be analyzed. Confidence in these estimates is higher for some climate variables (e.g. temperature) than for others (e.g. precipitation).

Some studies have been conducted to investigate the effects of multiple stream gauging stations' data on basin-wide calibration and validation in Turkey. The impact of size of sub-basins on the hydrologic parameters and their spatial variability in the estimation of the hydrologic parameters in the western Black Sea region were studied by Kocyigit et al. (2017) and Akay et al. (2018). The studies indicated that multi-site calibration and validation resulted in satisfactory outcomes for the direct flow hydrograph but not the peak flow prediction. These studies have only focused on the application of hydrological models, and very limited attempts were made in the model parametrization and calibration technique studies for basins in Turkey. Therefore, there is an obvious need to assess the impact of spacial dimension and calibration and validation techniques in hydrological modeling.

Consequently, the main objective of this study is to evaluate the impact of multi-site calibration on the performance of hydrological model calibration and validation. To achieve this, we used four different calibration techniques: (1) calibration using gauge data only from the downstream area (DC); (2) calibration using gauge data only from the upstream area (UC); (3) calibration using gauge data from both the upstream and downstream areas of the basin (Multi Gauge Calibration – MGC); and (4) calibration using gauge data first from the upstream area only, then using the data from the downstream area sequentially (UCDC). We analyzed the impact of these techniques on the hydrologic regime of the basin. Hydrological model simulations and statistical analysis, as well as climate change impact assessments, have been used in the study.

The Rize province is located in the northeastern Anatolian mainland of Turkey, so-called the eastern Black Sea region, with an area of 3,920 km². The average annual precipitation in the region is nearly 2,250 mm and total annual runoff is about 2,745 million m³ based on the reference period 1960–1996 (Sen & Kahya 2017). The study area is a mountainous region with elevations reaching higher than 3,000 m (Figure 1) and a very wet climate (Kahya et al. 2008a). In addition, Rize province falls into the homogeneous streamflow region (Kahya et al. 2008b).

The Digital Elevation Model (DEM), the soil map, and the land use map of the basin were obtained from Istanbul Technical University (ITÜ) – Hydraulics and Water Resource Engineering Laboratory. The digital elevation model (DEM) with 28-meter spatial resolution was used which is quite sufficient for hydrological modelling studies. A total of 22 sub-basins were identified after delimiting the basin boundary in ArcSWAT (Figure 1).

The observed precipitation data was accessed from the Turkish Meteorological Directorate (MGM). We used the maximum and minimum temperature data from the Climate Forecast System Reanalysis (CFSR) which was proven to have good accuracy with the observed temperature (Duana et al. 2019). The precipitation data in the vicinity of the Ikizdere basin was used for the SWAT simulations. In addition, the daily flow data of three flow gauging stations found in Ikizdere basin were received from the Turkish General Directorate of State Hydraulic Works (DSI) (for the years 1971–1996) to carry out the model calibration and validation. The flow gauging stations of Camlıkderesi or Ikizdere (station no. E22A015) and Tözköyderesi (station no. E22A033) were used to assess the impacts of spatial streamflow data on the overall model optimization performance. The third gauge of Cimil-Koknar (station no. D22A096) in Ikizdere basin was insufficiently small to be used for the study. After failing to simulate the Kaptanpaşa basin of Rize province within acceptable accuracy, due to a shortage of observed climate and streamflow data, the Ikizdere basin was selected for our study as sufficient data is available for hydrological modelling, which is also recommended by Kahya et al. (2015). Moreover, the mean monthly values of temperature, wind speed, solar radiation and humidity were generated from the Climate Forecast System Reanalysis (CFSR) data, which can be accessed from the website of US National Centers for Environmental Prediction (NCEP-CFSR). The climate change data for the three climate models and two scenarios were downloaded from the official IPCC data repository (www.ipcc-data.org/sim/gcm_clim/SRES_AR4/index.html)

SWAT is a river basin scale model developed to predict the impact of basin management on water, sediment, nutrients, and agricultural and chemical yields (Arnold et al. 1998). In this study, the ArcSWAT program was used in setting up the model for the study area. The observed precipitation data from the Turkish Meteorological agency stations was employed in the SWAT model calibration. The model uses only one climate station in each sub-basin when a single station is available within the sub-basin; otherwise, it demands the mean data of nearby stations. Once the delineation of the basin is completed, Hydrological Response Units (HRU) were categorized with respect to the soil, land use and slope maps overlay. Categorizing all the similar parcels of land with similar hydrology, SWAT employs a lumped calculation for parameters instead of a distributed calculation for the sake of saving time and resources. Once the HRU and sub-basin analysis is completed, the ArcSWAT interface then inputs climate data (such as precipitation and temperature) to finalize the SWAT model setup. The two streamflow series in Camlıkdere station located at the downstream and Tözköydere station located at the upstream of the study area were used for the model calibration and validation.

The SWAT model setup based on the study area was simulated using precipitation data for the period 1971–1996. The observed precipitation was assumed the only dominant climate variable that has direct and significant impacts on streamflow downstream. The simulations were run on a monthly time step and a warm-up period of three years for the model to adjust itself according to the data and basin conditions. The output of the model in the warm-up period was not used for further analysis. Following successful setting procedures for SWAT, it was calibrated and validated using the remaining historical streamflow data. More specifically, the first three years (1971–1973) were used in the model warm-up, and the next 12 years (1974–1990) in the model calibration, and finally the remaining six years (1991–1996) for the model validation. The model was calibrated using the following four methods: the first one using Ikizdere streamflow data (Downstream Calibration – DC), the second using Tözköydere streamflow data (Upstream Calibration – UC), the third using both Ikizdere and Tözköydere streamflow data (Multi-Gauge Calibration – MGC), and finally the fourth calibration using Tözköydere streamflow before the calibration using Ikizdere streamflow data (Upstream Calibration before Downstream Calibration – UCDC). From the review of relevant literatures, the technique employed in this research is quite different from what has been done in the past (at least for the Black Sea region). Thus, it will have a significant contribution in highlighting the significance of incorporating spatial dimension in hydrological model parametrization studies.

In the phase of model calibration and sensitivity analysis, we adopted the SWAT-CUP auto-calibration tool (Abbaspour 2015). The Latin Hypercube One-factor-At-a-Time design (LH-OAT) method proposed by Morris (1991) was implemented to carry out the sensitivity analysis method which uses a stratified sampling approach that better covers the sampling hypercube with fewer samples (White & Chaubey 2005; Van Griensven et al. 2006). The results are the relative sensitivities based on the linear approximations, providing only partial information about the sensitivity of objective function to the model parameters. In general, the larger the absolute value of t-stat and the smaller the p-value, the more sensitive is the parameter. The most sensitive model parameters are essential in model calibration as these parameters will affect the model more than less sensitive parameters. The four model calibration methods were tested for performance. In each method, a set of five calibration steps were conducted until the simulated streamflow at a desirable performance is attained as compared to the observed streamflow. It is important to note that each calibration step consists of 500 iteration runs. The final step was validation to evaluate the performance of the calibrated model with a separate precipitation forcing.

In this study, three regionally downscaled GCM climate models (MIMR, INCM3 and HADCM3) and two scenarios (A2 and A1B) were used. It is evident that the three climate models are not consistent in their estimation of future precipitation change. Except for the A1B scenario of INCM3, all the other climate change scenarios predicted the same or decreased precipitation across the seasons. The summary of IPCC climate change data used in our study is presented in Figure 2. The values indicate the projected precipitation in mm at the end of the century (2080–2099) with respect to the baseline period (1961–1990). We implemented the climate change analysis by directly substituting the percentage change in precipitation forecast by the climate model for each month of the year in the SWAT input/output files (TxtInOut folder). We copied the ‘swat.exe’ executable file to run the SWAT model from the TxtInOut folder after making the necessary change in all the sub-basin (.sub) files.

In this study, we used the SWAT model using the SWAT-CUP model optimization tool to assess our four model optimization techniques. In addition, we assessed the impact of climate change for the study area.

The SWAT model of the study area was set-up and calibrated using the observed flow data. The output of the model calibration, fitted parameter values and ranks of the sensitivity analysis are illustrated in Supplementary Material, Annex 1. Some of the parameters are significantly affected by the streamflow stations used in the model calibration. CH_N2 and CN parameters have shown to be sensitive, which indicated the direct impact of river channel characteristics on streamflow. In addition, SFTMP and SNO50COV were found to be sensitive which indicates the effect temperature on the flow regime of the basin. The majority of model parameters are less sensitive when it comes to the location of gauge data used for calibration.

Irrespective of the four model calibration techniques used, the best parameter values fitted for each calibration technique are quite similar along the calibration techniques used. The only exceptions are: CN2, CH_K2, LAT_TTIME, ALPHA_BF, GWQMN, PLAPS, SLSOIL and SOL_K. Most of these parameters are characterized by their effect on the river channel and soil characteristics.

Five statistical parameters: NS, R2, RSR, PBIAS and p-factor, were calculated to evaluate the performance of the model and the calibration techniques (Table 1). According to Moriasi et al. (2007), model simulation can be judged as satisfactory if NS > 0.50 and RSR < 0.70, and if PBIAS < ±25% for streamflow analysis. The SWAT model was calibrated using four methods (DC, UC, MGC, and UCDC). The calibration with Tözköydere flow data (UC) gave a good NS parameter (0.8 and 0.7 for calibration and validation respectively). That is due to low flow intensity upstream of the basin caused by a smaller contribution area for the river outlet. This makes calibration easier for the objective functions to achieve a large number in a short period of time. The calibration with only Ikizdere flow (DC) gave the NS to be 0.7 both for calibration and validation. The DC calibration technique resulted in a smaller NS value due to a larger contribution area for flow causing high flow at the outlet. Maximizing the objective function (NS) for a larger basin is more difficult than a basin of small area. The third technique (MGC) used both Ikizdere and Tözköydere flow for the model calibration. The objective function for both gauging stations were lower than that achieved in the previous steps (0.7 and 0.6 for Ikizdere calibration and validation respectively, 0.8 and 0.6 for Tözköydere calibration and validation respectively). The fourth method calibrated the model using Tözköydere flow before calibrating with Ikizdere flow values. The final model parameters best fits calibrated using the Tözköydere flows were not changed when calibrating using Ikizdere flows. That way, the interdependence of the upstream and downstream flows was maintained. Although the objective function was lower for the third technique (MGC), the mean, standard deviation and PBIAS values were found to perform better than the first two calibration techniques. The best objective function performance was observed for the fourth technique (UCDC), which gave far better NS values both for the upstream and downstream regions of the basin. The other observation was calibrating the model using the upstream gauge data (UC) which gave more accurate model performance when compared to the calibration using the downstream flow data (DC). The calibration of the model with the flow data of Ikizdere resulted in improved NS values at the Tözköydere location too. However, the reverse did not show a similar outcome (Table 1).

In the temporal analysis, the monthly average simulated flow values are plotted against time. The calibration of the model using the upstream flow (UC) (NS of 0.81 for calibration) compared to the calibration using the downstream flow (DC) (NS of 0.72 for calibration) demonstrated good accuracy (Figures 3 and 4). The main reason is that the upstream flow intensity is lower, making it less difficult for the algorithm to capture the distribution including the peak flow values. When the flow value increases downstream, it becomes difficult to estimate the flows, especially the peak flows resulting in lower calibration performance. When both the upstream and downstream flow data are used for model calibration simultaneously (MGC), the algorithm adjusts all model parameters providing the best optimization value (NS in our case) (Figure 5).

In addition, the temporal analysis proved the superiority of the fourth technique (UCDC) as compared to the first three techniques, since the hydrological model was optimized using two flow gauging stations one after the other (Figure 6). This allows the model to be optimized sequentially, starting from the upstream region down to the downstream region. The other techniques optimize the model with local gauge data contrary to the holistic approach implemented in UCDC. This results in loss of information about the characteristics of the basin which could otherwise be useful to improve the model performance. In the UCDS technique, we made sure not to change the parameters that were set using the upstream flow when calibrating using the downstream flow in order to maintain the inter-dependence of upstream–downstream flow.

The monthly average flow values for the simulated flow, as well as six climate change scenarios compared with the observed flow at Ikizdere sub-basin outlet, is plotted in Figure 7. There is no consistent trend in the forecasted flow for the various climate model precipitations. For the winter, all the climate models forecast a decrease in flow for the river. The other seasons forecasted a slight increase in flow. All in all, we did not observe an agreement amongst the climate projections forecast for the Rize basin.

It is quite challenging for a hydrological model to capture the complex processes in a basin, especially for mountainous regions with limited input data. Estimating the parameters that represent the various hydrological processes is one of the main issues hydrologists need to solve since it is impossible to measure all of the hydrological parameters. In order to assess the impact of calibration techniques on the model performance, four methods were tested, calibration using flow data: (1) only from upstream area (US), (2) only from downstream area (DC), (3) using both upstream and downstream areas of the basin–multi-gauge calibration (MGC), and (4) first the upstream then downstream flow data (UCDC). The study was aimed at assessing the performance of multi-gauge calibration over single gauge calibration techniques.

The results of our analysis have shown that using individual gauges for calibration improve the predictive ability of the model for the particular region. The DC calibration technique resulted in a lower performance compared to the UC, due to a larger drainage area for the outlet causing high flow. In calibration, maximizing the objective function (NS) for a large basin area is more difficult than for a basin of small area. Model calibration using the MGC technique significantly improved model performance as a whole, unlike the single gauge calibrations. This is due to the fact that the model obtains much wider information on the basin characteristics to adjust its parameters when more than one gauging station is used simultaneously for calibration. The UCDC calibration technique gave the best model performance since the model performance was improved step by step starting from the upstream area going to the downstream area of the basin. This technique not only improved the model for each region, but also the interdependence of the two regions was maintained since the model parameters calibrated by the upstream flow were not changed when calibrating the model with the downstream flow. Our study findings are in agreement with what has been reported by Cao et al. (2005), Santhi et al. (2008), Haas et al. (2016), Lu et al. (2015), Zhang et al. (2017), Bai et al. (2017) and Akay et al. (2018). They all reported improvements in model calibration performance with multi-site (multi-gauge) calibration compared with single-gauge calibrations for their studies. Evaluating the results of previous studies and our study, we can confidently say multi-gauge calibration gives much better model calibration performance than single-gauge calibration techniques for the Ikizdere basin, Black Sea region. However, we still recommend more studies to be carried out on other basins to reach a comprehensive conclusion about multi-gauge calibration techniques for the Black Sea region. Regarding the climate change assessment, we have not observed a consistent trend in the forecasted flow by the various climate model predictions. The IPCC climate projections forecast did not show a visible shift in the hydrologic regime for the Rize basin.

We would like to acknowledge TUBİTAK (the Turkish Scientific and Technological Research Council) for providing the funds to carry out this study. We are also indebted to the Turkish State Meteorological Service (MGM) and General Directorate of Hydraulic Works (DSI) for the observed weather and flow data access.

Data cannot be made publicly available; readers should contact the corresponding author for details.