*2.2. Hydrological Model*

*2.2. Hydrological Model* The Hydrologic Engineering Center's Hydrologic Modeling System (HEC‐HMS) is a hydrological model developed to simulate the precipitation–runoff processes of dendritic drainage basin systems [18]. The model is designed for both continuous and event‐based hydrological modeling and offers several different options for modeling the various com‐ ponents of the hydrological cycle. In event‐based modeling, storm precipitation is simu‐ lated during the simulation time interval ranging from a few hours to several days, de‐ pending on the basin size [45]. In continuous modeling, a continuous historical record of hydrological events, including dry and wet periods over several years, is simulated [46]. The main difference is that evapotranspiration and groundwater seepage can be neglected in event‐based modeling, while they cannot be ignored in continuous simulation [47]. HEC‐HMS can conduct hydrological simulation over a wide range with various simple modules to represent different components of the hydrological cycle. The selection of the appropriate model for each component depends on the experience of the modeler, the pur‐ The Hydrologic Engineering Center's Hydrologic Modeling System (HEC-HMS) is a hydrological model developed to simulate the precipitation–runoff processes of dendritic drainage basin systems [18]. The model is designed for both continuous and event-based hydrological modeling and offers several different options for modeling the various components of the hydrological cycle. In event-based modeling, storm precipitation is simulated during the simulation time interval ranging from a few hours to several days, depending on the basin size [45]. In continuous modeling, a continuous historical record of hydrological events, including dry and wet periods over several years, is simulated [46]. The main difference is that evapotranspiration and groundwater seepage can be neglected in eventbased modeling, while they cannot be ignored in continuous simulation [47]. HEC-HMS can conduct hydrological simulation over a wide range with various simple modules to represent different components of the hydrological cycle. The selection of the appropriate model for each component depends on the experience of the modeler, the purpose of the modeling, and the usability of the input data [48].

pose of the modeling, and the usability of the input data [48]. The HEC‐HMS modeling system has three main components: the basin model in which the topographic and physical characteristics of the basin are determined, the mete‐ The HEC-HMS modeling system has three main components: the basin model in which the topographic and physical characteristics of the basin are determined, the meteorological model in which the meteorological data are processed, and the control manager.

#### orological model in which the meteorological data are processed, and the control manager. *2.3. Basin Model*

*2.3. Basin Model* The HEC‐HMS basin model (Figure 3) simulates the process of the water falling to Earth by precipitation from the canopy to become groundwater, excluding bottom perco‐ lation. HEC‐HMS uses the soil moisture accounting (SMA) [31] algorithm to simulate the movement of water in soil under continuous simulations. This algorithm takes precipita‐ The HEC-HMS basin model (Figure 3) simulates the process of the water falling to Earth by precipitation from the canopy to become groundwater, excluding bottom percolation. HEC-HMS uses the soil moisture accounting (SMA) [31] algorithm to simulate the movement of water in soil under continuous simulations. This algorithm takes precipitation and evapotranspiration as inputs and computes surface runoff, groundwater runoff,

tion and evapotranspiration as inputs and computes surface runoff, groundwater runoff, evapotranspiration, and losses from bottom percolation (Figure 4; see USACE (2016) for

hydrograph simulation) and the monthly constant baseflow method was chosen for the baseflow calculation [49,50]. Initial parameters for the Clark method were obtained using the Kerby equation (Tc = G(L\*r/S0.5)0.467). The physical parameters of the sub‐basins at the exit of the three selected meteorological stations (for example, river length, drainage area, slope, etc.) were computed with Geographic Information Systems (GIS) by using the digital elevation model (DEM) obtained from the 10 m contour maps (Table 2). The initial values for the baseflow were taken as the current river flows because the beginning of the simulation was in the dry period, and they were distributed on the basis of the average

area‐based distribution in the sub‐basins.

evapotranspiration, and losses from bottom percolation (Figure 4; see USACE (2016) for further detail). The Clark Unit Hydrograph was chosen for the transformation method (or hydrograph simulation) and the monthly constant baseflow method was chosen for the baseflow calculation [49,50]. Initial parameters for the Clark method were obtained using the Kerby equation (Tc = G(L\*r/S0.5)0.467). The physical parameters of the sub-basins at the exit of the three selected meteorological stations (for example, river length, drainage area, slope, etc.) were computed with Geographic Information Systems (GIS) by using the digital elevation model (DEM) obtained from the 10 m contour maps (Table 2). The initial values for the baseflow were taken as the current river flows because the beginning of the simulation was in the dry period, and they were distributed on the basis of the average area-based distribution in the sub-basins. *Water* **2022**, *14*, x FOR PEER REVIEW 6 of 23

**Figure 3.** Conceptual model of the Kırkgöze–Çipak Basin in HMS, showing junctions, reaches, flow direction, and sub‐basins. **Figure 3.** Conceptual model of the Kırkgöze–Çipak Basin in HMS, showing junctions, reaches, flow direction, and sub-basins.

**Table 2.** Physical properties of the main sub‐basins. **MS1 MS2 MS3** Basin Slope (%) 0.173 0.167 0.215 Elevation (m) 2125.02 2175.37 2204.53 River Length (m) 21,917.56 12,094.25 11,208.97 Area (km2) 91.53 75.50 74.51 The basin model includes many parameters used for baseflow, hydrograph simula‐ tion, and SMA. For the estimation of these parameters in previous studies [10,45,51–53], it was found appropriate to use geodatabases, reducing the number of free parameters by starting the simulation during periods when initial conditions are easier to predict (i.e., the start of the water year), and use empirical equations or reliable sources. A combination of these methods was used in this study. As there was no map from which soil texture infor‐ mation of the study area could be obtained, the initial values of SMA were obtained from The basin model includes many parameters used for baseflow, hydrograph simulation, and SMA. For the estimation of these parameters in previous studies [10,45,51–53], it was found appropriate to use geodatabases, reducing the number of free parameters by starting the simulation during periods when initial conditions are easier to predict (i.e., the start of the water year), and use empirical equations or reliable sources. A combination of these methods was used in this study. As there was no map from which soil texture information of the study area could be obtained, the initial values of SMA were obtained from previous studies and then calibrated to match the observed streamflow. Canopy maximum retention and soil surface deposition were estimated by vegetation type and percentage of land slope, respectively [54,55]. The rate and amount of seepage in the soil profile and groundwater were estimated based on hydraulic conductivity [52]. Active soil depth was assumed to be 60 cm, considering the land cover. Fleming and Neary (2004) predicted HEC-HMS groundwater storage (groundwater 1 and groundwater 2), and seepage parameters [18] were based on recession analysis. These estimates from published literature were taken as initial values and they were calibrated during the simulations.

#### previous studies and then calibrated to match the observed streamflow. Canopy maximum *2.4. Meteorological Model*

retention and soil surface deposition were estimated by vegetation type and percentage of land slope, respectively [54,55]. The rate and amount of seepage in the soil profile and groundwater were estimated based on hydraulic conductivity [52]. Active soil depth was The Kırkgöze–Çipak basin is divided into a few sub-basins, as shown in Figure 3. The data obtained from three automatic meteorology and snow observation stations in

assumed to be 60 cm, considering the land cover. Fleming and Neary (2004) predicted HEC‐HMS groundwater storage (groundwater 1 and groundwater 2), and seepage param‐

taken as initial values and they were calibrated during the simulations.

the basin at altitudes of 2019 m (Kö¸sk), 2454 m (Güngörmez), and 2891 m (Radar) were used for the meteorological data required for the different parameter methods selected in the basin and for the meteorological model simulating the precipitation–runoff process. These stations provided time series of the maximum wind speed (m/s), wind direction, average air temperature (◦C), average humidity (% rh), air pressure (mbar), average soil temperature (◦C), solar radiation (W/m<sup>2</sup> ), average albedo, precipitation (mm), snow height (cm), snow density (gr/cm<sup>3</sup> ), and snow–water equivalent (cm) parameters over 15 min periods. References [35,36] showed that the climate data from the stations in the basin was sufficient, of good quality, and could be collected in real time. Measurements from the years 2008 to 2011 obtained from the Kö¸sk, Güngörmez, and Radar meteorology stations were used for the hydrological simulations to be conducted with HEC-HMS. *Water* **2022**, *14*, x FOR PEER REVIEW 7 of 23

**Figure 4.** Schematic of basin model in HEC‐HMS and its principal components. **Figure 4.** Schematic of basin model in HEC-HMS and its principal components.

used for the hydrological simulations to be conducted with HEC‐HMS.

The Kırkgöze–Çipak basin is divided into a few sub‐basins, as shown in Figure 3. The data obtained from three automatic meteorology and snow observation stations in the basin at altitudes of 2019 m (Köşk), 2454 m (Güngörmez), and 2891 m (Radar) were used for the meteorological data required for the different parameter methods selected in the basin and forthe meteorological model simulating the precipitation–runoff process. These stations provided time series of the maximum wind speed (m/s), wind direction, average air temperature (°C), average humidity (% rh), air pressure (mbar), average soil tempera‐

snow density (gr/cm3), and snow–water equivalent (cm) parameters over 15 min periods. References [35,36] showed that the climate data from the stations in the basin was suffi‐ cient, of good quality, and could be collected in real time. Measurements from the years 2008 to 2011 obtained from the Köşk, Güngörmez, and Radar meteorology stations were

*2.4. Meteorological Model*

*2.5. Precipitation Model*


**Table 2.** Physical properties of the main sub-basins.

#### *2.5. Precipitation Model*

A variety of different statistical techniques to distribute point observations over complex topography is given in published literature [56–64]. Although these studies improved high-resolution grid-type climate data estimations, uncertainties remained. In particular, it is more difficult to estimate the spatial distribution and the intensity of precipitation compared to other variables such as temperature, due to the regional, seasonal, and topographic characteristics [65]. The Kırkgöze–Çipak basin study area has a very large altitude range and other variable aspects, even though it is small in terms of scale. As a result of observations in the basin over a long time, it was determined that some convective precipitations took place independently from each other as in the northern aspects with quite high land altitudes where the Radar station is located and in the southern aspects where the Güngörmez station is located. Therefore, while snowmelt runoff simulations are performed throughout the basin, the emphasis is on how the precipitation is distributed regionally rather than how the precipitation may be distributed in the basin. The HEC-HMS program offers grid-based and polygonal-based solution alternatives to determine the precipitation distribution over the basin. This study was carried out on a polygonal basis, and the gage weights method was chosen for modeling the precipitation processes. The gage weights method is based on the Thiessen polygon method. The Thiessen polygon method, which is usually recommended for use in vast areas, does not distribute precipitation with respect to topographical effects and precipitation characteristics; instead, it performs it only over polygonal areas determined by the positions of the stations [66]. Therefore, the gage weights method used in HEC-HMS was modified for the study basin, which is heterogeneous in terms of altitude and exposure. While developing this polygonal area-based algorithm, in addition to the general behavior that is dependent on the topography of the region—the barrier effect (Figure 5, 4th elevation zone), the measurement data at the stream gauge stations and the ambient temperature, humidity, atmospheric pressure, rate of increase in cloudiness (observed), albedo, wind speed, and SWE values were all simultaneously examined. The area-based distribution of precipitation was simulated by six different zonal polygons shown in Figure 5.

In basins where there is a large altitude range, the use of data obtained from stations representing low levels may cause the precipitation input calculated for the entire basin to be lower than the actual value. It is recommended to extrapolate precipitation data to average hypsometric elevations for zones with elevation gradients [67], so that the point-based input values used in the modeling procedure can better represent a specific area. It is important for the precipitation data to align with other meteorological data with respect to time, so that the model can perform the necessary iterations accurately and reliably. For this reason, while making the area-based distribution of meteorological data, a general grouping based on altitude and exposure, taking into account station locations, was deemed appropriate so that simultaneous atmospheric homogeneity could be assured. For this reason, the meteorological variables in altitude zones 1 and 2 were based on the Güngörmez station, the meteorological variables in altitude zones 3 and 4 were based on the Kö¸sk station, and the meteorological variables in altitude zones 5 and 6 were based on the Radar station variables (Figure 5). In Table 3, the hypsometric elevations for each zone and the altitudes of the meteorological stations in these zones are given.

**Figure 5.** (**a**) Land survey of the elevation zones and (**b**) polygonal representation of the Kırkgöze– Çipak basin. **Figure 5.** (**a**) Land survey of the elevation zones and (**b**) polygonal representation of the Kırkgöze– Çipak basin.


As the hypsometric averages of zones 2, 3, and 6 are close to the altitudes of Gün‐ **Table 3.** Hypsometric values for the Kırkgöze–Çipak basin.

sition between stations exceeds the simulation time interval of 1 h, the precipitation is only distributed on the station's zone. Therefore, the predictive values may take zero values mathematically on the transition zones noticed on field trips. The daily total precipitation values in zones 1, 4, and 5 were analyzed according to the flow chart in Table 4, which was prepared by considering the station locations given in Figure 5 and calculated within the As the hypsometric averages of zones 2, 3, and 6 are close to the altitudes of Güngörmez, Kö¸sk, and the Radar stations, respectively, the average area-based precipitation calculations in these zones were taken directly at the respective stations. Based on the hypsometric elevations in other zones, a series of algorithms were run to obtain values in the direction of increasing or decreasing precipitation.

designated rules. Firstly, the 15 min precipitation series recorded at each station were converted into daily total precipitation series while removing possible measurement errors. This daily sum is a precaution for the following algorithm (Table 4), especially for modeling the natural distribution of the interpolated or extrapolated zonal values of the convective characteristic heavy snowfalls observed at the station points. Otherwise, if the precipitation transition between stations exceeds the simulation time interval of 1 h, the precipitation is only distributed on the station's zone. Therefore, the predictive values may take zero values mathematically on the transition zones noticed on field trips. The daily total precipitation values in zones 1, 4, and 5 were analyzed according to the flow chart in Table 4, which was prepared by considering the station locations given in Figure 5 and calculated within the designated rules.

> After calculating the daily precipitation altitudes for zones 1, 4, and 5, these altitudes were proportioned to the daily total precipitation altitude of zones 2, 3, and 6, respec-

tively, and precipitation coefficients were obtained. These coefficients were used for the conversion from the daily total precipitation altitude to the 15 min time interval. The first precipitation coefficient calculated to obtain the 15 min precipitation values for the 1st zone was multiplied by the 15 min precipitation series of the 2nd zone. This procedure was also performed for zones 4 and 3 and zones 5 and 6, respectively. Thus, while maintaining the atmospheric homogeneity, the precipitation altitudes and timings of the 1st, 4th, and 5th zones were adjusted with reference to the measurements taken from the 2nd, 3rd, and 5th zones, respectively. The results obtained were increased to a one-hour time interval selected as the HEC-HMS simulation time interval and entered into the program.

Complete reliable data could not be obtained from the pluviographs during the winter months because the diluted antifreeze in the rain gauge froze after a certain period of time under the effect of the cold weather and excessive precipitation at both the Radar station and the Güngörmez station; the movable scale shaft which measures the amount of precipitation discharged from the reservoir did not work due to freezing and jamming, even when the antifreeze did not freeze. The data from the Kö¸sk station showed that the freezing did not occur there due to the fact that the temperature was relatively higher than the other stations due to its lower altitude. As a result, much more reliable precipitation data were obtained there compared to the other stations during the winter months. Due to the problems encountered, especially at the Radar and Güngörmez stations, it was not found appropriate to use the data obtained from the rain gauges as direct precipitation data.

**Table 4.** The algorithm used for determining precipitation altitudes in the elevation zones where there was no station.


In winter, while the precipitation series were formed during the snow accumulation period, the differences in the 24-h averages of the snow–water equivalent altitudes (SWE) obtained from the snow pillows were taken. If the difference between these daily averages was positive, the SWE difference for that day was added to the station as precipitation. The timing of precipitation was adjusted in correlation with simultaneous albedo and humidity data, taking into account the effect of snow drift, while the distribution of precipitation during the day was determined by the amount of increase in the measured SWE during the day.

#### *2.6. Snowmelt Model*

In HEC-HMS (Version 4.2.1), there are two snowmelt modeling options. One of them is the gridded temperature index method and the second is the temperature index method, which is the method that was used in this study. The temperature index method is an extension of the degree-day approach to modeling a snowpack. A typical approach to the degree day is to have a fixed amount of snowmelt for each degree above freezing. This method includes a conceptual representation of the cold energy stored in the pack along with a limited memory of past conditions and other factors to compute the amount of melt for each degree above freezing. As the snowpack internal conditions and atmospheric conditions change, the melt coefficient also changes [18].

If the main source of energy in the spring is not the solar irradiance, snowmelt can be more effectively and simply computed using a temperature index model [10,68–73]. In hydrologics, an index is a meteorological or hydrological variable. Changes in the variable are associated with those in the parameter it is estimating, and which are more easily measured than the actual parameter. Either a coefficient (such as a degree-day factor) or a formula for more complex linear or curvilinear functions (such as the antecedent temperature index—melt rate function) may be used to describe this index relationship. Depending upon changing associated factors, it may be either constant or variable. Spatial and temporal basin value point measurements are represented by the index where average fixed relationships are known to exist between the measured values and basin values. However, snow accumulation and melting topics are complex, and the data required for physically-based energy budget calculations are comprehensive and challenging to obtain [68].

Some temperature index models require the snowpack's melt rate to be characterized [74,75]. This melt rate can be stated differently. One example is to express changes in the melting rate as a function of air temperature accumulation over several warm days for melting snow. This is achieved by using the ATIMR (antecedent temperature index—melt rate) function to determine the melt rate for a certain antecedent temperature index. Snow physics indicate that melting rates increase throughout the season due to both metamorphic processes causing ice crystal consolidation and the snowpack producing more water over time [24].

Past modeling studies have generally been based on a theoretical constant ATIMR curve generated by the USACE (1991) and used for characterization of melt rates [19–22,76,77]. The theoretical curve was included in the SSARR model in 1991. The SSARR guide [78] ATIMR values are shown in Table 5.


**Table 5.** Tabulation of melt rate as a function of ATIMR.

For the values in the customary U.S. system, please see the SSARR model guide in Appendix D, p. 17; the methodology used to calculate the metric system results is presented by ¸Sengül and ˙Ispirli (2021) in detail.

Modified melting rates have been used by many studies using the hypothetical ATIMR function of Table 5 during calibration; however, a commonly observed shortcoming in published literature is that no particular data is used to directly estimate the ATIMR curve. Sometimes the hypothetical ATIMR curve is taken as a starting point for snowpack simulations and different scenarios used to modify the curve to improve simulated results during calibration [19–22]. Sometimes the theoretical ATIMR curve is not modified, but an additional rate is applied to the melting rate obtained from the ATIMR curve in proportion to the varying rate over time [77,78]. However, the physical meaning of widely used ATIMR functions is important in hydrologic modeling studies [25].

It is necessary to refer to published literature or land data to understand how a generalized hypothetical ATIMR curve was generated. The values in Table 5 are thought to be a visualization of an ATIMR function generated from the authors' information—as a resulting of engineering decisions implemented in 1991 at the start of snowmelt modeling studies—from the documented results of land data or by undocumented means. However, now a review of this parameter is necessary to determine the reliability of regional snowmelt predictions [24]. Following a comprehensive review of published literature, no study was found carrying out a formal validation of the ATIMR parameter using observed data other than the studies by ¸Sengül and ˙Ispirli (2021) and Fazel et al. (2014). The first of these studies was a preliminary study of snowmelt modeling in this basin. The methodology determined by Fazel et al. (2014) represented only one year of data for certain single-point locations. When the HEC-HMS program performs full hydrological simulations on catchments it uses the temperature index methodology and restricts researchers to one ATIMR function for the whole basin. It is therefore necessary to develop an optimal area-based average ATIMR function later on and is hydrologically significant for modeling snowmelt-originated flows originating in complex mountainous terrains.

The HEC-HMS model program is capable of generating grid or polygonal area-based hydrological simulation models. The HEC-HMS program allows the creation of a meteorological model to represent the meteorological boundary conditions of a basin's physical behavior and some of the spatial and area-based variables distributed over that basin. However, published literature highlights a significant deficiency in the polygonal-based modeling of the HEC-HMS model program that is widely used and part of this study, in that only one meteorological model can be used for a basin model. Consequently, eighteen hydrological models must be created for eighteen sub-basins [79]. The meteorological model applies the climatic conditions represented by precipitation, evapotranspiration, and snowmelt, based upon the methods chosen. In basins where there are large differences in altitudes, it is impractical to apply one set of snowmelt parameters—such as the melt rate or snowmelt threshold temperature—over all the locations because of a range of factors that can include radiation effects, wind conditions, and others [70]. Snowmelt parameters would not be constant for a basin that exhibited a wide range of altitudes. This would include variables such as the water capacity of the snowpack and the threshold temperature at which precipitation occurs as snow or rain. In order to take this into consideration when the entire basin is modeled in one go using the polygonal-based method in HEC-HMS, it is necessary to enhance the temperature index model with area-based average ATIMR functions to cope with restricted availability of parameters and the increasing demand for accurate estimates for melt rates in both spatial and temporal terms.

The study conducted by Fazel et al. (2014) was originally the only approach to calculate the physical significance of the ATIMR curve beyond its manual calibration. As that study mentioned, although the ATI equation (the antecedent temperature index component of the ATIMR function) was provided, the SSARR guide did not describe the method used to generate the hypothetical ATIMR curve. The method provided by Fazel et al. (2014) for one snowmelt period at distinct station locations was subsequently developed and applied by ¸Sengül and ˙Ispirli (2021) to create ATIMR curves specific to the Kırkgöze– Çipak Basin using hourly temperatures and snow–water equivalent (SWE) data using error analysis methods recommended by Bombardelli and García (2003) obtained from the three meteorology and snow observation stations. The comparisons of both characteristics and statistical information from the snowmelt component simulation results of HEC-HMS, and the observed multivariate spatial–temporal SWE values of the region, shows a very high correlation between the generated ATIMR functions and the default SSARR values used in published literature [25].

