Next Article in Journal
Maritime Spatial Planning in the European Union on the Example of the Polish Part of the Baltic Sea
Next Article in Special Issue
Duality of Seasonal Effect and River Bend in Relation to Water Quality in the Chao Phraya River
Previous Article in Journal
Development of Multi-Objective Optimal Redundant Design Approach for Multiple Pipe Failure in Water Distribution System
Previous Article in Special Issue
Managing for Sustainability: The Development of Environmental Flows Implementation in China
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
State Key Laboratory of Desert and Oasis Ecology, Xinjiang Institute of Ecology and Geography, Chinese Academy of Sciences, Urumqi 830011, China
2
University of Chinese Academy of Sciences, Beijing 100049, China
3
National Institute of Water and Atmospheric Research, Christchurch 8011, New Zealand
*
Author to whom correspondence should be addressed.
Water 2019, 11(3), 554; https://doi.org/10.3390/w11030554
Submission received: 4 February 2019 / Revised: 12 March 2019 / Accepted: 13 March 2019 / Published: 17 March 2019
(This article belongs to the Special Issue Water Related Disaster and Water Environment Management)

Abstract

:
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.

1. 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: Section 2 and Section 3 introduce the study area and the methodology used; Section 4 provides the results; and Section 5 and Section 6 present the discussion and conclusions, respectively.

2. Study Area and Hydrological Model

2.1. 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 km2 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 known as the Kumaric River, is the main tributary of the Tarim River. The altitude of the SDKR basin ranges from 1435 to 7126 m above sea level, with an average elevation of 3750 m.
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 km2 (16.5% of the entire catchment).

2.2. 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 used to control the influence of its glacier temperature T g l a ( I d a y ) (°C) on day I d a y and on its glacier temperature T g l a ( I d a y 1 ) (°C) on day I d a y 1 , as follows:
T g l a ( I d a y ) = T g l a ( I d a y 1 ) · ( 1 L g l a ) + T a v · L g l a
where T a v is the mean air temperature (°C) on day I d a y . The amount of glacier melt in the elevation band on a given day is calculated using Equation (2):
G L A m e l t = { G F m l t · G L A c o v · [ T g l a + T m a x 2 G T m l t ] , T m a x > G T m l t 0 , T m a x G T m l t
where G T m l t (°C) is the melting temperature of the glacier, G L A m e l t is the amount of glacier melt on a given day (mm H2O), G F mlt is the melt factor for the day (mm H2O·°C−1·d−1), G L A c o v is the fraction of the glacier coverage in each HRU, T g l a is the temperature of the glacier in the elevation band, and T m a x is the maximum air temperature on a current day in the HRU elevation band (°C). When the air temperature is below G T m l t , the glacier will remain frozen and runoff is not generated. Conversely, when the air temperature is above G T m l t , the glacier starts to melt.
The melt factor G F mlt is thought to be relatively stable over the years and is simulated as a sine function of Julian day I d a y as follows:
G F m l t = g m f m x g m f m n 2 · s i n [ 2 π 365 · ( I d a y 128 ) ] + g m f m x + g m f m n 2
where g m f m x is the melt factor which is assumed to be the largest on 7 August (mm H2O·°C−1·d−1), and g m f m n is the melt factor for 7 February (mm H2O·°C−1·d−1), which is assumed to be the smallest.

2.3. 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 km2 [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).
Landuse data: The landuse dataset MODIS (Moderate resolution Imaging Spectroradiometer), featuring a 500 m spatial resolution, came from NASA (https://ladsweb.modaps.eosdis.nasa.gov/). The two major landuse types are bare land (74%) and water (24%, including glacier and snow), while the rest are forest and agricultural land.

3. 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.

3.1. Sensitivity Analysis Techniques

3.1.1. 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 1 j , , x ( i 1 ) j , x i j , , x n j ) is then calculated as
d i j ( X ) = f ( x 1 j , ,   x ( i 1 ) j ,   x i j + Δ , ,   x n j ) f ( x 1 j , ,   x ( i 1 ) j ,   x i j , ,   x n j ) Δ
where d i j ( 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 indices, μ ( μ i = j = 1 m | d i j | / m ) and σ ( σ i = j = 1 m ( d i j j = 1 m d i j m ) 2 / ( m 1 ) ), 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].

3.1.2. 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:
1 = i S i + i j > i S i j + + S 1 , 2 , , n
where Si indicates the main effect of factor xi, representing the average reduction of output variance when factor xi is fixed, and S i j implies the combined effect from factor xi and xj, and so on. In the Sobol method, the total effect STi ( S T i = S i + j S i j + j k S i j k + + S 1 , 2 , , n ,   1 i n ) is often used, which is the average output variance when the value of Xi remains unknown. When applying SDP, this strategy can accurately estimate Si and Sij in Equation (5), and has been tested in a couple of studies. In practice, S D i is defined as a ‘quasi-total effect’, which takes the form S D i = S i + j S i j . In general, however, it is considered an approximation of the total effect.

3.2. 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.
Among the applied multi-objective optimization algorithms, the ε-Nondominated Sorting Genetic Algorithm-II (ε-NSGAII) [54] extends the NSGAII (Nondominated Sorting Genetic Algorithm-II) by adding ε-dominance archiving and adaptive population sizing [55]. Essentially, the ε-NSGAII aims at a set of optimal Pareto solutions to handle these seemingly irreconcilable objective functions. Further details on multi-objective optimization and ε-NSGAII can be found in [46,54,55,56].
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):
l o g N S = 1 i = 1 N [ l o g ( Q o b s i ) l o g ( Q s i m i ) ] 2 i = 1 N [ l o g ( Q o b s i ) l o g ( Q o b s ) ¯ ] 2
W B I = 1 N i = 1 N | Q o b s A i Q s i m A i |
M A R D = 1 100 i = 1 N | l o g ( Q o b s P i ) l o g ( Q s i m P i ) |
where N is the total number of timesteps in the calibration period; Q o b s i and Q s i m i are observed and simulated outflow at time step i, respectively; l o g ( Q o b s ) ¯ is the average of logarithmic transformed observed outflows; Q o b s A i and Q s i m A i are the ith observed and simulated accumulated outflows, respectively; and Q o b s P i and Q s i m P i 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 (R2)—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. | P B I A S | denotes a deviation in the model simulations, so a PBIAS value of zero is ideal (meaning there is no deviation), while R2 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.

