Application of SWAT Using Snow Data and Detecting Climate Change Impacts in the Mountainous Eastern Regions of Turkey

: In recent years, the potential impacts of climate change on water resources and the hydrologic cycle have gained importance especially for snow-dominated mountainous basins. Within this scope, the Euphrates-Tigris Basin, a snow-fed transboundary river with several large dams, was selected to investigate the effects of changing climate on seasonal snow and runoff. In this study, two headwater basins of the Euphrates River, ranging in elevation between 1500–3500 m, were assigned and SWAT was employed as a hydrological modeling tool. Model calibration and validation were conducted in a stepwise manner for snow and runoff consecutively. For the snow routine, model parameters were adjusted using MODIS daily snow-covered area, achieving hit rates of more than 95% between MODIS and SWAT. Other model parameters were calibrated successively and later validated according to daily runoff, reaching a Nash-Sutcliffe efﬁciency of 0.64–0.82 in both basins. After the modeling stage, the focus was drawn to the impacts of climate change under two different climate scenarios (RCP4.5 and RCP8.5) in two 30-year projection periods (2041–2070 and 2071–2099). From the results, it is estimated that on average snow water equivalent decreases in the order of 30–39% and snow-covered days shorten by 37–43 days for the two basins until 2099. In terms of runoff, a slight reduction of at most 5% on average volume is projected but more notably, runoff center-time is expected to shift 1–2 weeks earlier by the end of the century. phase as compared to the accumulation. Not only time intervals but also aspect differences were spotted for each basin. The subbasins comprising of the north mountain ridges (facing south) have lower SWE-SCA thresholds compared to the south ridges (facing north) in both basins. These results indicate that besides snow season accumulation/depletion differences, aspect variations may also need to be considered directly or indirectly while producing the HRUs within the snow routine of the model. More details on this topic can be assessed in the future, when there is a good set of ground snow data representing different elevation bands along with the availability of ﬁner resolution satellite snow cover for a smaller-sized basin.