Calibration of the other parameters used in the meteorological model used in the temperature index method were performed by considering the values in published literature [48,78], namely (PX temperature = 2 ◦C, base temperature = 0 ◦C, wet melt rate = 3.2 mm/◦C-day, rain rate limit = 1.3 mm/day, ATI-melt rate coefficient = 0.98, cold

limit = 0 mm/day, ATI cold rate coefficient = 0, water capacity = 20%, ground melt = 0 mm/day). As a result, the SWE simulations necessary to arrive efficiently at the final water budget calculations were optimized throughout the basin [25].

The area-based common ATIMR function (Figure 6) is meant to represent all three point ATIMR functions, so point values should be examined together and in relation to their land and snow altitude. For example, at the low-altitude Kö¸sk Station, due to the low amount of snowpack and early melting, there was a limited ATI value, and the ATI values of stations at higher altitudes were increased in proportion to the area-based ATI values exceeding that threshold. The values of the point and final area-based ATIMR functions measured in the Kırkgöze–Çipak Basin are shown in Table 6. The HEC-HMS modeled SWE results using the common area-based ATIMR function for the different stations are shown in Figure 7. The area-based ATIMR value of 125 ◦C-day—the last value in Table 6—is the cumulative ATI value for which the snow observed over the specified period at all station locations had completely melted. For rainfall–runoff studies to be carried out across the basin, the value had to be increased and extrapolated to account for the greater snow depths observed at higher altitudes by modifying precipitation series [25]. *Water* **2022**, *14*, x FOR PEER REVIEW 14 of 23

**Figure 6.** Point and area‐based ATIMR curves. **Figure 6.** Point and area-based ATIMR curves.



mon‐areal ATIMR functions.

*2.7. Calibration Strategy*

**Figure 7.** Snow–water equivalent (SWE) simulations computed with typical SSARR, point and com‐

Manual parameter calibration was preferred in this study due to the karstic behavior of the basin. Manual calibration begins with an appropriate estimation of the initial pa‐ rameters to run the model. The Kırkgöze–Çipak model was developed on a daily timescale over a 3‐year period between the calibration (2008 to 2010) and validation (2010 to 2011)

#### *2.7. Calibration Strategy*

Manual parameter calibration was preferred in this study due to the karstic behavior of the basin. Manual calibration begins with an appropriate estimation of the initial parameters to run the model. The Kırkgöze–Çipak model was developed on a daily timescale over a 3-year period between the calibration (2008 to 2010) and validation (2010 to 2011) periods. In the study area, there were three stream gauge stations, namely Karasu-Çipak (DS˙I, 21-01), Büyükçay-Karagöbek (DS˙I, 21-168), and Kö¸sk Dere-Kö¸sk (DS˙I, 21-152); these stations regularly performed hydrometric measurements. Calibration was carried out using the records of these stream gauge stations to simulate the flow in the simulations performed with HEC-HMS. The locations of meteorological stations and stream gauge stations are shown in Figure 3. The simulation was initiated at the beginning of autumn when the soil was almost dry. Therefore, the initial storage was assumed to be empty. Initial storage has an effect on the simulated hydrograph from a few days to a maximum of a few months [80]. However, they are insignificant for long-term water resource planning. After running the simulation, the simulated results were compared with the data observed from the stream gauge stations. **Figure 6.** Point and area‐based ATIMR curves.

*Water* **2022**, *14*, x FOR PEER REVIEW 14 of 23

**Figure 7.** Snow–water equivalent (SWE) simulations computed with typical SSARR, point and com‐ mon‐areal ATIMR functions. **Figure 7.** Snow–water equivalent (SWE) simulations computed with typical SSARR, point and common-areal ATIMR functions.

#### *2.7. Calibration Strategy* **3. Results**

Manual parameter calibration was preferred in this study due to the karstic behavior of the basin. Manual calibration begins with an appropriate estimation of the initial pa‐ rameters to run the model. The Kırkgöze–Çipak model was developed on a daily timescale over a 3‐year period between the calibration (2008 to 2010) and validation (2010 to 2011) The Kırkgöze–Çipak Basin is one where there is a lot of snowfall. In modeling the transformation process from snowfall into runoff, primarily SWE simulations were performed. For this reason, using the temperature index method, the data from the three meteorology stations at different altitudes and exposures on the basin were taken as points, and simulations were conducted. The validation criteria for these simulations are presented in Table 7. However, the fact that the basin has altitude and exposure differences due

