Multi-Objective Calibration of a Distributed Hydrological Model in a Highly Glacierized Watershed in Central Asia

Understanding glacio-hydrological processes is crucial to water resources management, especially under increasing global warming. However, data scarcity makes it challenging to quantify the contribution of glacial melt to streamflow in highly glacierized catchments such as those in the Tienshan Mountains. This study aims to investigate the glacio-hydrological processes in the SaryDjaz-Kumaric River (SDKR) basin in Central Asia by integrating a degree-day glacier melt algorithm into the macro-scale hydrological Soil and Water Assessment Tool (SWAT) model. To deal with data scarcity in the alpine area, a multi-objective sensitivity analysis and a multi-objective calibration procedure were used to take advantage of all aspects of streamflow. Three objective functions, i.e., the Nash–Sutcliffe efficiency coefficient of logarithms (LogNS), the water balance index (WBI), and the mean absolute relative difference (MARD), were considered. Results show that glacier and snow melt-related parameters are generally sensitive to all three objective functions. Compared to the original SWAT model, simulations with a glacier module match fairly well to the observed streamflow, with the Nash–Sutcliffe efficiency coefficient (NS) and R2 approaching 0.82 and an absolute percentage bias less than 1%. Glacier melt contribution to runoff is 30–48% during the simulation period. The approach of combining multi-objective sensitivity analysis and optimization is an efficient way to identify important hydrological processes and recharge characteristics in highly glacierized catchments.