Introduction
Mountainous regions serve as a lifeline in socio-hydrological systems, as the origin of most rivers stems from such elevations [1][2][3][4]. Snow at high altitudes, sometimes called 'water towers' [3,5], contributes significantly to discharge downstream, and plays a major role in the hydrological cycle. Therefore, knowledge of the amount and spatial distribution of snow cover are essential [6]. However, snow measurements are particularly difficult in such areas, and to overcome the limitations of data scarcity, the utilization of satellite imagery has been an effective method for many hydrological studies since the late 1970s [7][8][9][10][11][12][13][14][15]. Because snow-dominated areas are very susceptible to increasing temperatures, snow is also a strong indicator of climate change. According to the Intergovernmental Panel on Climate Change (IPCC) report [16], there is no doubt that the climate system has been changing due to human influence; therefore, in recent years, many researchers have focused their attention on climate change impacts in snow-dominated mountainous basins [2,[17][18][19]. They commonly emphasize the need for analyzing future water availability in such regions, usually lacking sufficient ground observations. The Soil and Water Assessment Tool (SWAT) [20] has had widely successful applications in hydrological and environmental issues over many regions of the world [21][22][23][24]. Generally, hydrologic modelling in mountainous areas is challenging; moreover, [25] reported that SWAT had inadequate performance on high-elevated catchments until a new algorithm was adopted by [26], incorporating elevation bands that discretize the snowmelt process based on basin topographic controls.
Along with this development, different studies have been applied at various parts of the world on the purpose of observing orographic effects using snow-melt parameters [27][28][29][30]. In some snow-dominated basins, runoff and snow processes of different models are compared with SWAT [31,32]. Additionally, remote sensing data like MODIS images are used along with SWAT in various mountainous regions [14,33]. In these studies, ground or satellite snow data are exploited in order to validate the model representation of the snow extent both spatially and temporally. Furthermore, SWAT is also made use of as an effective hydrologic simulator when assessing the impact of climate change on different regions of the globe [34][35][36][37]. Overall, SWAT produces acceptable results and has been shown to be an efficient tool for the planning and management of water resources, including snow-dominated catchments, when compared to other models [24,[38][39][40][41].
In the current study, SWAT is utilized for the simulation of two contiguous headwater basins located in the mountainous eastern Turkey, with ArcSWAT, a GIS-based interface. There have been previous applications of SWAT in Turkey covering different fields: water quality [42,43], optimal water-management strategies [43,44], the effect of climate change on groundwater [45], and water resource potential [46,47], but the novelty of this study is its use of remotely sensed and in situ measured snow data for SWAT model calibration/validation. Hence, the novel contribution and main objectives of the current research are: (i) the application of SWAT in snow-dominated mountainous basins as a pioneer study in Turkey, (ii) utilizing remote sensing and ground snow data for model calibration and validation besides runoff, and (iii) evaluating the impacts of climate change scenarios on the mountainous headwater basins.  [48] Land Cover dataset, the major land cover types for each basin are rangelands (55.4% Murat and 52.4% Karasu) and agricultural areas (36.7% Murat and 38.8% Karasu). As for the soil map produced by the Food and Agriculture Organization of the United Nations Educational, Scientific and Cultural Organization (FAO-UNESCO) [49], different types of kastanozems and leptosols are the major soil groups in the two basins (100% Murat and 93% Karasu).

Study Area
Annual average discharges (Q) recorded at stream gauging stations of the basins are 55.03 m 3 /s (294.2 mm) for Murat and 22.38 m 3 /s (244.2 mm) for Karasu. Agri (1648 m) and Erzurum (1758 m) are two meteorological stations within the basin boundaries ( Figure 1). Annual total precipitation (P) and mean temperature (T) are around 490 mm and 6.7 • C for Agri station and 470 mm and 5.4 • C for Erzurum station. Yearly hydro-meteorological conditions for both basins are presented in Figure 2.

SWAT Model
The Soil and Water Assessment Tool (SWAT) is a comprehensive, semi-distributed, continuous, process-based hydrological model developed in the early 1990s [20]. Initially, the objective in model development was to predict the impact of management on water, sediment, and agricultural chemical yields in large ungauged basins. The current version of SWAT [50] has been extensively developed since its origin based on the Simulator for Water Resources in Rural Basins (SWRRB) model [51]. A detailed history of SWAT model development can be found in [52][53][54].
In recent years, SWAT has become a well-known and convenient model for many applications, from hydrological to environmental issues, and its area of utilization is widespread throughout the world [22,23,40,41]. On the other hand, the model's setup and calibration procedures are demanding because of several physical model inputs and parameters.
In SWAT, a basin is divided into subbasins based on topography; moreover, hydrological response units (HRUs) for each subbasin are defined as unique combinations of slope, land use, and soil. In this study, ArcSWAT GIS-based interface is used to identify subbasins and HRUs, as well as handling model input files. Major hydrological processes that include surface runoff, lateral flow, evapotranspiration, recharge, revap, return flow, infiltration, and percolation are computed in each HRU within a water-balance calculation based on four reservoirs (snow, soil, and shallow and deep aquifer) [55][56][57].
It has been proven in several studies that using elevation bands with associated snow parameters and lapse rates increases the model's success for snow-dominated mountainous basins [6,14,26,[58][59][60][61]. Among these studies, there is no common idea or recommendation on the number of elevation bands selected, which can be defined as a function of elevation or area. For this study, ten elevation bands are preferred in modelling to elaborate the model snow outputs in detail. Table 1 shows model inputs under four major classes: HRU definition, climate, calibration/validation, and climate projection data. SWAT needs some physical data in order to properly define HRUs. Shuttle Radar Topography Mission (SRTM) data [62] is utilized as a digital elevation model (DEM) in delineating the basin and river network. CORINE and FAO data are employed for generating land-use and soil type, respectively. For the climate data, daily precipitation and temperature were obtained from the Turkish Meteorological Service (MGM) gauge measurements. The optional relative humidity, solar radiation, and wind speed data were gathered from global weather data derived from Climate Forecast System Reanalysis (CFSR) [63], which is provided as an option under SWAT.

Model Inputs
For model calibration/validation, discharge and snow water equivalent (SWE) data were collected from the Turkish Hydraulic Works (DSI). Daily discharge at basin outlets and monthly snow tube SWE measurements from snow stations located in the vicinity of the basins are utilized. As for remote sensing data, MODerate Resolution Imaging Spectroradiometer (MODIS) daily snow-cover images (M*D10A1) were downloaded from NASA [64]. MODIS (Terra/Aqua) snow-cover extent was processed through a series of cloud-filtering algorithms (including a combination of Terra/Aqua, spatio-temporal, elevation, and seasonal filters) to obtain the daily cloud-free snow extent [65].
Climate projection data were gathered from the Turkish Meteorological Service (MGM) using the output of the MPI-ESM-MR Global Circulation Model whose climate parameters are dynamically downscaled to 20 km domain.

Base Model Setup
For setting up the base model without any calibration, ArcSWAT 2012.10_4.19 desktop version was used as an interface. Firstly, in the watershed delineation step, the minimum subbasin threshold area (STA) values were defined credibly [66,67]. When a certain STA was selected for each basin (

Snow Parameters Fitting Procedure
SWAT has several model parameters for calibration and for such processes; an external automatic calibration software called SWAT-CUP (SWAT-Calibration and Uncertainty Program) [68] is preferred by SWAT-users for ease and efficiency [57]. When all the model parameters were set for calibration with runoff only in SWAT-CUP, unrealistic values were achieved, especially for some of the snow parameters, as also experienced by [61]. These abnormalities were especially seen in SMFMN (minimum melt factor in the Northern Hemisphere) parameter being greater than SMFMX (maximum melt factor in the Northern Hemisphere), as well as in high/low PLAPS/TLAPS rates. These results could be attributed to either calibrating too many parameters at the same time or fitting the model parameters to flow only. As stated by [69], according to the parametrization procedure, snow parameters and lapse rates should be fitted first before flow calibration. For this reason, a temporal methodology was undertaken in this study to compare the SWE-SCA (SWAT derived Snow Water Equivalent in mm and MODIS derived Snow-Covered Area in percent) relationship.
SWE calculations in SWAT are HRU-based; therefore, values were first transformed into subbasins and later converted as basin-wide SWE outputs. On the other hand, SCA values were derived from MODIS daily cloud-free snow extent for each basin. Although these two snow components (SWE vs SCA) have different units, a time series comparison could be carried out [12,14,[70][71][72][73][74][75][76].
In our fitting procedure, the consistency between SWE and SCA was achieved by manually adjusting the SWAT snow parameters within physical limits using previous investigations on the basins [77][78][79]. Figure 3 illustrates an example of the SWE-SCA comparison for one year (2004) in each basin, including default and calibrated snow parameter sets. In this comparison, the important elements are the timing match between the beginning and end dates of snow cover as well as the accumulation and depletion phases. Figure 4 depicts SWAT-derived SWE using default/calibrated snow parameters and MODIS-derived SCA in both Murat and Karasu basins for the whole modeling period (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011). A simple hit rate performance measure was utilized to assess the comparison success as given in Equation (1). In order not to degrade the match performance from small fluctuations (due to improper satellite snow-cloud discrimination or ground partial snow cover at the beginning and end periods of the snow season), a limit of 1 mm for SWE and 5% for SCA was set. Hit where 'a' is the number of snow days for both SWE and SCA, 'd' is the number of no snow days for both SWE and SCA, 'b' is the number of snow days for SWE and no snow days for SCA, and 'c' is the number of no snow days for SWE and snow days for SCA.
Water 2021, 13, 1982 6 of 25 SWE-SCA (SWAT derived Snow Water Equivalent in mm and MODIS derived Snow-Covered Area in percent) relationship. SWE calculations in SWAT are HRU-based; therefore, values were first transformed into subbasins and later converted as basin-wide SWE outputs. On the other hand, SCA values were derived from MODIS daily cloud-free snow extent for each basin. Although these two snow components (SWE vs SCA) have different units, a time series comparison could be carried out [12,14,[70][71][72][73][74][75][76].
In our fitting procedure, the consistency between SWE and SCA was achieved by manually adjusting the SWAT snow parameters within physical limits using previous investigations on the basins [77][78][79]. Figure 3 illustrates an example of the SWE-SCA comparison for one year (2004) in each basin, including default and calibrated snow parameter sets. In this comparison, the important elements are the timing match between the beginning and end dates of snow cover as well as the accumulation and depletion phases. Figure 4 depicts SWAT-derived SWE using default/calibrated snow parameters and MODIS-derived SCA in both Murat and Karasu basins for the whole modeling period (2002-2011). A simple hit rate performance measure was utilized to assess the comparison success as given in Equation (1). In order not to degrade the match performance from small fluctuations (due to improper satellite snow-cloud discrimination or ground partial snow cover at the beginning and end periods of the snow season), a limit of 1 mm for SWE and 5% for SCA was set.
where 'a' is the number of snow days for both SWE and SCA, 'd' is the number of no snow days for both SWE and SCA, 'b' is the number of snow days for SWE and no snow days for SCA, and 'c' is the number of no snow days for SWE and snow days for SCA.   Table 3 displays the SWE-SCA comparison hit-rate scores for the calibrated model in each basin and water year. Overall, hit rate scores are above 95% both in the calibration and validation periods, indicating a good fit in terms of temporal comparison. Table 4 tabulates the suitable snow parameters and lapse rates fitted after applying the SWE-SCA temporal methodology. These values were manually calibrated, and tended to suit both basins (since the two basins are quite similar in terms of climate, topography, and land use), hence they may not be global optimums but rather a physically meaningful and trustable parameter set.   Table 3 displays the SWE-SCA comparison hit-rate scores for the calibrated model in each basin and water year. Overall, hit rate scores are above 95% both in the calibration and validation periods, indicating a good fit in terms of temporal comparison. Table 4 tabulates the suitable snow parameters and lapse rates fitted after applying the SWE-SCA temporal methodology. These values were manually calibrated, and tended to suit both basins (since the two basins are quite similar in terms of climate, topography, and land use), hence they may not be global optimums but rather a physically meaningful and trustable parameter set.

Flow-Sensitivity Analysis, Calibration, Validation and Water Balance
After manually fitting the snow parameters, the calibration procedure for the remaining model parameters was conducted utilizing SWAT-CUP according to observed runoff. Before starting the calibration, to determine the critical parameters, a sensitivity analysis is performed by SWAT-CUP using the Sequential Uncertainty Fitting Algorithm-version2 (SUFI-2) algorithm, which relies on the Latin hypercube sampling [80]. A global (all-ata-time) sensitivity analysis procedure as well as a large number of runs (around 1000 were proposed in order to identify the impact of each parameter on the objective function) was preferred to achieve more reliable results described by [23]. Table 5 tabulates the parameters considered in sensitivity analysis and their statistical results for each basin. There are two statistics employed for the assessment of sensitive parameters in SUFI-2 algorithm, named t-stat and p-value. The larger the absolute value of t-stat, the more sensitive a parameter becomes; on the other hand, the significance of parameter sensitivity increases as p-value approaches zero. At the end of the sensitivity analysis, ten common model parameters: CN2, ALPHA_BF, SOL_Z, SOL_AWC, SOL_K, SOL_BD, GW_DELAY, RCHRG_DP, GWQMN, and REVAPMN were selected as being more sensitive at each basin and set forth for the calibration process, as shown in Table 5.
For the automatic calibration phase on streamflow, the SUFI-2 algorithm of SWAT-CUP was again utilized with the selected sensitive model parameters. Suitable change methods, relative (r_) or replace (v_) were assigned for each parameter along with the upper and lower limits. The term "r" is used for the relative adjustment of a parameter within a given range and the term "v" means directly replacing the parameter value with the assigned value. Two iterations with 1500 runs each were performed with Nash-Sutcliffe Efficiency (NSE) set as the objective function. The results of the calibration process are provided in Table 6 showing the parameter value/range for each basin. During the SWAT-CUP calibration process, the degree of uncertainty was also evaluated using two statistical measures referred to as p-factor and r-factor. The p-factor is defined as the percentage of measure data captured by 95% prediction uncertainty (95 PPU), calculated at the 2.5% and 97.5% levels of cumulative distribution of an output variable obtained through Latin hypercube sampling, and the r-factor shows the average thickness of the 95 PPU band. The p-factor and r-factor are 0.94, 0.38 for Murat and are 0.79, 0.21 for Karasu basin, respectively.
In order to evaluate the goodness of fit, Nash-Sutcliffe Efficiency (NSE), Coefficient of Determination (R 2 ) and Percent Bias (PBIAS) values were determined. Table 7 shows performance statistics for SWAT modeling considering daily calibration/validation results in each basin. Although there are no explicit standards, model outcomes could be assessed according to the ratings adapted by [81]. Accordingly, almost all categories are within a "good" to "very good" performance rating, confirming a trustable model setup. For Murat basin, calibration results are slightly better than the validation period, whereas this is reversed for Karasu basin. A pictorial representation of the hydrographs for each basin is shown in Figure 5. All periods considering the calibration (2002-2007) and validation (2008)(2009)(2010)(2011) are seen in common charts for Murat and Karasu basins and the (dis)agreement can easily be detected between observed and simulated flow data for each year.
SWAT provided water balance components as model outputs on monthly and annual basis. These outputs were first compared annually in terms of total water yield (WYLD) between the simulated and observed values considering 2002-2011 water years for each basin (Table 8) where the error margin seems to be within 1-2%. Afterwards, each modelled water balance component was analyzed annually as shown in Table 9. Most of the input components appear to be similar between the two basins, although output constituents show important differences. Finally, Figure 6 depicts some of the principal water balance components (precipitation, water yield, and evapotranspiration) in monthly terms for Murat and Karasu basins. Since they are mountainous basins, most of the streamflow is generated in the spring months due to snowmelt, and evapotranspiration is high during summer. In order to evaluate the goodness of fit, Nash-Sutcliffe Efficiency (NSE), Coefficient of Determination (R 2 ) and Percent Bias (PBIAS) values were determined. Table 7 shows performance statistics for SWAT modeling considering daily calibration/validation results in each basin. Although there are no explicit standards, model outcomes could be assessed according to the ratings adapted by [81]. Accordingly, almost all categories are within a "good" to "very good" performance rating, confirming a trustable model setup. For Murat basin, calibration results are slightly better than the validation period, whereas this is reversed for Karasu basin. A pictorial representation of the hydrographs for each basin is shown in Figure 5

Snow Validation with Ground and Satellite Data
After the model had been calibrated and validated using flow data, additional validation was carried out with different types of snow data. The first of these utilized manually measured discrete snow observations. There are six snow stations located in and around the two basins under study (Figure 1). Manual snow observations with a snow tube were conducted at these snow stations once/twice a month. The station locations, names, elevations and the corresponding subbasin numbers are shown in Figure 7. Being located at different elevations, each snow station corresponded to a certain elevation zone of the nearest subbasin. In order to test the representativeness of the snow station versus elevation band of the subbasin, Figure 8 illustrates discrete SWE measure-

Snow Validation with Ground and Satellite Data
After the model had been calibrated and validated using flow data, additional validation was carried out with different types of snow data. The first of these utilized manually measured discrete snow observations. There are six snow stations located in and around the two basins under study (Figure 1). Manual snow observations with a snow tube were conducted at these snow stations once/twice a month. The station locations, names, elevations and the corresponding subbasin numbers are shown in Figure 7. Being located at different elevations, each snow station corresponded to a certain elevation zone of the nearest subbasin. In order to test the representativeness of the snow station versus elevation band of the subbasin, Figure 8 illustrates discrete SWE measurements on the continuous model SWE for selected years. Looking at Figure 8, it is seen that in Murat basin, Eleskirt and Haciomer stations show a good correspondence with their elevation band, whereas Dogangun station gives a higher observed SWE compared to its model. In Karasu basin, Hacimahmut station gives a fine agreement to its elevation; on the other hand, Guzelyayla and Yesildere stations seem to slightly overestimate their corresponding elevation range in the snow season. These outcomes are commonly observed in many situations when point measurements are assumed to represent a certain spatial/elevation extent. Overall, it can be deduced that the match between discretely measured and continuously modelled SWE may show differences from year to year and station to station, although the general trend agrees quite well. Looking at Figure 8, it is seen that in Murat basin, Eleskirt and Haciomer stations show a good correspondence with their elevation band, whereas Dogangun station gives a higher observed SWE compared to its model. In Karasu basin, Hacimahmut station gives a fine agreement to its elevation; on the other hand, Guzelyayla and Yesildere stations seem to slightly overestimate their corresponding elevation range in the snow season. These outcomes are commonly observed in many situations when point measurements are assumed to represent a certain spatial/elevation extent. Overall, it can be deduced that the match between discretely measured and continuously modelled SWE may show differences from year to year and station to station, although the general trend agrees quite well.
The second validation assessment using snow data compared SWAT-modelled SWE with MODIS satellite SCA in spatial terms. In the calibration process, a temporal analysis was conducted, resulting in a very high correspondence (over 95% match). This time a spatial analysis was undertaken to relate modelled SWE and remotely sensed SCA during a snow season. In the literature, this relation (SWE-SCA threshold, meaning the amount of modelled SWE in mm over which there is complete snow cover by satellite) has been tested by different authors, sometimes keeping the threshold constant [14,71] and variable at others [73,76] during the modeling period. In this study, HRU-based SWE calculations of SWAT were converted into subbasins and compared with MODIS cloud-free SCA for each subbasin unit. This methodology was evaluated for an example snow season of 2006 where most of the snow conditions (SWAT SWE and MODIS SCA accumulation, peak and depletion value/timing) were on average terms (Figure 4). Water 2021, 13,1982 14 of 25  The second validation assessment using snow data compared SWAT-modelled SWE with MODIS satellite SCA in spatial terms. In the calibration process, a temporal analysis was conducted, resulting in a very high correspondence (over 95% match). This time a spatial analysis was undertaken to relate modelled SWE and remotely sensed SCA during a snow season. In the literature, this relation (SWE-SCA threshold, meaning the amount of modelled SWE in mm over which there is complete snow cover by satellite) has been tested by different authors, sometimes keeping the threshold constant [14,71] and variable at others [73,76] during the modeling period. In this study, HRU-based SWE calculations of SWAT were converted into subbasins and compared with MODIS cloud-free SCA for each subbasin unit. This methodology was evaluated for an example snow season of 2006 where most of the snow conditions (SWAT SWE and MODIS SCA accumulation, peak and depletion value/timing) were on average terms (Figure 4).
Since it is not practical to visualize the results in daily steps, example periods representing snow accumulation (04 December 2005) and snow depletion (08 April 2006) were selected for each basin as depicted in Figure 9. Two different assumptions are tested here, one keeping the SWE-SCA threshold constant throughout the year, and the other varying the threshold for accumulation and depletion periods until reaching the best possible agreement between SWE and SCA grids. Figure 9 parts (a,b) and (e,f) represent the two basins in the accumulation and depletion periods, respectively, with a 15 mm SWE-SCA constant threshold as employed by [14]. On the other hand, Figure 9 parts (c,d) and (g,h) show the same intervals utilizing variable thresholds in each subbasin to have the best possible visual match between SWAT-modelled SWE and MODIS-derived SCA as proposed by [76]. The results for the daily 2006 snow season are summarized in Table 10. It is seen that even though the SWE-SCA threshold ranges are a little different for the two basins, most probably due to topographic and climatic factors, there is a clear Since it is not practical to visualize the results in daily steps, example periods representing snow accumulation (4 December 2005) and snow depletion (8 April 2006) were selected for each basin as depicted in Figure 9. Two different assumptions are tested here, one keeping the SWE-SCA threshold constant throughout the year, and the other varying the threshold for accumulation and depletion periods until reaching the best possible agreement between SWE and SCA grids. Figure 9 parts (a,b) and (e,f) represent the two basins in the accumulation and depletion periods, respectively, with a 15 mm SWE-SCA constant threshold as employed by [14]. On the other hand, Figure 9 parts (c,d) and (g,h) show the same intervals utilizing variable thresholds in each subbasin to have the best possible visual match between SWAT-modelled SWE and MODIS-derived SCA as proposed by [76]. The results for the daily 2006 snow season are summarized in Table 10. It is seen that even though the SWE-SCA threshold ranges are a little different for the two basins, most probably due to topographic and climatic factors, there is a clear variation between the accumulation and depletion periods. The thresholds were higher during the depletion phase as compared to the accumulation. Not only time intervals but also aspect differences were spotted for each basin. The subbasins comprising of the north mountain ridges (facing south) have lower SWE-SCA thresholds compared to the south ridges (facing north) in both basins. These results indicate that besides snow season accumulation/depletion differences, aspect variations may also need to be considered directly or indirectly while producing the HRUs within the snow routine of the model. More details on this topic can be assessed in the future, when there is a good set of ground snow data representing different elevation bands along with the availability of finer resolution satellite snow cover for a smaller-sized basin. north mountain ridges (facing south) have lower SWE-SCA thresholds compared to the south ridges (facing north) in both basins. These results indicate that besides snow season accumulation/depletion differences, aspect variations may also need to be considered directly or indirectly while producing the HRUs within the snow routine of the model. More details on this topic can be assessed in the future, when there is a good set of ground snow data representing different elevation bands along with the availability of finer resolution satellite snow cover for a smaller-sized basin.

Climate Change Impacts
The climate is a pivotal factor that directly affects all elements of the hydrological cycle. As mentioned in the Fifth Assessment Report of the IPCC [16], changing precipitation and temperatures are altering hydrological systems, especially in mountainous regions. In the same report, four different scenarios referred to as Representative Concentration Pathways (RCPs-RCP2.6, RCP4.5, RCP6, and RCP8.5) are emphasized depending on the greenhouse gas concentration trajectories.
In the recent years, the Turkish Meteorological Office (MGM) has conducted climate projections to reveal the impacts of climate change in different regions of the country under the project "Climate Projections with New Scenarios for Turkey" [82]. Three Global Circulation Models (GCMs), HadGEM2-ES, MPI-ESM-MR, and GFDL-ESM2M, were exploited, and climate parameters were further dynamically downscaled using RegCM4.3.4 (Regional Climate Model) based on the RCP4.5 and RCP8.5 scenarios. While 1971-2000 was used as the reference period, the projection period covers 2016-2099.
In this study, the precipitation and temperature values of MPI-ESM-MR projections using the two scenarios (RCP4.5 and RCP8.5) were utilized in SWAT, as they were seen to have the best match with the reference observation datasets [83]. The projections are treated in two 30-year sub-periods: 2041-2070 (1st period) and 2071-2099 (2nd period). According to the MPI-ESM-MR projections based on two emission scenarios during the selected two periods, the model input variables showed changes as given in Table 11 with respect to the reference period . For the two basins under study, the variations in temperature (∆Tmp) are all in the positive direction (+1.8-+4.5 • C), increasing with the more pessimistic scenario as well as period. When the change in precipitation (∆Pcp) is considered, there are more fluctuations within periods and basins. Some of the positive changes may be attributed to the Black Sea effect (several climate change projections indicate an increase or no change in precipitation especially in the eastern Black Sea region) close to where Karasu basin is situated, as seen in Figure 1. For Murat basin, precipitation shows a small decrease almost unrelated to the type of scenario. Using these projections, climate change impacts were assessed with the calibrated SWAT model over the two basins in terms of snow and discharge. For snow comparison, SWAT-averaged SWE outputs were plotted for reference, first and second periods of the two emission scenarios (RCP4.5 and RCP8.5) in Murat and Karasu basins as shown in Figure 10. It is clearly seen that for long-term averages, the change in mean SWE values (∆SWE) is negative, which means less snow will accumulate on the ground compared to what currently occurs. From the figure, it is also detectable that by the end of the century, the snow season will narrow as snow starts to accumulate in mid-November instead of late October, and it will deplete quicker and disappear at the end of May rather than June. In terms of numbers, the change in the snow-covered days (∆Snow days) is anticipated to decrease by 23 to 43 days in Murat basin and 22 to 37 days for Karasu basin in the two projection periods, respectively (Table 12). For average SWE values, a decrease of 21% and 30% is expected to occur in Murat basin during the first and second periods, respectively, whereas these values reach 26% and 39% for Karasu basin. Figure 11 depicts SWE change rates for different elevation zones in the second period (2071-2099) of the RCP8.5 scenario, when the impact is most evident. It is clearly observed that the highest variations are in the lower elevations below 1800 m.    Considering the impacts on discharge (∆Q), the volume of runoff does not seem to be affected as much as snow, showing a decrease of at most 5% (Table 13). This is probably closely related to the change in precipitation between the periods and scenarios (Table 11) that show varying trends in each basin. For discharges, a center-time (CT) concept [17] is used which marks the day of the water year when 50% of annual total runoff is reached. There is a close relation between CT and snow melt in snow-dominated basins, because a temporal shift in the CT may indicate the impact of climate change. Table 13 shows the days representing CT for the reference time, first and second projection intervals of RCP4.5 and RCP8.5 scenarios in Murat and Karasu basins. Although the reference CTs for the two basins are almost a week apart, from the numbers it is evident that CT is shifting to earlier dates in the future (Figure 12). Since the region of interest is mountainous headwaters and the climate change scenarios project a gradual increase in temperatures, as expected more dramatical shifts are foreseen in the second period of the RCP8.5 scenario. With these results, CT for Murat is predicted to shift towards mid-April, whereas in Karasu this time shift shows the end of April. Overall, a 9-to 11-day CT translation is expected for the area of study by 2099. These results agree with the outcomes of other studies in the literature [18,19,84] stating that a 1-2 week CT translation is expected for regions above 2000 m and 2-4 week shifts are anticipated for areas between 1000-2000 m, thus predicting that no permanent snow cover will exist below 1000 m by the end of the century. Considering the impacts on discharge (ΔQ), the volume of runoff does not seem to be affected as much as snow, showing a decrease of at most 5% (Table 13). This is probably closely related to the change in precipitation between the periods and scenarios (Table 11) that show varying trends in each basin. For discharges, a center-time (CT) concept [17] is used which marks the day of the water year when 50% of annual total runoff is reached. There is a close relation between CT and snow melt in snow-dominated basins, because a temporal shift in the CT may indicate the impact of climate change. Table 13 shows the days representing CT for the reference time, first and second projection intervals of RCP4.5 and RCP8.5 scenarios in Murat and Karasu basins. Although the reference CTs for the two basins are almost a week apart, from the numbers it is evident that CT is shifting to earlier dates in the future (Figure 12). Since the region of interest is mountainous headwaters and the climate change scenarios project a gradual increase in temperatures, as expected more dramatical shifts are foreseen in the second period of the RCP8.5 scenario. With these results, CT for Murat is predicted to shift towards mid-April, whereas in Karasu this time shift shows the end of April. Overall, a 9-to 11-day CT translation is expected for the area of study by 2099. These results agree with the outcomes of other studies in the literature [18,19,84] stating that a 1-2 week CT translation is expected for regions above 2000 m and 2-4 week shifts are anticipated for areas between 1000-2000 m, thus predicting that no permanent snow cover will exist below 1000 m by the end of the century.

Conclusions
This research investigates a pioneer application of SWAT in the mountainous eastern regions of Turkey. Since the study comprises high elevated areas, snow is the most significant element of the hydrological cycle. However, it is well-known that hydrological modeling in snow-dominated basins is challenging due to the scarcity of observed data, considering complex topography and harsh climatic conditions.
In this context, SWAT was chosen to model two snow-dominated headwater basins of the transboundary Euphrates River in Turkey. A stepwise calibration procedure was applied whereby firstly, model snow parameters were manually adjusted within physical limits exploiting a satellite-assisted temporal methodology. Over 95% hit rate was achieved comparing MODIS daily cloud-free snow cover images with SWAT snow outputs for the modeling period. Such a calibration method showed a valuable example of adjusting model snow parameters using snow observations rather than runoff which could be applicable for poorly gaged areas. The second part of the calibration procedure continued by applying the SUFI-2 algorithm within SWAT-CUP software. Initially, sensitive model parameters were determined and later, they were automatically calibrated according to NSE objective function on daily runoff terms. NSE calibration/validation results ranged between 0.64-0.82 for daily timescales, corresponding to 0.74-0.89 for monthly terms, respectively. These values indicate that SWAT could produce acceptable results also in mountainous catchments.
Apart from validating the model with runoff only, validation with respect to snow was also undertaken. A satisfactory trend was observed between continuous model and discrete point snow measurements during the snow season, although some point observations deviated from their corresponding elevation bands at specific periods. This could be improved by collecting more manual points or even better by setting up automatic stations measuring continuous SWE or snow depth for a more valuable direct comparison. In terms of spatially distributed snow validation, it was displayed that instead of a single SWE-SCA threshold, varying limits were demonstrated especially for accumulation and depletion periods. Furthermore, a noteworthy aspect effect was highlighted as a detailed future study that could be assessed while producing HRUs within the model setup.
Once a trustable hydrological model was attained, the impacts of climate change on the basins were investigated. Regionally downscaled MPI-ESM-MR model clearly indicated an increasing temperature trend with different emission scenarios (RCP4.5 and RCP8.5) and projection periods (2041-2070 and 2071-2099), but a variable precipitation change. Accordingly, a more drastic variation is expected to occur in terms of snow amount (30-39%) and shortening of snow days (37-43 days) rather than runoff volume. The greatest SWE change is foreseen to occur in the lower elevations (<1800 m) of the region. Although runoff volume seems to be affected less, the shift in streamflow timing to earlier dates (1-2 weeks) is a notable indication of changing climate for the headwater basins.
The results of this study indicate that climate change will have a significant impact on snow amounts and runoff regimes for the two headwater basins which may continue to increase downstream to the large storage reservoirs by the end of the century. Multihydrologic models including multi-climate projections expanding to the water-scarce lower areas should be conducted in order to draw more concrete conclusions on the outcomes. Overall, the consequences of climate change promote decision-makers to consider new plans and more sustainable management techniques in terms of water resources, agriculture, livestock, and winter tourism in the region.  Data Availability Statement: All data source details are given in "Model Inputs Section" and summarized in Table 1.