to its complex topographical structure, and because the HEC-HMS has to simulate the entire basin with a single meteorological model, SWE simulations needed to be performed spatially while normally point simulations are performed for SWE simulations. This is because if the snow–water equivalent simulations are not performed properly, the snowmelt runoff process cannot be modeled properly. For this reason, a common melting parameter was developed for the three stations and the validity of the results is shown. The typical SSARR function and the common-areal melting parameter developed to be used in this and similar studies, along with the snow–water equivalent simulations at station locations, are compared in Figure 7.

After obtaining accurate snow–water equivalent simulations (Figure 7), basin-wide snowmelt runoff simulations were performed. As the stream gauge stations Kö¸sk Dere-Kö¸sk (DS˙I 21-152), Büyük Çay-Karagöbek (DS˙I 21-168), and Karasu-Çipak (DS˙I 21-01) in the basin were at the lower and main exit points of the basin, the runoff series taken from these stations were used for calibration. The calibrated parameters of the sub-basins of the Kırkgöze–Çipak Basin are summarized in Table 8 on the scale of the main sub-basins. The improved final runoff simulation results obtained from the stream gauge station points are presented in Figure 8 along with the observed values. The similarity between the runoff obtained as a result of the hydrological simulation and the runoff values obtained from the stream gauge stations in the basin where the simulation is carried out is very important in terms of validating simulation accuracy and reliability.

**Table 7.** The R<sup>2</sup> , root mean square error (RMSE), Nash–Sutcliffe efficiency (NSE) and Kling–Gupta efficiency (KGE) values for modeled and actual SWE (using estimated point and area-based ATIMR functions with the default SSARR curve for the snow accumulations and melting periods from 2008 to 2011).


**Table 8.** Initial and calibrated parameters for the three main sub-basins of the Kırkgöze–Çipak Basin.



#### **Table 8.** *Cont.*

It was observed that the simulations from the three-year time period simulated the real values very well both in temporal terms and statistically on the basin basis. R<sup>2</sup> , RMSE, NSE, and KGE values of runoff simulations at stream gauge stations are presented in Table 9.

**Table 9.** R 2 , RMSE, NSE, and KGE values of flow rate simulations at the stream gauge stations.


**Figure 8.** Observed and simulated flow rates at the (**a**) Karasu‐Çipak (DSİ 21‐01), (**b**) Köşk Dere‐ Köşk (DSİ 21‐152), and (**c**) Büyük Çay‐Karagöbek (DSİ 21‐168) stream gauge stations. **Figure 8.** Observed and simulated flow rates at the (**a**) Karasu-Çipak (DS˙I 21-01), (**b**) Kö¸sk Dere-Kö¸sk (DS˙I 21-152), and (**c**) Büyük Çay-Karagöbek (DS˙I 21-168) stream gauge stations.

(**c**)

#### **4. Discussion**

The observation of the stream hydrographs and the statistical analyses showed that the flowrate simulations at the basin outlet (DS˙I 21-01) were better than at the other stream gauge stations (DS˙I 21-152, DS˙I 21-168). It was observed that the simulations for the peak flow rates during the melting periods at these two stream gauge stations inside the study basin were higher than the measured values. However, the fact that the water budget calculations for the main basin outlet can only be obtained with the exaggerated simulations of related hydrographs indicates the presence of karstic formations in the basin. As a matter of fact, the presence of many springs observed in the study basin due to the general hydrogeological formations in the basin, and the fact that the groundwater model does not exhibit linear behavior, confirms that the land has a karstic character [23,81]. For this reason, while performing HEC-HMS model calibrations, the automatic calibration process was initially followed up to a point but later abandoned. Still, a manual calibration process was used in the study to reveal the actual behavior of the basin in general. In hydrological modeling studies, a model that reflects the basin characteristics well is expected to have good statistical indicators such as R<sup>2</sup> and RMSE or metric scores such as the Nash–Sutcliffe efficiency (NSE) and the Kling–Gupta efficiency (KGE). However, the reverse may not always be accurate [82–84]. Considering the modeling limitations for karstic behavior in the HEC-HMS model, this study may not provide characteristic flow simulations at the main basin outlet by automatically optimizing the parameters while evaluating the statistical and metric results of flow rate simulations with observation gauges inside the basin. In this case, a physically meaningful manual calibration of all the natural events that may cause the change in the flow simulations for the main basin outlet is required. However, it should not be ignored that the model results obtained with a more intense effort can improve the results much more; it has been concluded that the calibration and validation period results are sufficient. Undoubtedly, many clues can be brought about the actual basin behavior with the hydrological models established in the computer environment. The studies carried out will help the development of new techniques that can fully model natural behavior at every point of the watershed and for selecting or combining the appropriate models.

Similarly, in the manual calibration process of snow–water equivalent simulations, each characteristic detail of the SWE curve, especially during the melting period, was primarily modeled in a physically meaningful way to analyze the events for future studies. For example, the R<sup>2</sup> values calculated on the point scale for the Radar station in Table 7 are slightly lower than the values calculated on the areal-based values, unexpectedly, because the manually selected parameters also try to better simulate the critical rain on snow events during the melting period on high elevation zones [25].

#### **5. Conclusions**

In this study, basin characterization preprocessing was conducted with GIS-based HEC-GeoHMS, and basin and meteorological models were created. The outputs obtained were used as inputs for the hydrological simulation program HEC-HMS. The simulations for the years 2008 to 2011 were carried out with the model developed for the runoff of the Kırkgöze–Çipak basin and its sub-basins, where a significant part of the annual total runoff (70 to 80%) is formed by snowmelt.

The boundaries of the chosen Kırkgöze–Çipak basin study area were determined using the HEC-GeoHMS program, and its characterization was carried out and the model inputs were obtained for the HEC-HMS application. When determining the boundaries of a basin and its sub-basins, the outer basin boundary, and the surface stream network, the longest flow path, etc., are determined and then the whole basin is divided into sub-basins. After that, the physical parameters of these sub-basins are determined. In the next stage, a meteorological model definition is created for the climate characterization of the sub-basins.

The snowmelt rate function, which is the most effective parameter for the simulation of the snow–water equivalent during the implementation of the snowmelt model with the basin temperature index method, was primarily obtained from the locations of the three meteorology stations in the study basin. Then, these curves that were originally obtained as points were reduced to a single curve representing the study basin in general. Considering the characteristic behavior of point and areal snow–water equivalent simulations, as well as R <sup>2</sup> and RMSE values, the parameters required for snow–water equivalent simulation—one of the most important steps in simulating the flow rate for a basin where there is a lot of snowfall—were successfully integrated into the runoff model. With this physically-based approach, it has also been shown that regional studies on snow can be carried out more reliably and quickly.

After the basin-wide snow–water equivalent simulations were successfully performed, a hydrological model was created with HEC-HMS, and the runoff outputs of this model were correlated with the observed data from Kö¸sk Dere-Kö¸sk (DS˙I 21-152), Büyük Çay-Karagöbek (DS˙I 21-168), and Karasu-Çipak (DS˙I 21-01) stream gauge stations, and thus model calibration was performed.

Although the Kırkgöze–Çipak Basin is small in terms of surface area, it is a basin with a large altitude range. It is inevitable that meteorological variability will be high in such a basin. In hydrological model studies, it is very important for the accuracy and reliability of the simulations that the meteorological data distribution across the basin is in line with the real values in the field. Having the Kö¸sk, Güngörmez, and Radar meteorology stations, which are located at the appropriate altitude and location in the basin, ensures that the meteorological variable distribution was as close to reality as possible, and it also maximized the reliability of the hydrological model parameters obtained from the HEC-HMS.

As a result, in this study, it has been shown that, with the HEC-HMS hydrological model, flow rate simulations can be performed with very good R<sup>2</sup> and RMSE values and also NSE and KGE scores at the outlet of the snow-dominated, mountainous Kırkgöze Basin, which has a very complex topography. It is believed that the model parameters obtained and the methodology used will be a source for hydrological model studies to be carried out in similar mountain basins and help authorities to use water resources well, not only regionally, but also nationally and internationally.

**Author Contributions:** All authors contributed to the study conception and design. Conceptualization: S.¸S. and M.N.˙I., Methodology: S.¸S. and M.N.˙I., Software: S.¸S. and M.N.˙I., Validation and Formal analysis: S.¸S. and M.N.˙I., Writing—Original Draft: S.¸S. and M.N.˙I., Writing—Review and Editing: S.¸S., Visualization: S.¸S., Supervision: S.¸S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study has not received any funding from any institution/organization.

**Data Availability Statement:** All the data are provided in this paper, and no additional data are available to provide.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