Introduction
The arid region in Central Asia is heavily dependent on rainfall, snow, and glacier melt water from the Tienshan Mountains, which is considered as the 'water tower' [1][2][3].The water resources formed in the Tienshan Mountains are crucial not only for agricultural practice, but also for mining [4], hydropower generation, urban development, and several other industries [5].Generally acknowledged as sensitive indicators of climate change [6,7], glaciers have been retreating at an accelerated rate [8,9] due to global warming.Understanding the snow and glacier processes in the alpine mountains is critical for water resources management in water-limited regions.
Recently, many studies [6,[10][11][12][13][14][15] on current and future water resources in Central Asia have been conducted to gain a deeper understanding of the complex hydrological processes of the region.These studies also intend to gauge, through both statistical analysis on historical data and hydrological modeling, how the area's water resources might change in the coming years.Using statistical analysis, Zhang et al. [11] looked at the increase in streamflow in the Tienshan Mountains, which the researchers attributed to a significant increase in precipitation.Shen et al. [12], also employing statistical analysis, showed that mean streamflow has experienced a significant increase since the 1990s that is consistent with temperature and precipitation fluctuations at both annual and seasonal scales.In addition to studying hydrological processes, hydrological models can be used to predict future water resources when coupled with climate models.Sun et al. [13] applied the Hydrologiska Byrans Vattenavdelning (HBV) model to the Urimqi River Catchment in the northern slope of the Tienshan Mountains to investigate potential changes in future runoff under different glacier change scenarios, based on the RegCM3 regional climate model.Fang et al. [14], in their research on the Kaidu River, applied the Regional Climate Model (RCM) to drive the SWAT model to predict changes in runoff in the 21st century, while Immerzeel et al. [15] used the Aral Mountain model coupled with the General Circulation Model (GCM) to project average water supply variations for the Amu Darya and Syr Darya rivers.
The Aksu River, which is the largest tributary of the Tarim River, provides 70~80% of total water resources to the main stream [16,17].According to the Xinjiang Statistical Yearbook 2016 (Statistical Bureau of Xinjiang Uygur Autonomous Region, 2016), approximately 95% of the water resources are used for irrigation and support 67% of the total population.The SaryDjaz-Kumaric River (SDKR) is a major tributary of the Aksu River.Due to its importance, extensive hydrological modeling research has been conducted on the SDKR.For example, Duethmann et al. [6] assessed the effects of changes in meteorological factors on streamflow trends by applying the semi-distributed hydrological model, WASA.Wang et al. [18] applied a glacier-enhanced SWAT model to investigate the influence of temperature and precipitation alterations on snow and glacier melt and discharge, and Zhao et al. [19] applied the Variable Infiltration Capacity (VIC) model to analyze glacier variation and the responses of hydrological processes to climate change by coupling the degree-day glacier melt method with the energy-balance based macro-scale hydrological model.However, most of these hydrological model applications did not include a sensitivity analysis and were calibrated with a single-objective function.
The importance of sensitivity analysis has been addressed by many researchers [20][21][22][23][24][25] and issued in regulatory prescriptions as guidelines for modeling (e.g., European Commission, 2005; the U.S. EPA Council for Regulatory Environmental Modeling, 2003).Furthermore, since single-objective calibration can only reflect one aspect of the hydrograph (e.g., baseflow), its application only partially reflects the catchment behaviors (e.g., response of baseflow to precipitation).Hence, these applications might limit a full understanding of more complex hydrological processes.
To address this gap in the research, the present study applies the hydrological model SWAT to the SDKR basin by integrating a glacier dynamic module.A multi-objective sensitivity analysis is performed and a multi-objective calibration method is employed.The objective of this study is to understand the hydrological processes in the data-scarce alpine catchment and quantify the contributions of rainfall, snowmelt, and glacier melt to the streamflow.The paper is arranged as follows: Sections 2 and 3 introduce the study area and the methodology used; Section 4 provides the results; and Sections 5 and 6 present the discussion and conclusions, respectively.

Study Area
The SaryDjaz-Kumaric River (SDKR) basin is a typical hydrological catchment in the mid-latitude alpine regions, recharged by snow and glacier melt water and rainfall water.Located in the southern flank of the central Tienshan Mountains, the SDKR drains an area of 12,887 km 2 above the Xiehela hydrological station (Figure 1).The catchment extends westward from Kyrgyzstan into China, with most of the basin area (about 82%) situated in Kyrgyzstan.In China, the SaryDjaz-Kumaric River, also The climate of the SDKR basin is characterized as dry continental climate, with average annual temperatures ranging between −1 and −8 °C.Precipitation in the region is greatly affected by the diverse topography and sophisticated terrain, with most of the precipitation (about 86.5%) occurring from April to September [17].A great majority of the annual runoff (over 88%) occurs from May to September [26].Summer runoff is strongly influenced by precipitation as well as glacial and snow melt.According to the Randolph Glacier Inventory [27], there are 1545 glaciers in the study area, covering a total area of 2128 km 2 (16.5% of the entire catchment).

SWAT Model and Glacier Module
The Soil and Water Assessment Tool (SWAT) model [28] is a semi-distributed hydrological model.It has been widely applied to assess a variety of water management options, such as the impact of land management practices and the short-and long-term effects of climate change on hydrological processes and sediment and agricultural chemical yields [29][30][31][32].The spatial disaggregation in SWAT includes the partition of the study watershed into subbasins and into hydrological response units (HRUs), the latter of which represent the smallest basic computing element comprising a unique combination of landuse, soil, slope, and management conditions.Processes simulated by SWAT include in-subbasin processes (snow, water, sediment, water quality, etc.) and river routing (water, sediment, and in-stream water quality) to the watershed outlet through the river network.Additional details on SWAT applications can be found in [33][34][35].
To simulate glacier processes in our study area, a glacier melt module needs to be developed.The degree-day model, widely used for snowmelt simulation [19,[36][37][38][39], is thus integrated into the SWAT model.Although the glacier dynamic module includes glacier melt and glacier accumulation, we focus only on glacier melt, as glacier accumulation occurs over a long time period and is beyond the scope of the present research [40,41].
Glacier melt is controlled by the glacier surface temperature and melting rate.In our study, glacier melt is simulated for each HRU, and the meltwater is then routed to a subbasin level.Glaciers are assumed to be located at the highest elevation band.At the given HRU, a lagging factor Lgla is The climate of the SDKR basin is characterized as dry continental climate, with average annual temperatures ranging between −1 and −8 • C. Precipitation in the region is greatly affected by the diverse topography and sophisticated terrain, with most of the precipitation (about 86.5%) occurring from April to September [17].A great majority of the annual runoff (over 88%) occurs from May to September [26].Summer runoff is strongly influenced by precipitation as well as glacial and snow melt.According to the Randolph Glacier Inventory [27], there are 1545 glaciers in the study area, covering a total area of 2128 km 2 (16.5% of the entire catchment).

SWAT Model and Glacier Module
The Soil and Water Assessment Tool (SWAT) model [28] is a semi-distributed hydrological model.It has been widely applied to assess a variety of water management options, such as the impact of land management practices and the short-and long-term effects of climate change on hydrological processes and sediment and agricultural chemical yields [29][30][31][32].The spatial disaggregation in SWAT includes the partition of the study watershed into subbasins and into hydrological response units (HRUs), the latter of which represent the smallest basic computing element comprising a unique combination of landuse, soil, slope, and management conditions.Processes simulated by SWAT include in-subbasin processes (snow, water, sediment, water quality, etc.) and river routing (water, sediment, and in-stream water quality) to the watershed outlet through the river network.Additional details on SWAT applications can be found in [33][34][35].
To simulate glacier processes in our study area, a glacier melt module needs to be developed.The degree-day model, widely used for snowmelt simulation [19,[36][37][38][39], is thus integrated into the SWAT model.Although the glacier dynamic module includes glacier melt and glacier accumulation, we focus only on glacier melt, as glacier accumulation occurs over a long time period and is beyond the scope of the present research [40,41].
Glacier melt is controlled by the glacier surface temperature and melting rate.In our study, glacier melt is simulated for each HRU, and the meltwater is then routed to a subbasin level.Glaciers are assumed to be located at the highest elevation band.At the given HRU, a lagging factor L gla is used to control the influence of its glacier temperature T gla(I day ) ( • C) on day I day and on its glacier temperature T gla(I day −1) ( • C) on day I day − 1, as follows: where T av is the mean air temperature ( • C) on day I day .The amount of glacier melt in the elevation band on a given day is calculated using Equation (2): where GT mlt ( • C) is the melting temperature of the glacier, GLA melt is the amount of glacier melt on a given day (mm H 2 O), GF mlt is the melt factor for the day (mm ), GLA cov is the fraction of the glacier coverage in each HRU, T gla is the temperature of the glacier in the elevation band, and T max is the maximum air temperature on a current day in the HRU elevation band ( • C).When the air temperature is below GT mlt , the glacier will remain frozen and runoff is not generated.Conversely, when the air temperature is above GT mlt , the glacier starts to melt.The melt factor GF mlt is thought to be relatively stable over the years and is simulated as a sine function of Julian day I day as follows: where gm f mx is the melt factor which is assumed to be the largest on 7 August (mm and gm f mn is the melt factor for 7 February (mm ), which is assumed to be the smallest.

Data Collection
SWAT runs on a daily step and requires specific information about meteorological data, topographical features, soil properties, and land management practices.
Meteorological input: There are three meteorological stations (Aksu, Koilu, and TienShan stations) located within or near the watershed.The daily meteorological data (precipitation, maximum and minimum temperatures, average wind speed, and relative humidity) for the Aksu station are from the National Meteorological Information Center (http://data.cma.cn/), while daily climate data for the Koilu and TienShan stations were interpolated based on two datasets: APHRODITE (Asian Precipitation Highly Resolved Observational Data Integration Towards the Evaluation of Water Resources, http://www.chikyu.ac.jp/precip) and CRU (Climatic Research Unit, https://crudata.uea.ac.uk/cru/data/hrg/).In the interpolation, as APHRODITE is only available to 2007, CRU data were used to extend the APHRODITE dataset to 2011, using the following steps.First, the average monthly CRU and APHRODITE precipitation datasets from 1961 to 2007 were calculated.Second, the correction factor of monthly precipitation from 2008 to 2011 was obtained by correlating monthly APHRODITE and CRU datasets from 1961 to 2007, based on the assumption that the monthly means of the two datasets were equal.Finally, the APHRODITE monthly precipitation data from 2008 to 2011 were calculated based on the correction factor.This procedure was also applied to obtain daily maximum and minimum temperatures.Further details on this process can be found in Fang et al. [14] and Wang et al. [18].
Hydrological data: The observed discharge data, including daily discharge originally derived from water level data from 1997 to 2011 at the Xiehela hydrological station (1435 m, Figure 1), were obtained from the Tarim River Basin Management Bureau of Xinjiang.
Glacier data: The glacier map was selected from the second Chinese glacier inventory [27] for the Chinese region and the most recently updated Randolph Glacier Inventory [42] for the Kyrgyzstan region.There are 1536 glaciers in the SDKR basin, with a total area of 2868.2 km 2 [19]; these glaciers are mainly located at 3000-6500 m [18].According to Su et al. [43], the average thickness of the glaciers in the Aksu River basin is 181.21 m.
Digital elevation model (DEM): The 90 m DEM was downloaded from the Shuttle Radar Topography Mission (http://www2.jpl.nasa.gov/srtm/,released in 2000) of the National Aeronautics and Space Administration (NASA).The DEM was used to determine drainage area, flow direction, and basin boundary.
Soil data: The soil map for the study region in China was obtained from the Xinjiang Institute of Ecology and Geography, Chinese Academy of Sciences, while the soil map for the study region in Kyrgyzstan was obtained from the Food and Agriculture Organization of the United Nations.All soil properties (e.g., soil type, texture, depth, etc.) were derived from the Xinjiang Agricultural Bureau and Soil Survey Office.Key soil types in the study area include alpine meadow soil, subalpine meadow soil, and rocky soil (The soil map can be found in supplementary materials).

Methodology
Given the large number of distributed parameters, we applied an efficient sensitivity technique by combining the Morris and SDP methods to distinguish insensitive parameters from sensitive ones, and quantitatively analyzed important hydrological processes.We also adopted the multi-objective optimization method ε-NSGAII for model calibration.

Morris Method
The Morris method [44], which is a qualitative approach for screening out insensitive factors, is based on the concept of individual randomized One-at-a-Time.For an n-dimensional random variable X = (x 1 , . . . ,x i−1 , x i , . . . ,x n ), where each factor x i is within [0, 1], let each x i take discrete values in the set 0, 1 p−1 , 2 p−1 , . . ., 1 where p is an even integer.Then the local sensitivity of random variable X at its jth sample x 1j , . . ., x (i−1)j , x ij , . . ., x nj is then calculated as where d ij (X) is the elementary effect of the sample point jth for factor x i (where i = 1, 2, . . . ,n and j = 1, 2, . . ., m), and ∆ is a predefined increment ∆ = p 2(p−1) .For a random variable X, m samples will generate m local sensitivity indices for each factor x i .The two main sensitivity ), are normally used as sensitivity indices in the Morris method, where µ represents the degree of parameter sensitivity and σ represents the level of interaction between the factor and other factors.Further details can be found in Saltelli et al. [45] and [46][47][48][49][50][51].

State-Dependent Parameter Method (SDP)
SDP [52] is a sensitivity method that uses recursive filtering and smoothing estimation to emulate the widely used quantitative Sobol indices [53].In the Sobol method, the total sensitivity of the model (as unit 1) can be decomposed as different orders of Sobol indices of factors, as follows: where S i indicates the main effect of factor x i , representing the average reduction of output variance when factor x i is fixed, and S ij implies the combined effect from factor x i and x j , and so on.In the Sobol method, the total effect S Ti (S Ti = S i + ∑ j S ij + ∑ j ∑ k S ijk + . . .+ S 1,2,...,n , 1 ≤ i ≤ n) is often used, which is the average output variance when the value of X i remains unknown.When applying SDP, this strategy can accurately estimate S i and S ij in Equation ( 5), and has been tested in a couple of studies.In practice, S Di is defined as a 'quasi-total effect', which takes the form S Di = S i + ∑ j S ij .In general, however, it is considered an approximation of the total effect.

Multi-Objective Calibration
In hydrological modeling, single-objective calibration only considers one aspect of the hydrograph, e.g., the Water Balance Index only considers the general water balance at a catchment scale, and the Nash-Sutcliffe tends to favor flow peaks.In fact, due to the nature of hydrological modeling, matching both water balance and flow peaks with single-objective calibration can be challenging.Because of this, there has been a shift towards multi-objective calibration in hydrological modeling.
In the present study, three objective functions were used in sensitivity analysis and calibration to compare observed and simulated flows.The first one is the Nash-Sutcliffe efficiency coefficient of logarithms (LogNS) of the observed and simulated daily streamflow.It is worth mentioning that we selected LogNS instead of NS as the objective function to avoid overfitting discharge peaks.The second objective is the water balance index (WBI), which is formulated by using the mean absolute error between the simulated and observed flow accumulation curves, as recommended by Yang et al. [46].The third objective is the mean absolute relative difference (MARD), which is calculated as the logarithms of simulated and observed flow duration curves.The objective functions are shown in Equations ( 6)-( 8): where N is the total number of timesteps in the calibration period; Q i obs and Q i sim are observed and simulated outflow at time step i, respectively; log(Q obs ) is the average of logarithmic transformed observed outflows; Q A i obs and Q A i sim are the ith observed and simulated accumulated outflows, respectively; and Q P i obs and Q P i sim are the ith percentiles of observed and simulated outflow duration curves, respectively.
Model evaluation was performed after the sensitivity analysis and multi-objective optimization.Three indices-the Nash-Sutcliffe efficiency coefficient (NS), the percentage bias (PBIAS), and the correlation coefficient (R 2 )-were selected to evaluate and compare the capability of SWAT and SWAT_Glacier to simulate hydrological processes in the watershed under study.The standard equations can be found in Zhao et al. [19] and Yin et al. [57].
For these three measures, NS is favored for large flows [46] and not greater than 1.0, while PBIAS is preferred for appraising the average deviation of the simulated from the observed, with a negative result indicating an underestimation and vice-versa.|PBI AS| denotes a deviation in the model simulations, so a PBIAS value of zero is ideal (meaning there is no deviation), while R 2 represents the collinearity of observed and simulated data.
To minimize the number of distributed parameters for calibration, we adopted the aggregation method for SWAT applied in Yang et al. [58].In this method, the concept of a 'factor' is a way to change a group of distributed parameters in the SWAT model.For example, r__Sol_k is a relative change to all Sol_k (soil hydraulic conductivity) parameters, while v__Alpha_bf indicates a replacement of all Alpha_bf (baseflow decay factor) parameters.
In our study, both SWAT and SWAT_Glacier were set up using the same procedure (note that there is no glacier set up for the SWAT application).The period 1997-2001 was used for model warm-up, 2002-2007 was used for calibration, and 2008-2011 was used for validation.The parameter sensitivity analysis and optimization methods were then applied to the hydrological model, and all results were based on the SWAT_Glacier model.To minimize the number of distributed parameters for calibration, we adopted the aggregation method for SWAT applied in Yang et al. [58].In this method, the concept of a 'factor' is a way to change a group of distributed parameters in the SWAT model.For example, r__Sol_k is a relative change to all Sol_k (soil hydraulic conductivity) parameters, while v__Alpha_bf indicates a replacement of all Alpha_bf (baseflow decay factor) parameters.

Sensitivity Analysis
In our study, both SWAT and SWAT_Glacier were set up using the same procedure (note that there is no glacier set up for the SWAT application).The period 1997-2001 was used for model warmup, 2002-2007 was used for calibration, and 2008-2011 was used for validation.The parameter sensitivity analysis and optimization methods were then applied to the hydrological model, and all results were based on the SWAT_Glacier model.

Sensitivity Analysis
Gmfmn: Glacier melt factor on 7 February (mm -denotes the corresponding parameter is not involved in model calibration due to its low sensitivity, although it is involved in sensitivity analysis.The range for each parameter according to the SWAT model theoretical documentation [59]. In the SWAT_Glacier application, we can see that v__Alpha_bf and v__Gmtmp are the most sensitive factors for LogNS.Alpha_bf is a baseflow recession constant, which controls the resident time in the groundwater system.This indicates that the groundwater process is extremely important to LogNS.Gmtmp is the surface threshold temperature for glacier melt to occur.This factor affects the total runoff simulation by controlling the glacier melting rate, which means the glacier process is highly important to runoff.The following five factors were also sensitive: snow melt-related factors (v__Sftmp, v__Smtmp, and v__Smtmp), indicating that snow melt plays an important role in the hydrological process, and Gw_delay and Ch_k2, which characterize the groundwater flow mechanism and groundwater-stream water interaction.
For WBI, the most sensitive parameter is Tlaps.Because the water balance of measured and simulated daily flow series is described using WBI, Tlaps influences the temperature input in each elevation band and determines the melting rate of snow and glaciers.Thus, it is the dominant factor for streamflow.For this reason, the temperature factor v__Tlaps has a strong effect on water balance, while other factors (v__Gmtmp, v__Gmfmx, v__Sftmp, and v__Smtmp) are also sensitive due to conspicuous interaction with v__Tlaps, as revealed by the high σ values of these factors.
The conclusion regarding the dominant effect of temperature is consistent with previous research [12], while, for MARD, the sensitivity of the parameters is similar to that in LogNS.It is worth noting that the precipitation lapse rate (Plaps) is not as sensitive as glacier and snow melt parameters.In our study of the SDKR basin, we found that temperature played a more important role in discharge than precipitation, as is consistent with Duethmann et al.'s [6] findings.
In using the Morris method to compare the SWAT_Glacier and SWAT applications, similar important hydrological processes are identified.These include groundwater processes and snow melt, as indicated by the relevant factors, e.g., v__Alpha_bf and v__smtmp, with the exception of spatial variation of precipitation (v__Plaps) and glacier melt (v__Gmtmp and v__Gmftmp).However, v__Plaps is more sensitive in the SWAT application than in the SWAT_Glacier application, as it partially compensates the impact of glacier melt on streamflow.
A quantitative sensitivity analysis was performed using the SDP method.Figure 3 shows the sensitivity results for LogNS, WBI, and MARD for both SWAT_Glacier and SWAT applications.
In the SWAT_Glacier application, LogNS, S i accounted for 69.3%, with S Di contributing up to 84.1% of parameter uncertainty.The two dominating sensitivity factors (i.e., v__Tlaps and v__Gmtmp) accounted for around 50% of LogNS uncertainty.This is actually quite high and indicates that temperature and glacier melt played an important role in the streamflow simulation.For WBI, the main effect accounted for 57.8%, with S Di contributing up to 81.3% of WBI uncertainty.The remaining 18.7% was associated with higher interactions among parameters.The sensitivity results based on MARD were similar to those for LogNS.
Among the three objective functions, the sensitivity of Plaps is significantly higher in the SWAT application than in the SWAT_Glacier model.However, Plaps indicated that rainfall contributed significantly to water balance without considering the glacier hydrological process.Another interesting phenomenon is that the S i of Tlaps in the SWAT model is much lower than its value in the SWAT_Glacier model, which is complementary to the high value of the sensitivity factor Plaps.These results point to the existence of an important compensatory effect between precipitation and glacier melt in the data-sparse alpine basin.

Multi-Objective Optimization
Based on the sensitivity analysis, fourteen factors were selected for calibration by applying ε-NSGAII with the three objective functions LogNS, WBI, and MARD.
Figure 4 provides a visual overview of the Pareto solutions and the relationship among the three objective functions.The upper left corner shows the Pareto solutions scattered in a 3-D space, while the remaining three subplots illustrate projections of these solutions in 2-D subspaces.Unsurprisingly, two high and negative correlation coefficients have emerged.MARD, at r = −0.63,shows a significant negative correlation with WBI, which indicates strong trade-off interactions along the Pareto surface.Note that, along the Pareto surface, a better MARD will lead to a worse WBI, and vice-versa; the same is true for the relationship between negative LogNS and MARD.For Pareto sets, during the calibration period the average negative LogNS is −0.911 (between −0.903 and −0.917), the average WBI is 0.020 (between 0.016 and 0.029), and the average MARD is 0.044 (between 0.028 and 0.079), giving an NS value of 0.802 (between 0.781 and 0.825).The larger the LogNS value and the smaller the WBI and MARD values, the better the model performance.All of these Pareto solutions are located in the space where NS > 0.75 and WBI < 0.1, which Moriasi proposed as being the main criterion for satisfactory modeling [60].

Multi-Objective Optimization
Based on the sensitivity analysis, fourteen factors were selected for calibration by applying ε-NSGAII with the three objective functions LogNS, WBI, and MARD.
Figure 4 provides a visual overview of the Pareto solutions and the relationship among the three objective functions.The upper left corner shows the Pareto solutions scattered in a 3-D space, while the remaining three subplots illustrate projections of these solutions in 2-D subspaces.Unsurprisingly, two high and negative correlation coefficients have emerged.MARD, at r = −0.63,shows a significant negative correlation with WBI, which indicates strong trade-off interactions along the Pareto surface.Note that, along the Pareto surface, a better MARD will lead to a worse WBI, and vice-versa; the same is true for the relationship between negative LogNS and MARD.For Pareto sets, during the calibration period the average negative LogNS is −0.911 (between −0.903 and −0.917), the average WBI is 0.020 (between 0.016 and 0.029), and the average MARD is 0.044 (between 0.028 and 0.079), giving an NS value of 0.802 (between 0.781 and 0.825).The larger the LogNS value and the smaller the WBI and MARD values, the better the model performance.All of these Pareto solutions are located in the space where NS > 0.75 and WBI < 0.1, which Moriasi proposed as being the main criterion for satisfactory modeling [60].

Model Performance
As mentioned above, all of the study's Pareto solutions are satisfactory and have a very narrow range of objective functions.We chose here one solution that has the lowest negative LogNS (or highest LogNS), with corresponding parameter values provided in Table 1.  2.  1).The calibrated Tlaps in the SWAT model (−4.04 °C km −1 ) was also significantly lower than that in the

Model Performance
As mentioned above, all of the study's Pareto solutions are satisfactory and have a very narrow range of objective functions.We chose here one solution that has the lowest negative LogNS (or highest LogNS), with corresponding parameter values provided in Table 1.  2.  1).The calibrated Tlaps in the SWAT model (−4.04 • C km −1 ) was also significantly lower than that in the SWAT_Glacier model (−6.65 • C km −1 ).In contrast, the calibrated Plaps (49.98 mm km −1 ) in the SWAT model is much higher than 10.61 mm km −1 in the SWAT_Glacier model.This is due to a high precipitation compensation for SWAT, as discussed in Section 4.1.In addition, the lower Gwqmn value in the SWAT model means less groundwater storage capacity and thus faster groundwater discharge (as baseflow) as a response to soil leakage.Note that the SWAT values of Ch_k1 and Ch_k2 are higher than those calibrated by the SWAT_Glacier model, indicating that the interaction between discharge and groundwater is stronger in the SWAT model.Additionally, in terms of water balance, although the estimated volumes of general water resources from these two methods are similar, the results from SWAT are based on artificially added precipitation, whereas those from SWAT_Glacier are based on glacier melt.
In Figure 5, for SWAT_Glacier, the observed and simulated hydrographs highly match each other, and both the peaks and the baseflow are well-simulated.NS and R 2 values are higher than SWAT for both the calibration and validation periods.In the daily discharge simulations, NS and R 2 achieved 0.82 simultaneously, and PBIAS was only 0.94% during the calibration period.SWAT_Glacier also obtained good simulations during the validation period, with NS and R 2 above 0.79, and |PBI AS| far below 10% for the calibration period.
In contrast, the original SWAT model was unable to account for the hydrological processes in this watershed.A large error was found in the SWAT simulations, as illustrated by the relatively low NS values (0.74) and R 2 (0.75) and high absolute PBIAS values (10.64%) for the calibration period, and NS = 0.67, R 2 = 0.73, and PBIAS = −21.49%for the validation period.A simulated hydrograph by SWAT_Glacier could mimic the observations in the SDKR basin.According to Moriasi et al. [60], NS greater than 0.75 and PBIAS less than 10% are indicative of excellent model calibration and validation of river outflow on different time scales daily and monthly, while a PBIAS >20% is regarded as unsatisfactory.
A comparison of multi-objective and single-objective optimization methods shows that the single-objective function achieved similarly good results in one criterion but worse results in other criteria.Therefore, for instance, taking logNS as the single-objective function, we achieved 0.82 for NS (which is as good as multi-objective optimization) but a worse PBIAS and R 2 .Thus, in terms of performance, SWAT_Glacier with multi-objective optimization can be considered effective and efficient for hydrological modeling in the highly glacierized catchments.
For the winter season, the SWAT and SWAT_Glacier models simulations performed well in reproducing observations.However, for the summer, the models showed significant differences in peak flows.The simulated peaks by SWAT are far less than the observed values.This is likely due to the lack of a glacier melt process causing an underestimation of the daily discharge during the glacier melt season (May to September), which is also the rainy season.This corresponds to an underestimation of flow (blue line) in Figure 5, as rainfall and glacier melt occur during the same season.Without a glacier module, SWAT tends to overestimate rainfall to match observed flow discharge, as is demonstrated by higher Plaps in SWAT than in SWAT_Glacier.However, there is still a slight underestimation of peak flows, since one of our objective functions is LogNS.
A comparison of the results revealed that SWAT_Glacier model gave an acceptable performance for the data-sparse alpine basin.The use of multi-objective optimization in the calibration procedure effectively reduced the parameter uncertainty arising from the choice of objective functions.Figure 6 shows 95% uncertainty bands arising from conflict among the three chosen objective functions in flow simulations by SWAT_Glacier.For calibration and validation periods, the simulated daily discharge agreed quite well with observations.Overall, the variations in the simulations were similar to those in the observations for both baseflow and peaks.Here, the results also show narrower uncertainty bands, which indicate lower uncertainty.As can be seen in Figure 6c,d, a regression analysis suggests that a simulated monthly discharge shows the best agreement with the observed values.The value of R 2 for the calibration period (R 2 = 0.93) was slightly higher than the validation period (R 2 = 0.91).Finally, correlations statistics indicate that SWAT_Glacier performed well in capturing the natural monthly discharge variability adequately for both the calibration and validation periods.A comparison of the results revealed that SWAT_Glacier model gave an acceptable performance for the data-sparse alpine basin.The use of multi-objective optimization in the calibration procedure effectively reduced the parameter uncertainty arising from the choice of objective functions.Figure 6 shows 95% uncertainty bands arising from conflict among the three chosen objective functions in flow simulations by SWAT_Glacier.For calibration and validation periods, the simulated daily discharge agreed quite well with observations.Overall, the variations in the simulations were similar to those in the observations for both baseflow and peaks.Here, the results also show narrower uncertainty bands, which indicate lower uncertainty.As can be seen in Figure 6c,d, a regression analysis suggests that a simulated monthly discharge shows the best agreement with the observed values.The value of R 2 for the calibration period (R 2 = 0.93) was slightly higher than the validation period (R 2 = 0.91).Finally, correlations statistics indicate that SWAT_Glacier performed well in capturing the natural monthly discharge variability adequately for both the calibration and validation periods.

Glacier Melt Contribution
In catchments that are snow and glacier melt-dominated, temperature is the dominant driving factor affecting the melting of the snow and glaciers.Figure 7 shows different contributions of rainfall,

Glacier Melt Contribution
In catchments that are snow and glacier melt-dominated, temperature is the dominant driving factor affecting the melting of the snow and glaciers.Figure 7 shows different contributions of rainfall, snow melt, and glacier melt to streamflow during the melt season.At the beginning of the melting season, snow melts first.Then, as the temperature rises, the proportion of snow meltwater increases significantly, reaching its highest value (40%) in June.In July and August, glacier melt water accounts for 33% and 48% of the streamflow, respectively.In the SDKR basin, the main contributions of glacier and snow melt occur from June to August, which is consistent with the agricultural irrigation season in the downstream Aksu region and thus contributes to the development of downstream oasis agriculture [61].It is clear that the amount and spatiotemporal distribution of water resources have become the primary prerequisites for sustainable development of agriculture in arid areas.Glacier and snow melt water are important water sources in the region, which further underscores that glacio-hydrological processes cannot be neglected in the hydrological simulations.
The simulation results also reveal that, although glaciers cover only around 20% of the SDKR drainage area, glacier melt water contributed up to 30-48% of the water source during the simulation period.These results are consistent with a number of recent studies.For example, Zhao et al. [26] concluded that, based on the VIC model, the glacier melt contribution was 43.8%, and Chen et al. [62], in conducting hydrograph separation based on end-member mixing analysis, suggested a contribution of 59% of glacier melt water during the glacier melt period (July-August).In a related study, Yang [63] estimated a contribution of 40% of glacier melt to streamflow based on empirical coefficient methods.

On Objective Functions
In hydrological modeling, the choice of objective functions is subject to the modeling purpose.In this study, we chose three objective functions for model calibration (LogNS, WBI and MARD), as our purpose is to study the general hydrological process (i.e., sources of flow and flow frequency) that serve as water resources management in the catchment.In our estimation, these three objectives In the SDKR basin, the main contributions of glacier and snow melt occur from June to August, which is consistent with the agricultural irrigation season in the downstream Aksu region and thus contributes to the development of downstream oasis agriculture [61].It is clear that the amount and spatiotemporal distribution of water resources have become the primary prerequisites for sustainable development of agriculture in arid areas.Glacier and snow melt water are important water sources in the region, which further underscores that glacio-hydrological processes cannot be neglected in the hydrological simulations.
The simulation results also reveal that, although glaciers cover only around 20% of the SDKR drainage area, glacier melt water contributed up to 30-48% of the water source during the simulation period.These results are consistent with a number of recent studies.For example, Zhao et al. [26] concluded that, based on the VIC model, the glacier melt contribution was 43.8%, and Chen et al. [62], in conducting hydrograph separation based on end-member mixing analysis, suggested a contribution of 59% of glacier melt water during the glacier melt period (July-August).In a related study, Yang [63] estimated a contribution of 40% of glacier melt to streamflow based on empirical coefficient methods.

On Objective Functions
In hydrological modeling, the choice of objective functions is subject to the modeling purpose.In this study, we chose three objective functions for model calibration (LogNS, WBI and MARD), as our purpose is to study the general hydrological process (i.e., sources of flow and flow frequency) that serve as water resources management in the catchment.In our estimation, these three objectives can satisfy our purpose: LogNS avoids overfitting discharge peaks, WBI focuses on the general water balance, and MARD considers the frequency distribution of discharge.However, if flooding is the main study focus, one might choose difference in peaks and time difference between peaks as objective functions.A complete list of objective functions and their purpose can be found in Yang et al. [46] and Gupta et al. [64].

About SWAT_Glacier and Input Data
Currently, our SWAT_Glacier model can simulate discharge but it cannot predict when glacial lake outburst floods will occur in glacierized catchments.However, such floods are likely to have a significant adverse impact on agricultural production and downstream residents.As the global climate continues to warm, the risk of glacial lake outburst floods is increasing.Therefore, improving the SWAT_Glacier model to make it predictable will be our upcoming work, despite the challenges this task presents.
The SWAT_Glacier model applies the degree-day approach to calculate glacier melt, as it only requires temperature data.This contrasts with other energy-based models that need additional data such as solar radiation, which is not available in our study.To our knowledge, the suitability of degree-day model is case-dependent and should be verified prior to application.The related uncertainty in glacier melt contribution will be studied in a subsequent research endeavor.
In order to further improve the simulation performance of the model, it is also necessary to obtain more meteorological observations, e.g., precipitation in mountainous areas that are currently poorly observed due to rugged terrain.In the SDKR basin, over 80% of the catchment area is located in Kyrgyzstan and less than 20% in China [65], making it difficult to interpret precipitation characteristics for the entire basin using only China-based meteorological observation station.Due to the sparse meteorological station density, it is not possible to obtain daily precipitation observation data with high temporal resolution.APHRODITE, which is generally considered the best dataset available for the region [1], was employed in this study, but it only goes to 2007, which might limit many applications.Therefore, we need to further compare other datasets and integrate high-resolution meteorological data to drive hydrological models.

Conclusions
This study investigated the glacio-hydrologic processes in the SaryDjaz-Kumaric River (SDKR) basin in Central Asia.A degree-day melt module was integrated into the semi-distributed hydrological model SWAT.Its application in this data-sparse glacierized catchment was analyzed with a multi-objective sensitivity analysis and optimization and compared to the original SWAT application.
Generally, the SWAT_Glacier model did a better job of reproducing daily flows.The NS and R 2 values were all over 0.80, and the |PBIAS| was only 0.94% in the calibration period.These results contrasted sharply with the SWAT model, which showed low NS (0.74) and R 2 (0.75) values and a high absolute |PBIAS| value (10.64%) in the calibration period.Similar results also occurred for the validation period, with the SWAT_Glacier model showing NS = 0.79, R 2 = 0.80, and |PBIAS| = −5.4%,whereas the SWAT model showed lower NS (0.67) and R 2 (0.73) values (both below 0.75) and a higher |PBIAS| value (−21.49%).This suggests that the SWAT_Glacier model with a multi-objective sensitivity analysis and optimization is useful for studying glacio-hydrological processes in alpine regions where there are scarce meteorological data.In the glacierized catchment, the SWAT_Glacier

Figure 1 .
Figure 1.Location of the SaryDjaz-Kumaric River (SDKR) basin in the Tienshan Mountains and the Tarim River Basin (bottom right).

Figure 1 .
Figure 1.Location of the SaryDjaz-Kumaric River (SDKR) basin in the Tienshan Mountains and the Tarim River Basin (bottom right).

Figure 2
Figure2illustrates the sensitivity results for the initial 26 factors (listed in Table1) for the SWAT_Glacier application and the 23 factors for the SWAT application, based on the Morris method and incorporating both insensitive factors and sensitive factors.The horizontal (µ) and vertical axes (σ) represent parameter sensitivity and interaction with other parameters, respectively.

Figure 2
Figure2illustrates the sensitivity results for the initial 26 factors (listed in Table1) for the SWAT_Glacier application and the 23 factors for the SWAT application, based on the Morris method and incorporating both insensitive factors and sensitive factors.The horizontal (μ) and vertical axes (σ) represent parameter sensitivity and interaction with other parameters, respectively.

Figure 2 .
Figure 2. Multi-objective sensitivity analysis of hydrological parameters based on the Morris method for SWAT_Glacier (left column) and SWAT (right column), with μ indicating factor sensitivity and σ indicating factor nonlinearity and interactions with other factors.

Figure 2 .
Figure 2. Multi-objective sensitivity analysis of hydrological parameters based on the Morris method for SWAT_Glacier (left column) and SWAT (right column), with µ indicating factor sensitivity and σ indicating factor nonlinearity and interactions with other factors.

Figure 3 .
Figure 3. Multi-objective sensitivity analysis results based on the SDP method for SWAT_Glacier (left column) and SWAT (right column).

Figure 3 .
Figure 3. Multi-objective sensitivity analysis results based on the SDP method for SWAT_Glacier (left column) and SWAT (right column).

Figure 4 .
Figure 4.The Pareto solutions scattered in 3-D space and the projections in 2-D subspaces, with r indicating corresponding correlation coefficients of these solutions in 2-D subspace projections.
Figure 5 shows comparisons of daily observed discharge simulated with SWAT and SWAT_Glacier for the calibration (from January 2002 to December 2007) and validation (from January 2008 to December 2011) periods.The associated statistical evaluation values are listed in Table

Figure 4 .
Figure 4.The Pareto solutions scattered in 3-D space and the projections in 2-D subspaces, with r indicating corresponding correlation coefficients of these solutions in 2-D subspace projections.
Figure 5 shows comparisons of daily observed discharge simulated with SWAT and SWAT_Glacier for the calibration (from January 2002 to December 2007) and validation (from January 2008 to December 2011) periods.The associated statistical evaluation values are listed in Table

Figure 5 .
Figure 5. Observed and simulated hydrographs in the SDKR basin for the calibration and validation periods.Discharge is represented using anomaly values due to data confidentiality.Black lines indicate observed daily discharge; blue lines indicate daily simulations by the SWAT model; red lines indicate daily simulations by the SWAT_Glacier model; the black bar shows precipitation (Separate hydrographs can be found in supplementary materials).

Figure 5 . 21 Figure 6 .
Figure 5. Observed and simulated hydrographs in the SDKR basin for the calibration and validation periods.Discharge is represented using anomaly values due to data confidentiality.Black lines indicate observed daily discharge; blue lines indicate daily simulations by the SWAT model; red lines indicate daily simulations by the SWAT_Glacier model; the black bar shows precipitation (Separate hydrographs can be found in supplementary materials).Water 2019, 10, x FOR PEER REVIEW 15 of 21

Figure 6 .
Figure 6.Times series of daily observed and simulated flows for the SDKR basin for calibration (a) and validation (b) periods with the shaded area indicating a 95% uncertainty band arising from three conflicting objective functions.The relationship between observed and simulated monthly discharge is given by (c,d).

Water 2019 , 21 Figure 7 .
Figure 7.The contributions of rainfall, snow melt, and glacier melt to streamflow during the melting season (April-September).

Figure 7 .
Figure 7.The contributions of rainfall, snow melt, and glacier melt to streamflow during the melting season (April-September).

Table 1 .
Selected factors and their ranges for sensitivity analysis and calibration.Estimated values for SWAT and SWAT_Glacier are also displayed.

Table 2 .
SWAT and SWAT_Glacier applications with multi-objective optimization for both calibration and validation periods; the model performance of SWAT_Glacier based on single-objective optimization is shown for comparison purposes.

Table 2 .
SWAT and SWAT_Glacier applications with multi-objective optimization for both calibration and validation periods; the model performance of SWAT_Glacier based on single-objective optimization is shown for comparison purposes.
Calibrated parameters are different for SWAT and SWAT_Glacier applications (Table