4. Results

4.1. Sensitivity Analysis

Figure 2 illustrates the sensitivity results for the initial 26 factors (listed in Table 1) 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.
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, Si accounted for 69.3%, with SDi 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 SDi 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 Si 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.

4.2. 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].

4.3. 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. 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 2.
Calibrated parameters are different for SWAT and SWAT_Glacier applications (Table 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 R2 values are higher than SWAT for both the calibration and validation periods. In the daily discharge simulations, NS and R2 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 R2 above 0.79, and | P B I A S | 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 R2 (0.75) and high absolute PBIAS values (10.64%) for the calibration period, and NS = 0.67, R2 = 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 R2. 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 R2 for the calibration period (R2 = 0.93) was slightly higher than the validation period (R2 = 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.

5. Discussion

5.1. 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.

5.2. 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].

5.3. 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.

6. 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 R2 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 R2 (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, R2 = 0.80, and |PBIAS| = −5.4%, whereas the SWAT model showed lower NS (0.67) and R2 (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 model can also be used to calculate the proportion of glacial meltwater contribution to runoff, which is not available in SWAT.
Sensitivities of complex parameters were also quantified by combining multi-objective sensitive analysis Morris and SDP to identify important hydrological processes. In this study, the groundwater process, the snow process, and the glacier melt process were identified as important hydrological processes in the SDKR basin. SWAT without a glacier module used overestimated precipitation to compensate for the crucial glacier process. Multi-objective optimization technology that reflects a compromise of flow time series, basic water balance, and flow duration curve (i.e., with objective functions negative LogNS, WBI, and MARD) for calibration clearly contributes to the improvement of model simulations.

Supplementary Materials

The following are available online at https://www.mdpi.com/2073-4441/11/3/554/s1, Figure S1: The soil map of the SaryDjaz-Kumaric River (SDKR) basin, Figure S2: Observed streamflows in the SDKR basin for the calibration and validation period, Figure S3: Simulated streamflows by SWAT_Glacier model in the SDKR basin for the calibration and validation periods, Figure S4: Simulated streamflows by SWAT in the SDKR basin for the calibration and validation periods.

Author Contributions

H.J. built up the hydrological model and wrote the first draft. G.F. did the model sensitivity analysis. J.Y. proposed the main structure of this study. Y.C. provided useful advice and revised the manuscript. All authors contributed to the final manuscript.

Funding

This research was funded by the National Natural Science Foundation of China (Grant Nos. 41630859, 41701043) and the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant No. XDA19030204).

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Lutz, A.F.; Immerzeel, W.W.; Shrestha, A.B.; Bierkens, M.F.P. Consistent increase in High Asia’s runoff due to increasing glacier melt and precipitation. Nat. Clim. Chang. 2014, 4, 587–592. [Google Scholar] [CrossRef]
  2. Immerzeel, W.W.; van Beek, L.P.H.; Bierkens, M.F.P. Climate change will affect the Asian water towers. Science 2010, 328, 1382–1385. [Google Scholar] [CrossRef] [PubMed]
  3. Chen, Y.; Li, W.; Deng, H.; Fang, G.; Li, Z. Changes in Central Asia’s water tower: Past, present and future. Sci. Rep. 2016, 6. [Google Scholar] [CrossRef] [PubMed]
  4. Liu, M.; Wen, J.; Tan, G.; Liu, G.; Wu, B. Experimental studies and pilot plant tests for acid leaching of low-grade copper oxide ores at the Tuwu Copper Mine. Hydrometallurgy 2016, 165, 227–232. [Google Scholar] [CrossRef]
  5. Bolch, T. Hydrology: Asian glaciers are a reliable water source. Nature 2017, 545, 161. [Google Scholar] [CrossRef]
  6. Duethmann, D.; Bolch, T.; Farinotti, D.; Kriegel, D.; Vorogushyn, S.; Merz, B.; Pieczonka, T.; Jiang, T.; Su, B.; Güntner, A. Attribution of streamflow trends in snow and glacier melt-dominated catchments of the Tarim River, Central Asia. Water Resour. Res. 2015, 51, 4727–4750. [Google Scholar] [CrossRef] [Green Version]
  7. Gan, R.; Luo, Y.; Zuo, Q.T.; Sun, L. Effects of projected climate change on the glacier and runoff generation in the Naryn River Basin, Central Asia. J. Hydrol. 2015, 523, 240–251. [Google Scholar] [CrossRef] [Green Version]
  8. Farinotti, D.; Longuevergne, L.; Moholdt, G.; Duethmann, D.; Mölg, T.; Bolch, T.; Vorogushyn, S.; Güntner, A. Substantial glacier mass loss in the Tien Shan over the past 50 years. Nat. Geosci. 2015, 8, 716–722. [Google Scholar] [CrossRef]
  9. Liu, Q.; Liu, S.Y. Response of glacier mass balance to climate change in the Tianshan Mountains during the second half of the twentieth century. Clim. Dyn. 2016, 46, 303–316. [Google Scholar] [CrossRef]
  10. Meng, F.; Liu, T.; Huang, Y.; Luo, M.; Bao, A.; Hou, D. Quantitative detection and attribution of runoff variations in the Aksu River Basin. Water 2016, 8, 338. [Google Scholar] [CrossRef]
  11. Zhang, Q.; Xu, C.Y.; Tao, H.; Jiang, T.; Chen, Y.D. Climate changes and their impacts on water resources in the arid regions: A case study of the Tarim River basin, China. Stoch. Environ. Res. Risk Assess. 2010, 24, 349–358. [Google Scholar] [CrossRef]
  12. Shen, Y.-J.; Shen, Y.; Fink, M.; Kralisch, S.; Chen, Y.; Brenning, A. Trends and variability in streamflow and snowmelt runoff timing in the southern Tianshan Mountains. J. Hydrol. 2018, 557, 173–181. [Google Scholar] [CrossRef]
  13. Sun, M.; Li, Z.; Yao, X.; Zhang, M.; Jin, S. Modeling the hydrological response to climate change in a glacierized high mountain region, northwest China. J. Glaciol. 2015, 61, 127–136. [Google Scholar] [CrossRef]
  14. Fang, G.H.; Yang, J.; Chen, Y.N.; Zammit, C. Comparing bias correction methods in downscaling meteorological variables for a hydrologic impact study in an arid area in China. Hydrol. Earth Syst. Sci. 2015, 19, 2547–2559. [Google Scholar] [CrossRef] [Green Version]
  15. Immerzeel, W.W.; Lutz, A.; Droogers, P. Climate Change Impacts on the Upstream Water Resources of the Amu and Syr Darya River Basins; FutureWater: Wageningen, The Netherlands, 2012. [Google Scholar]
  16. Fan, Y.; Chen, Y.; Li, W. Increasing precipitation and baseflow in Aksu River since the 1950s. Quat. Int. 2014, 336, 26–34. [Google Scholar] [CrossRef]
  17. Krysanova, V.; Wortmann, M.; Bolch, T.; Merz, B.; Duethmann, D.; Walter, J.; Huang, S.; Tong, J.; Buda, S.; Kundzewicz, Z.W. Analysis of current trends in climate parameters, river discharge and glaciers in the Aksu River basin (Central Asia). Hydrol. Sci. J. 2015, 60, 566–590. [Google Scholar] [CrossRef]
  18. Wang, X.; Luo, Y.; Sun, L.; Zhang, Y. Assessing the effects of precipitation and temperature changes on hydrological processes in a glacier-dominated catchment. Hydrol. Process. 2015, 29, 4830–4845. [Google Scholar] [CrossRef]
  19. Zhao, Q.D.; Zhang, S.Q.; Ding, Y.J.; Wang, J.; Han, H.D.; Xu, J.L.; Zhao, C.C.; Guo, W.Q.; Shangguan, D.H. Modeling hydrologic response to climate change and shrinking glaciers in the highly glacierized Kunma Like River Catchment, Central Tian Shan. J. Hydrometeorol. 2015, 16, 2383–2402. [Google Scholar] [CrossRef]
  20. Saltelli, A.; Ratto, M.; Tarantola, S.; Campolongo, F. Update 1 of: Sensitivity analysis for chemical models. Chem. Rev. 2012, 112, 1–21. [Google Scholar] [CrossRef]
  21. Annoni, P.; Bruggemann, R.; Saltelli, A. Random and quasi-random designs in variance-based sensitivity analysis for partially ordered sets. Reliab. Eng. Syst. Saf. 2012, 107, 184–189. [Google Scholar] [CrossRef]
  22. Ferretti, F.; Saltelli, A.; Tarantola, S. Trends in sensitivity analysis practice in the last decade. Sci. Total Environ. 2016, 568, 666–670. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  23. Havel, A.; Tasdighi, A.; Arabi, M. Assessing the hydrologic response to wildfires in mountainous regions. Hydrol. Earth Syst. Sci. 2018, 22, 2527–2550. [Google Scholar] [CrossRef] [Green Version]
  24. Pathiraja, S.; Anghileri, D.; Burlando, P.; Sharma, A.; Marshall, L.; Moradkhani, H. Time-varying parameter models for catchments with land use change: The importance of model structure. Hydrol. Earth Syst. Sci. 2018, 22, 2903–2919. [Google Scholar] [CrossRef]
  25. Singh, V.; Kumar Goyal, M.; Surampalli, R.Y.; Munoz-Arriola, F. Sub catchment assessment of snowpack and snowmelt change by analyzing elevation bands and parameter sensitivity in the high Himalayas. Hydrol. Earth Syst. Sci. Discuss. 2017, 1–31. [Google Scholar] [CrossRef]
  26. Zhao, Q.D.; Ye, B.S.; Ding, Y.J.; Zhang, S.Q.; Yi, S.H.; Wang, J.; Shangguan, D.H.; Zhao, C.C.; Han, H.D. Coupling a glacier melt model to the Variable Infiltration Capacity (VIC) model for hydrological modeling in north-western China. Environ. Earth Sci. 2013, 68, 87–101. [Google Scholar] [CrossRef]
  27. Guo, W.; Liu, S.; Xu, J.; Wu, L.; Shangguan, D.; Yao, X.; Wei, J.; Bao, W.; Yu, P.; Liu, Q.; et al. The second Chinese glacier inventory: Data, methods and results. J. Glaciol. 2015, 61, 357–372. [Google Scholar] [CrossRef]
  28. Arnold, J.; Srinivasan, R.; Muttiah, R.; Williams, J. Large area hydrologic modeling and assessment—Part 1: Model development. J. Am. Water Resour. Assoc. 1998, 34, 73–89. [Google Scholar] [CrossRef]
  29. Sahoo, S.; Dhar, A.; Debsarkar, A.; Kar, A. Impact of water demand on hydrological regime under climate and LULC change scenarios. Environ. Earth Sci. 2018, 77. [Google Scholar] [CrossRef]
  30. Xu, H.; Brown, D.G.; Steiner, A.L. Sensitivity to climate change of land use and management patterns optimized for efficient mitigation of nutrient pollution. Clim. Chang. 2018. [Google Scholar] [CrossRef]
  31. Bajracharya, A.R.; Bajracharya, S.R.; Shrestha, A.B.; Maharjan, S.B. Climate change impact assessment on the hydrological regime of the Kaligandaki Basin, Nepal. Sci. Total Environ. 2018, 625, 837–848. [Google Scholar] [CrossRef]
  32. Pesce, M.; Critto, A.; Torresan, S.; Giubilato, E.; Santini, M.; Zirino, A.; Ouyang, W.; Marcomini, A. Modelling climate change impacts on nutrients and primary production in coastal waters. Sci. Total Environ. 2018, 628–629, 919–937. [Google Scholar] [CrossRef]
  33. Me, W.; Abell, J.M.; Hamilton, D.P. Effects of hydrologic conditions on SWAT model performance and parameter sensitivity for a small, mixed land use catchment in New Zealand. Hydrol. Earth Syst. Sci. 2015, 19, 4127–4147. [Google Scholar] [CrossRef] [Green Version]
  34. Arnold, J.G.; Allen, P.M.; Morgan, D.S. Hydrologic model for design and constructed wetlands. Wetlands 2001, 21, 167–178. [Google Scholar] [CrossRef]
  35. Guo, T.; Gitau, M.; Merwade, V.; Arnold, J.; Srinivasan, R.; Hirschi, M.; Engel, B. Comparison of performance of tile drainage routines in SWAT 2009 and 2012 in an extensively tile-drained watershed in the Midwest. Hydrol. Earth Syst. Sci. 2018, 22, 89–110. [Google Scholar] [CrossRef] [Green Version]
  36. Hock, R. Temperature index melt modelling in mountain areas. J. Hydrol. 2003, 282, 104–115. [Google Scholar] [CrossRef] [Green Version]
  37. Radić, V.; Hock, R. Regionally differentiated contribution of mountain glaciers and ice caps to future sea-level rise. Nat. Geosci. 2011, 4, 91–94. [Google Scholar] [CrossRef]
  38. Griessinger, N.; Seibert, J.; Magnusson, J.; Jonas, T. Assessing the benefit of snow data assimilation for runoff modeling in Alpine catchments. Hydrol. Earth Syst. Sci. 2016, 20, 3895–3905. [Google Scholar] [CrossRef] [Green Version]
  39. Ayala, A.; Pellicciotti, F.; MacDonell, S.; McPhee, J.; Burlando, P. Patterns of glacier ablation across North-Central Chile: Identifying the limits of empirical melt models under sublimation-favorable conditions. Water Resour. Res. 2017, 53, 5601–5625. [Google Scholar] [CrossRef]
  40. Kang, S.; Wang, F.; Morgenstern, U.; Zhang, Y.; Grigholm, B.; Kaspari, S.; Schwikowski, M.; Ren, J.; Yao, T.; Qin, D.; et al. Dramatic loss of glacier accumulation area on the Tibetan Plateau revealed by ice core tritium and mercury records. Cryosphere 2015, 9, 1213–1222. [Google Scholar] [CrossRef] [Green Version]
  41. Kehrwald, N.M.; Thompson, L.G.; Tandong, Y.; Mosley-Thompson, E.; Schotterer, U.; Alfimov, V.; Beer, J.; Eikenberg, J.; Davis, M.E. Mass loss on Himalayan glacier endangers water resources. Geophys. Res. Lett. 2008, 35. [Google Scholar] [CrossRef] [Green Version]
  42. Pfeffer, W.; Arendt, A.; Bliss, A.; Bolch, T.; Cogley, J.; Gardner, A.; Sharp, M. The Randolph Glacier Inventory a globally complete inventory of glaciers. J. Glaciol. 2014, 60, 537–552. [Google Scholar] [CrossRef]
  43. Su, Z.; Ding, L.; Liu, C. Tianshan glacier thickness and its reserve calculation. Xinjiang Geogr. 1984, 7, 37–44. (In Chinese) [Google Scholar] [CrossRef]
  44. Morris, M.D. Factorial sampling plans for preliminary computational experiments. Technometrics 1991, 33, 161–174. [Google Scholar] [CrossRef]
  45. Saltelli, A.; Ratto, M.; Andres, T.; Campolongo, F.; Cariboni, J.; Gatelli, D.; Saisana, M.; Tarantola, S. Global Sensitivity Analysis: The Primer; John Wiley: Hoboken, NJ, USA, 2008. [Google Scholar]
  46. Yang, J.; Castelli, F.; Chen, Y. Multiobjective sensitivity analysis and optimization of distributed hydrologic model MOBIDIC. Hydrol. Earth Syst. Sci. 2014, 18, 4101–4112. [Google Scholar] [CrossRef] [Green Version]
  47. Yang, J.; Liu, Y.; Yang, W.; Chen, Y. Multi-objective sensitivity analysis of a fully distributed hydrologic model WetSpa. Water Resour. Manag. 2011, 26, 109–128. [Google Scholar] [CrossRef]
  48. Hamdia, K.M.; Ghasemi, H.; Zhuang, X.; Alajlan, N.; Rabczuk, T. Sensitivity and uncertainty analysis for flexoelectric nanostructures. Comput. Methods Appl. Mech. Eng. 2018, 337, 95–109. [Google Scholar] [CrossRef]
  49. Kalcic, M.M.; Chaubey, I.; Frankenberger, J. Defining Soil and Water Assessment Tool (SWAT) hydrologic response units (HRUs) by field boundaries. Int. J. Agric. Biol. Eng. 2015, 8, 69–80. [Google Scholar] [CrossRef]
  50. Sebek, J.; Albin, N.; Bortel, R.; Natarajan, B.; Prakash, P. Sensitivity of microwave ablation models to tissue biophysical properties: A first step toward probabilistic modeling and treatment planning. Med Phys. 2016, 43. [Google Scholar] [CrossRef]
  51. Yang, H.; Wen, J.; Wang, S.; Li, Y. Thermal design and optimization of plate-fin heat exchangers based global sensitivity analysis and NSGA-II. Appl. Therm. Eng. 2018, 136, 444–453. [Google Scholar] [CrossRef]
  52. Ratto, M.; Pagano, A.; Young, P. State Dependent Parameter metamodelling and sensitivity analysis. Comput. Phys. Commun. 2007, 177, 863–876. [Google Scholar] [CrossRef]
  53. Sobol, I.M. Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates. Math. Comput. Simul. 2001, 55, 271–280. [Google Scholar] [CrossRef]
  54. Deb, K.; Agrawal, S.; Pratap, A.; Meyarivan, T. A Fast Elitist Non-Dominated Sorting Genetic Algorithm for Multi-Objective Optimization: NSGA-II; Springer: Berlin/Heidelberg, Germany, 2000; pp. 849–858. [Google Scholar]
  55. Tang, Y.; Reed, P.M.; Kollat, J.B. Parallelization strategies for rapid and robust evolutionary multiobjective optimization in water resources applications. Adv. Water Resour. 2007, 30, 335–353. [Google Scholar] [CrossRef]
  56. Shah, R.; Reed, P. Comparative analysis of multiobjective evolutionary algorithms for random and correlated instances of multiobjective d-dimensional knapsack problems. Eur. J. Oper. Res. 2011, 211, 466–479. [Google Scholar] [CrossRef]
  57. Yin, Z.; Feng, Q.; Liu, S.; Zou, S.; Li, J.; Yang, L.; Deo, R. The spatial and temporal contribution of glacier runoff to watershed discharge in the Yarkant River Basin, Northwest China. Water 2017, 9, 159. [Google Scholar] [CrossRef]
  58. Yang, J.; Reichert, P.; Abbaspour, K.C.; Yang, H. Hydrological modelling of the chaohe basin in china: Statistical model formulation and Bayesian inference. J. Hydrol. 2007, 340, 167–182. [Google Scholar] [CrossRef]
  59. Neitsch, S.L.; Arnold, J.G.; Kiniry, J.R.; Williams, J.R. Soil and Water Assessment Tool Theoretical Documentation Version 2009; Texas Water Resources Institute: College Station, TX, USA, 2009. [Google Scholar]
  60. Moriasi, D.N.; Arnold, J.G.; Van Liew, M.W.; Bingner, R.L.; Harmel, R.D.; Veith, T.L. Model evaluation guidelines for systematic quantification of accuracy in watershed simulations. Trans. ASABE 2007, 50, 885–900. [Google Scholar] [CrossRef]
  61. Ai, Z.; Yang, Y.; Wang, Q.; Manevski, K.; Wang, Q.; Hu, Q.; Eer, D.; Wang, J. Characteristics and influencing factors of crop coefficient for drip-irrigated cotton under plastic-mulched condition in arid environment. J. Agric. Meteorol. 2018, 74, 1–8. [Google Scholar] [CrossRef] [Green Version]
  62. Chen, H.; Chen, Y.; Li, W.; Li, Z. Quantifying the contributions of snow/glacier meltwater to river runoff in the Tianshan Mountains, Central Asia. Glob. Planet. Chang. 2019, 174, 47–57. [Google Scholar] [CrossRef]
  63. Yang, Z. Glacier water resources in China. Nat. Resour. 1987, 46–55, 68. (In Chinese) [Google Scholar]
  64. Gupta, H.V.; Sorooshian, S.; Yapo, P.O. Toward improved calibration of hydrologic models: Multiple and noncommensurable measures of information. Water Resour. Res. 1998, 34, 751–763. [Google Scholar] [CrossRef] [Green Version]
  65. Shen, Y.; Wang, G.; Ding, Y.; Mao, W.; Liu, S.; Wang, S.; Mamatkanov, D.M. Changes in glacier mass balance in watershed of Sary Jaz-Kumarik Rivers of Tianshan Mountains in 1957–2006 and their impact on water resources and trend to end of the 21th century. J. Glaciol. Geocryol. 2009, 31, 790r800. [Google Scholar]
Figure 1. Location of the SaryDjaz-Kumaric River (SDKR) basin in the Tienshan Mountains and the Tarim River Basin (bottom right).
Figure 1. Location of the SaryDjaz-Kumaric River (SDKR) basin in the Tienshan Mountains and the Tarim River Basin (bottom right).
Water 11 00554 g001
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. 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.
Water 11 00554 g002
Figure 3. Multi-objective sensitivity analysis results based on the SDP method for SWAT_Glacier (left column) and SWAT (right column).
Figure 3. Multi-objective sensitivity analysis results based on the SDP method for SWAT_Glacier (left column) and SWAT (right column).
Water 11 00554 g003
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 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.
Water 11 00554 g004
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. 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 11 00554 g005
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).
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 11 00554 g006
Figure 7. The contributions of rainfall, snow melt, and glacier melt to streamflow during the melting season (April–September).
Figure 7. The contributions of rainfall, snow melt, and glacier melt to streamflow during the melting season (April–September).
Water 11 00554 g007
Table 1. Selected factors and their ranges for sensitivity analysis and calibration. Estimated values for SWAT and SWAT_Glacier are also displayed.
Table 1. Selected factors and their ranges for sensitivity analysis and calibration. Estimated values for SWAT and SWAT_Glacier are also displayed.
No.FactorUnderlying SWAT ParametersSWAT Parameter RangeEstimated Parameter Values for SWAT ApplicationEstimated Factor Values for SWAT_Glacier Application
General SWAT parameter
1v__Alpha_bfAlpha_bf: Baseflow alpha factor[0, 1]0.440.64
2v__TlapsTlaps: Temperature lapse rate (°C km−1)[−10, −2]−4.04−6.65
3v__PlapsPlaps: Precipitation lapse rate (mm km−1)[0, 200]49.9810.61
4v__Ch_k2Ch_k2: Effective hydraulic conductivity in main channel alluvium (mm h−1)[0, 500]493.90248.78
5v__Gw_delayGw_delay: Groundwater delay time (day)[0, 500]497.95334.42
6r__SlsubbsnSlsubbsn: Average slope length (m)[−0.5, 0.5]−0.46−0.49
7v__Ch_k1Ch_k1: Effective hydraulic conductivity in tributary channel alluvium (mm h−1).[0, 300]246.43179.86
8r__Sol_kSol_kl: Saturated hydraulic conductivity (mm h−1)[−0.5, 0.5]0.490.49
9r__CN2CN2: SCS runoff curve number for moisture condition[−0.5, 0.5]−0.06−0.36
10v__GwqmnGwqmn: Threshold water level in shallow aquifer for baseflow (mm)[0, 1000]180.78240.48
11v__Gw_revapGw_revap: Groundwater ‘revap’ coefficient[−0.02, 0.2]--
12v__Ch_n2Ch_n2: Manning’s ‘n’ for main channel (-)[0, 0.3]--
13r__Sol_zSol_z: Depth from soil surface to bottom of layer (mm)[−0.5, 0.5]--
14v__RevapmnRevapmn: Threshold depth of water in shallow aquifer for revap (mm)[0, 500]--
15r__Sol_awcSol_awc: Available water capacity of the soil layer (-)[−0.5, 0.5]--
16v__EscoEsco: Soil evaporation compensation factor (-)[0, 1]--
17v__OV_NOV_N: Manning’s ‘n’ for overland flow (-)[0, 30]--
18v__SurlagSurlag: Surface runoff lag time (day)[0, 24]--
19v__SmtmpSmtmp: Snow melt base temperature(°C)[−10, 10]3.263.41
20v__SftmpSftmp: Snowfall temperature (°C)[−10, 10]3.332.59
21v__SmfmxSmfmx: Snowmelt factor on 21 June (mm °C−1·d−1)[5, 10]8.92-
22v__SmfmnSmfmn: Snowmelt factor on 21 December (mm °C−1·d−1)[0, 5]--
23v__SnocovmxSnocovmx: Water content of snow cover (mm H2O)[1, 500]479.94462.72
Glacier Module Parameters
24v__GmtmpGmtmp: Glacier melt base temperature (°C)[0, 10]-0.61
25v__GmfmxGmfmx: Glacier melt factor on 7 August (mm °C−1·d−1)[5, 10]--
26v__GmfmnGmfmn: Glacier melt factor on 7 February (mm °C−1·d−1)[0, 5]--
Note: - 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].
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.
ModelPeriodFunctionsDailyMonthly
NSPBIASR2NSPBIASR2
SWATCalibrationMulti-objective0.74−10.64%0.750.88−10.54%0.90
SWAT_GlacierCalibrationMulti-objective0.820.94%0.830.931.07%0.93
SWAT_GlacierCalibrationLogNS0.82−15.24%0.740.80−15.11%0.89
SWAT_GlacierCalibrationWBI0.69−1.48%0.710.80−1.25%0.84
SWAT_GlacierCalibrationMARD0.48−19.62%0.510.58−19.39%0.63
SWATValidationMulti-objective0.67−21.49%0.730.78−21.43%0.86
SWAT_GlacierValidationMulti-objective0.79−5.40%0.800.91−5.36%0.91
SWAT_GlacierValidationLogNS0.71−17.24%0.680.73−16.72%0.82
SWAT_GlacierValidationWBI0.57−8.40%0.590.69−7.79%0.73
SWAT_GlacierValidationMARD0.39−29.17%0.490.49−28.87%0.61

Share and Cite

MDPI and ACS Style

Ji, H.; Fang, G.; Yang, J.; Chen, Y. Multi-Objective Calibration of a Distributed Hydrological Model in a Highly Glacierized Watershed in Central Asia. Water 2019, 11, 554. https://doi.org/10.3390/w11030554

AMA Style

Ji H, Fang G, Yang J, Chen Y. Multi-Objective Calibration of a Distributed Hydrological Model in a Highly Glacierized Watershed in Central Asia. Water. 2019; 11(3):554. https://doi.org/10.3390/w11030554

Chicago/Turabian Style

Ji, Huiping, Gonghuan Fang, Jing Yang, and Yaning Chen. 2019. "Multi-Objective Calibration of a Distributed Hydrological Model in a Highly Glacierized Watershed in Central Asia" Water 11, no. 3: 554. https://doi.org/10.3390/w11030554

APA Style

Ji, H., Fang, G., Yang, J., & Chen, Y. (2019). Multi-Objective Calibration of a Distributed Hydrological Model in a Highly Glacierized Watershed in Central Asia. Water, 11(3), 554. https://doi.org/10.3390/w11030554

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop