water

Article
Variation of Runoff and Runoff Components of the Upper
Shule River in the Northeastern Qinghai–Tibet Plateau under
Climate Change
Jinkui Wu 1,2 , Hongyuan Li 1,3, *, Jiaxin Zhou 4 , Shuya Tai 2,3,5 and Xueliang Wang 1,3,6

1 Key Laboratory of Ecohydrology of Inland River Basin, Northwest Institute of Eco-Environment and
Resources, Chinese Academy of Sciences, Lanzhou 730000, China; jkwu@lzb.ac.cn (J.W.);
wangxueliang@nieer.ac.cn (X.W.)
2 State Key Laboratory of Cryospheric Sciences, Northwest Institute of Eco-Environment and Resources,
Chinese Academy of Sciences, Lanzhou 730000, China; taishuya@nieer.ac.cn
3 University of Chinese Academy of Sciences, Beijing 100049, China
4 College of Geography and Environmental Engineering, Lanzhou City University, Lanzhou 730070, China;
zhoujiaxin@lzb.ac.cn
5 College of Civil Engineering, Lanzhou Jiaotong University, Lanzhou 730070, China
6 Pingliang Hydrological Station of Gansu Province, Pingliang 744000, China
* Correspondence: lihongyuan@nieer.ac.cn

Abstract: Quantifying the impact of climate change on hydrologic features is essential for the
scientific planning, management and sustainable use of water resources in Northwest China. Based
on hydrometeorological data and glacier inventory data, the Spatial Processes in Hydrology (SPHY)



model was used to simulate the changes of hydrologic processes in the Upper Shule River (USR) from


1971 to 2020, and variations of runoff and runoff components were quantitatively analyzed using
Citation: Wu, J.; Li, H.; Zhou, J.; Tai,
the simulations and observations. The results showed that the glacier area has decreased by 21.8%
S.; Wang, X. Variation of Runoff and
with a reduction rate of 2.06 km2 /a. Significant increasing trends in rainfall runoff, glacier runoff
Runoff Components of the Upper
(GR) and baseflow indicate there has been a consistent increase in total runoff due to increasing
Shule River in the Northeastern
rainfall and glacier melting. The baseflow has made the largest contribution to total runoff, followed
Qinghai–Tibet Plateau under Climate
Change. Water 2021, 13, 3357.
by GR, rainfall runoff and snow runoff, with mean annual contributions of 38%, 28%, 18% and
https://doi.org/10.3390/w13233357 16%, respectively. The annual contribution of glacier and snow runoff to the total runoff shows a
decreasing trend with decreasing glacier area and increasing temperature. Any increase of total
Academic Editors: Wenchuan Wang, runoff in the future will depend on an increase of rainfall, which will exacerbate the impact of drought
Zhongkai Feng and Mingwei Ma and flood disasters.

Received: 2 November 2021 Keywords: SPHY model; runoff components; glacier runoff; rainfall runoff; contribution
Accepted: 23 November 2021
Published: 26 November 2021

Publisher’s Note: MDPI stays neutral 1. Introduction

with regard to jurisdictional claims in
Rivers originating in the high mountains are among the most meltwater-dependent
published maps and institutional affil-
river systems in arid and semi-arid regions, where the water resources are a key restricting
iations.
factor on social, economic and ecosystem development. Glaciers, snow meltwater and
rainfall in the high mountains are the major contributors of freshwater to most of the
main river basins in Northwest China [1–5]; large human populations and agricultural
production activities depend on these resources downstream [6–8].
Copyright: © 2021 by the authors.
The areas covered by glaciers and snow are expected to be the most vulnerable to
Licensee MDPI, Basel, Switzerland.
climate change and global warming [9–11]. Glaciers are sensitive indicators of global
This article is an open access article
warming, and glacier melting is one of the major contributors to the increase of river
distributed under the terms and
runoff [12–15]. As a result of global warming, glaciers in the Qilian Mountains are retreating
conditions of the Creative Commons
Attribution (CC BY) license (https://
fast and may lead to serious alterations in the hydrologic processes of many inland river
creativecommons.org/licenses/by/
basins [16–18]. Snow cover is significantly affected by climate change. With global warming,
4.0/).
the extent, depth of snow cover and number of snow cover days in the Qilian Mountains

Water 2021, 13, 3357. https://doi.org/10.3390/w13233357 https://www.mdpi.com/journal/water

Water 2021, 13, 3357 2 of 21

tend to decrease [19–21]. However, with the increase of precipitation, the extent, depth and
days of snow cover tend to increase again [22], resulting in a complex impact of snow runoff
on river runoff. Under climate change and global warming, the precipitation in Northwest
China also shows an increasing trend [23,24]. The increase in precipitation is altitude-
dependent: the higher the altitude, the greater the rate of increase in precipitation [25],
and rapid changes in precipitation in the high mountains in Northwest China have been
observed [26]. However, the precipitation pattern in Northwest China is fundamentally
different from that in eastern China [27,28]. The increase of precipitation is mainly due
to the increase of extreme and short-term convective rainfall [29–32]. Brief, heavy rainfall
often forms flash floods, causing damage to the safety of people and property [33].
Due to the presence of glaciers, snow cover and short-term heavy rainfall, the hydro-
logic processes in the high mountains are very complex, and the accelerated shrinking and
melting of glaciers and snow, and increasing short-duration heavy rainfall, will inevitably
change the evolution of river runoff [34]. Across the high mountains in Northwest China,
there is great variation in the contributions of glaciers, snow meltwater and rainfall to river
runoff, which is poorly quantified [34,35]. The lack of understanding of the hydrologic
processes of the alpine basins in Northwest China is one of the main sources of uncertainty
in assessing the regional hydrological impacts of climate change [35].
Understanding and exploring the contribution of glacier meltwater to rivers is the
premise and basis for predicting and coping with future changes in glacier and snow
resources and their impacts in Northwest China [36]. The study of runoff generation in
alpine regions has always been one of the hot and difficult issues in hydrology studies.
Methods commonly used to quantify the proportions of runoff components in the total
runoff include the isotope tracing method [37–39], numerical simulation method [40] and
hydrological model method [34]. However, the isotope tracing method does not effectively
distinguish shallow soil water from deep groundwater, and the numerical simulation
method often only distinguishes the baseflow. Distributed hydrological models have an
advantage in runoff segmentation because of their modular approach to hydrological
processes. Common distributed hydrological models such as the VIC model [41–43], SWAT
model [44] and SRM model [45] do not refine the segmentation of runoff components,
while the spatial processes in hydrology (SPHY) model can divide the total runoff into
glacier runoff (GR), snow runoff, rainfall runoff (overland flow and shallow lateral flow)
and the baseflow [46,47]. The detailed segmentation of total runoff by the SPHY model is
helpful to reveal the essential driving factors of long-time series runoff evolution.
The Upper Shule River (USR) is commonly considered to be a typical inland river
recharged by glacier and snow meltwater. The contribution of glacier meltwater to the river
has been a hot issue that researchers have tried to quantify, especially with the accelerated
glacier shrinking due to global warming. To address this issue, simulations of total runoff
and GR based on the HBV model [48] and VIC model [41,43] have been carried out in the
Upper Shule River, but these two models only separate GR from total runoff and do not
distinguish total runoff in more detail, thus failing to reveal the essential driving factors of
the increase in total runoff. To explore the response of hydrological processes to climate
change and reveal the essential driving factors of runoff evolution in the USR, there is an
urgent requirement to quantitatively assess the impact of glaciers, snowmelt and rainfall
on the evolution of river runoff under climate change. Therefore, this study combined
observed hydrometeorological data and multi-period glacier inventory data to simulate the
variation of runoff and runoff components in the USR using the SPHY model, with three
specific objectives: (1) to explore the response of hydrologic features to climate change and
the evolution of runoff, (2) to quantify the contributions of runoff components to the total
runoff and the evolution trend of runoff components and (3) to reveal the essential driving
factors of runoff evolution.
Water 2021,13,
Water2021, 13,3357
x FOR PEER REVIEW 3 3ofof2121

2. Data and Methods

2.
In this
this study, anomaly
anomaly analysis of temperature,
analysis of temperature,precipitation
precipitationandandrunoff
runoffininthe
theUSR
USR
was carried
was carried out based on station-observed
station-observed temperature,
temperature,precipitation
precipitationand andrunoff
runoffdata.
data.
Driving
Driving data
dataincluding
includingair
airtemperature
temperatureand precipitation,
and along
precipitation, with
along glacier
with thickness
glacier data
thickness
based on the
data based onGlabtop2 model,
the Glabtop2 werewere
model, prepared and obtained
prepared according
and obtained to the
according torequirements
the require-
of the SPHY
ments of the model. Based Based
SPHY model. on theon
calibration and validation
the calibration of theofSPHY
and validation model,
the SPHY using
model,
the observed
using runoff
the observed and and
runoff multi-period inventoried
multi-period inventoriedglacier
glacierarea,
area,the
the runoff and
andrunoff
runoff
components
components were simulated for the periodperiod of
of 1971
1971 to
to2020.
2020.

2.1.
2.1. Study
Study Area
Area
The USR is located in the the western
western partpartofofthe
theQilian
QilianMountains
Mountainson onthe
thenortheast
northeastedge
edge
of the Qinghai–Tibetan Plateau. It is a runoff formation and conservation
of the Qinghai–Tibetan Plateau. It is a runoff formation and conservation area for the area for the Shule
River, one ofone
Shule River, the of
three
the major inland
three major riversrivers
inland of theofHexi Corridor
the Hexi in Northwest
Corridor China.
in Northwest The
China.
USR is located above the National Hydrological Station of Changmabao,
The USR is located above the National Hydrological Station of Changmabao, which co- which covers an
area 2 ◦ 0 ◦ 0 ◦ 0 ◦ 0
vers of
an10,961
area ofkm (92 km
10,961 54 ~99
2 14 ~E,99°14′
(92°54′ 38 36E,~41 34 ~N),
38°36′ withN),
41°34′ an with
altitude that ranges
an altitude thatbetween
ranges
2038 and 2038
between 5792 and
m a.s.l.
5792(Figure
m a.s.l.1).
(Figure 1).

Figure 1.
Figure 1. Overview of the
the Upper
Upper Shule
ShuleRiver
River(USR)
(USR)and
andthe
thelocations
locationsofofhydrological
hydrologicaland
andmeteoro-
meteoro-
logical stations.
logical stations.

Due to the special geographical

geographical location,
location, the
the climate
climatecharacteristics
characteristicsofofthe theUSR
USRare are
characterized by both an alpine climate
characterized climate and
and aa temperate
temperate continental
continentalclimate.
climate.TheThemean
mean
annual temperature
annual temperature isisapproximately
approximately−1.62 −1.62°C,◦ and
C, andthe the
mean meanannual precipitation
annual is ap-is
precipitation
proximately 235.2
approximately mm.
235.2 mm.TheTheUSRUSR reaches its highest
reaches temperature
its highest in Julyinand
temperature Julylowest tem-
and lowest
perature in January,
temperature in January,and the
andprecipitation is extremely
the precipitation uneven
is extremely duringduring
uneven the year, with
the 91.4%
year, with
of the of
91.4% precipitation concentrated
the precipitation betweenbetween
concentrated May and September
May during the
and September warmthe
during season.
warm
The mean
season. Theannual runoff of
mean annual the USR
runoff of thefrom
USR1961
from to1961
2020tois2020
10.54is×10.54
108 m×3/a, 8 m3the
10and /a, runoff
and the
has significant seasonal characteristics, with 68.4% of runoff concentrated
runoff has significant seasonal characteristics, with 68.4% of runoff concentrated between between June
and September.
June The USR
and September. TheisUSR
a typical river fed
is a typical by glaciers
river and snow
fed by glaciers and meltwater. The glacier
snow meltwater. The
area ofarea
glacier the USR
of theisUSR
the largest in thein
is the largest Qilian Mountains,
the Qilian and the
Mountains, andexisting glaciers
the existing covercover
glaciers an
area
an of 409
area kmkm
of 409 2 based
2 based on the
on theSecond Chinese
Second ChineseGlacier Inventory
Glacier Inventory[49–51], accounting
[49–51], for
accounting
about
for 3.73%
about of the
3.73% ofbasin area. area.
the basin Due toDue the to
high
thealtitude and cold
high altitude climate,
and the mainthe
cold climate, vegeta-
main
tion types in
vegetation the USR
types in theare
USRalpine meadows,
are alpine alpine alpine
meadows, grasslands and desert
grasslands and grasslands.
desert grasslands.
Water 2021, 13, 3357 4 of 21

2.2. Data Sources

The hydrometeorological data used to analyze the variation of regional climate change
and evolution of runoff in the USR, and the forcing data for the SPHY model, primarily
included the temperature, precipitation and runoff. Among them, annual and monthly
mean temperature, annual and monthly precipitation and annual and monthly runoff were
used to analyze the inter-annual and intra-annual variation of temperature, precipitation
and runoff in the USR. Daily mean, maximum and minimum temperatures and daily
precipitation were the primary forcings for the SPHY model, and the monthly runoff was
used to calibrate and validate the performance of the SPHY model.
Temperature and precipitation data were obtained from one national hydrological
station, two national meteorological stations and six automatic meteorological stations,
and runoff data were obtained from the national hydrological station data recorded in
the Hydrological Yearbook of the People’s Republic of China. The details of each me-
teorological and hydrological station are given in Table 1 and Figure 1. Since the time
series of temperature and precipitation data were not consistent among meteorological and
hydrological stations, the required data were fitted and unified into the same time series
based on the monthly lapse rate of temperature and the monthly gradient of precipitation
in the USR in this study [41].

Table 1. Information on meteorological and hydrological stations in the USR and adjacent areas.

Station Name Lat/N Lon/E Elevation/m Data Type Station Type Period
CMB 39◦ 490 96◦ 510 2080 DR/Prec/Temp NHS 2000~2020
YM 40◦ 160 97◦ 020 1526 Prec/Temp NMS 1961~2020
TL 38◦ 520 98◦ 220 3360 Prec/Temp NMS 1961~2015
LHG 39◦ 300 96◦ 300 4180 Prec/Temp AWS 2009~2020
LQZ 39◦ 310 96◦ 580 2433 Prec/Temp AWS 2009~2020
YEH 39◦ 260 97◦ 060 2868 Prec AWS 2006~2020
GH 38◦ 500 97◦ 430 3443 Prec/Temp AWS 2009~2020
SMT 38◦ 250 98◦ 180 3885 Prec/Temp AWS 2009~2020
MK 38◦ 340 98◦ 220 4156 Prec AWS 2009~2020
Abbreviations: DR, daily runoff; Prec, precipitation; Temp, temperature; NHS, national hydrological station; NMS,
national meteorological station; AWS, automatic weather station.

In addition to temperature and precipitation as the basic forcing data, the forcing
data of the SPHY model also included topography data of a digital elevation model
(DEM), physical parameters of soil, land use data and glacier outlines and thickness. The
topography data were HydroSHEDS with a spatial resolution of 90 m based on SRTM DEM
from the United States Geological Survey (USGS), including hydrological information such
as river networks, basin boundaries and confluence directions [52]. Physical parameters of
soil include the field capacity (mm/mm), saturated water content (mm/mm), permanent
wilting point (mm/mm), wilting point (mm/mm) and saturated hydraulic conductivity
(mm/d) for the rootzone layer, and field capacity (mm/mm), saturated water content
(mm/mm) and saturated hydraulic conductivity (mm/d) for the subzone layer. Physical
parameter maps of soil with a 500 m spatial resolution were obtained from the Chinese soil
dataset based on Harmonized World Soil Database version 1.1 (HWSD v1.1), provided by
the National Cryosphere Desert Data Center [53]. Land use data were derived from the
Globcover dataset with various crop coefficients; this dataset is available freely and can
be obtained from the ESA website in the GeoTIFF format and the World Geodetic System
1984 (WGS 84) coordinate system at 10 arc-seconds (~300 m) of spatial resolution [54].
Glacier outline data were derived from the First and Second Chinese Glacier Invento-
ries provided by the National Cryosphere Desert Data Center [51], and a dataset of glacier
outlines over the Qilian Mountain area (1980–2015) was provided by the National Tibetan
Plateau Data Center [55]. The First Chinese Glacier Inventory in the USR was based on
aerial photographs taken in 1956, 1957 and 1966, and the results showed that the glacier
area in the USR was about 473 km2 , accounting for 4.32% of the USR area. LANDSAT
Water 2021, 13, 3357 5 of 21

TM/ETM + and ASTER images in 2006, 2007 and 2009 were used in the Second Chinese
Glacier Inventory in the USR, which showed that the glacier area was about 409 km2 ,
accounting for 3.73% of the USR area. The estimation of glacier thickness will be detailed
in Section 2.3.

2.3. Estimation of Glacier Thickness

The ice thickness of a glacier is a fundamental prerequisite for the SPHY model, to
simulate glacier retreat, ice redistribution and GR. In this study, Glacier Bed Topography
model version 2 (GlabTop2) was used to estimate the ice thickness of the glacier. The ice
thickness is inferred from the surface slope by parameterizing the basal shear stress (τ) as a
function of the elevation range (∆H) of each glacier [56–58]. A detailed formula explanation
of GlabTop2 is depicted as follows:
τ
h= (1)
f ρg sin(α)

0.5 + 159.8∆H − 43.5(∆H )2 ∆H ≤ 1.6 km


τ= (2)
150 ∆H > 1.6 km
where h (m) is the ice thickness of randomly picked DEM cells within the glacierized area,
τ (kPa) is the basal shear stress parameterized with the glacier elevation range [59], f is
the shape factor, set to 0.8 for mountain glaciers [60], g (9.81 m/s2 ) is the gravitational
acceleration and α refers to the surface slope. The calculation of ice thickness is grid-based
and requires a DEM and glacier mask as the inputs, and the ice thickness distribution
for all glacier cells is interpolated from the ice thickness at random cells and from the ice
thickness at the glacier margins, known to be zero. For more details of the GlabTop2 model,
refer to the literature [56].

2.4. Description of SPHY Model

2.4.1. Model Principle
The model used in the simulation of runoff and runoff components in the USR was
the SPHY model, a fully spatially distributed, high-resolution cryospheric–hydrological
leaky-bucket type of model, applied on a cell-by-cell basis on flexible spatial and temporal
scales, under various physiographical and hydroclimatic conditions [46,47]. The SPHY
model is based on the PCRaster dynamic modeling framework [61–63] and describes the
main terrestrial hydrologic processes in a conceptual way so that changes in storage and
flux can be assessed over time and space [47].
The structure of the SPHY model can be described from two dimensions: horizontal
and vertical. For the horizontal dimension, the SPHY model divides the land surface into
a glacierized area and non-glacierized area. The snow and rain falling on the glacier will
be taken into glacier processes, while the snow and rain falling on the non-glacierized
area will be taken into snow processes and rainfall-runoff processes, respectively. For the
vertical dimension, the soil column structure is similar to the VIC model [64,65], with two
upper soil storage areas and a third for groundwater storage in the non-glacierized area,
while there is only one for groundwater storage in the glacierized area. The horizontal
and vertical structure means that the SPHY model can simulate the formation processes of
different runoff components with flexible modules.
Depending on the modules being switched on in this study, four runoff components of
GR, snow runoff, rainfall runoff, and baseflow were simulated. Rainfall runoff is the sum
of surface runoff and lateral flow from the first soil layer, with the groundwater module
used. After simulating the four runoff components, the total runoff is eventually calculated
by adding together the four runoff components:

TR = GR + SR + RR + BF (3)
 sum of surface runoff and lateral flow from the first soil layer, with the groundwater mod-
ule used. After simulating the four runoff components, the total runoff is eventually cal-
Water 2021, 13, 3357 culated by adding together the four runoff components: 6 of 21

𝑇𝑅 = 𝐺𝑅 + 𝑆𝑅 + 𝑅𝑅 + 𝐵𝐹 (3)
where TR, GR, SR, RR and BF are the total runoff, GR, snow runoff, rainfall runoff and
where TR, GR, SR, RR and BF are the total runoff, GR, snow runoff, rainfall runoff and
baseflow, respectively.
baseflow, respectively.
2.4.2.
2.4.2.Mass
MassConserving
Conservingof ofGlacier
Glacier
The
TheSPHY
SPHYmodelmodeladopts
adoptstwo twoschemes
schemes to to simulate
simulate thethe glacier
glacier processes
processes based
based on on the
the
degree-day factor method. In the first scheme, the glaciers are
degree-day factor method. In the first scheme, the glaciers are implemented as a fixed implemented as a fixed
mass
massgenerating
generating glacier
glacier melt,
melt, where
where thethe snow
snow and and rainfall
rainfall on
on the
the glacier
glacier are
are not
not taken
taken into
into
account. In the second scheme, the glacier process is simulated based
account. In the second scheme, the glacier process is simulated based on mass conservation on mass conserva-
tion
withwith glacier
glacier retreat
retreat and and ice redistribution,
ice redistribution, in which
in which the snow
the snow and rainfall
and rainfall on the onglacier
the glac-
are
ier are into
taken taken into account.
account. In thisthe
In this study, study,
second thescheme
secondconsidering
scheme considering
the glacierthe glacier was
dynamics dy-
namics
used towas used to
simulate thesimulate the glacier process.
glacier process.
The
The mass-conserving
mass-conserving module of the the glacier
glacier divides
dividesthe theglacier
glacierinto
intoseveral
severalrectangular
rectangu-
lar grids according to the spatial resolution of the SPHY model,
grids according to the spatial resolution of the SPHY model, and divides the glacier and divides the glacier
into
into an accumulation area and ablation area based on the equilibrium
an accumulation area and ablation area based on the equilibrium line altitude (ELA) line altitude (ELA)of
of
thethe gridded
gridded glacier
glacier (Figure
(Figure 2). 2). When
When thethe glacier
glacier mass mass
lossloss of any
of any gridgrid
in thein ablation
the ablation
area
area exceeds
exceeds the thickness
the thickness of theofgrid
the grid glacier,
glacier, the glacier
the grid grid glacier will disappear,
will disappear, i.e.,
i.e., the the glac-
glacier will
ier will retreat,
retreat, and when andthewhentotalthe totalbalance
mass mass balance of the glacier
of the whole whole glacier for a hydrological
for a hydrological year is
year is negative,
negative, the accumulated
the accumulated mass of mass
the of the accumulation
accumulation areabewill
area will be equally
equally redistrib-
redistributed to
uted to the retreated
the retreated ablationablation
area at the areaend
at the endhydrological
of the of the hydrological
year. year.

Figure
Figure2.2.Division
Divisionof
ofaccumulation
accumulationarea
areaand
andablation
ablationarea
area of
of glacier.
glacier. Figure
Figure cited
cited from
from [66].
[66].

Dependingon
Depending onthe
thetemperature
temperatureand andaacritical
criticaltemperature
temperaturethreshold,
threshold,precipitation
precipitation can can
eitherfall
either fallas
as rainfall
rainfallor
or as
as snow
snow onon the
the glacier
glacier [47].
[47]. The
The rainfall
rainfall falling
falling on
on the
the glacier
glacier will
will
participate in
participate inthe
themelting
meltingprocess
process of
of the
theglacier
glacierin
inthe
theform
formofofmelting
meltingwater,
water, and
andthethesnow
snow
storageisisupdated
storage updatedwith withsnow
snow accumulation
accumulation and/or
and/or snowmelt.
snowmelt. The The melting
melting of glaciers
of glaciers con-
contributes to the total runoff by means of a slow and fast component, i.e., (i)
tributes to the total runoff by means of a slow and fast component, i.e., (i) percolation to percolation
to the
the groundwater
groundwater layer
layer that
that eventually
eventually becomes
becomes baseflow,and
baseflow, and(ii)
(ii)direct
directGR,
GR,respectively.
respectively.

2.4.3. Model Calibration and Validation

2.4.3. Model Calibration and Validation
The data source used for the First Chinese Glacier Inventory of the USR was aerial
The data source used for the First Chinese Glacier Inventory of the USR was aerial
photographs taken in the 1950s and 1960s, which reflect the distribution of glaciers in the
photographs taken in the 1950s and 1960s, which reflect the distribution of glaciers in the
USR before the 1970s. So, in this study, the year 1971 was used as the starting year for
USR before the 1970s. So, in this study, the year 1971 was used as the starting year for the
the simulation of runoff and runoff components, which minimized the simulation errors
simulation of runoff and runoff components, which minimized the simulation errors re-
resulting from inaccuracy of glacier inventory data. Meanwhile, this study set 1971 to 1990
sulting from inaccuracy of glacier inventory data. Meanwhile, this study set 1971 to 1990
as the calibration period and calibrated the SPHY model with the monthly runoff observed
at the outlet of the USR to obtain the sensitivity parameters (Table 2) for the SPHY model,
using these parameters for the runoff simulation from 1991 to 2020. Some of the sensitive
parameters of the SPHY model could be obtained from the relevant literature. In addition to
the calibration and validation of the simulated monthly runoff using the measured monthly
runoff, this study also calibrated and validated the simulated glacier area of the SPHY
model using the glacier area of the USR extracted from the First and Second Chinese Glacier
Water 2021, 13, 3357 7 of 21

Inventories and the multi-period Qilian Mountains Glacier Inventory. These parameters are
not interrelated and have different effects on different hydrological processes, which allows
them to be calibrated one by one. These parameters were calibrated using a stepwise
approximation method, i.e., we first assigned a maximum value to a parameter, then
a minimum value to that parameter, and then gradually decreased the range of values
assigned until the best NSE, R2 and RE were obtained.
To evaluate the performance of the SPHY model for simulating monthly runoff, the
Nash-Sutcliffe Efficiency Coefficient (NSE, Equation (4)), the coefficient of determination
(R2 , Equation (5)) and the relative error (RE, Equation (6)) were used in this study to
evaluate the efficiencies of the calibration and validation. If the monthly NSE > 0.5 and
monthly R2 > 0.6, the model performance was considered to be acceptable [67].
2
∑1n ( Qobs,i − Qsim,i )
NSE = 1 − 2
(4)
∑1n Qobs,i − n1 (∑1n Qobs,i )

2
∑1n


2 Qobs,i − Qobs,i Qsim,i − Qsim,i
R =
2
2 (5)
∑1n Qobs,i − Qobs ∑1n Qsim,i − Qsim
 n
∑1 Qsim,i

E= − 1 × 100% (6)
∑1n Qobs,i
where n is the time-step of monthly observed runoff and monthly simulated runoff; Qobs,i
(m3 /s) and Qsim,i (m3 /s) are the monthly observed runoff and monthly simulated runoff at
time step i, respectively; Qobs,i (m3 /s) and Qsim,i (m3 /s) are the mean monthly observed
runoff and mean monthly simulated runoff, respectively.

Table 2. Sensitive parameters of the SPHY model used for the simulation during the calibration
period and validation period.

Parameters Description Unit Value Literatures

Root depth Rootzone soil thickness mm 300
Sub depth Subzone soil thickness mm 1200
Temperature threshold for ◦C
Tcrit 0 [68]
precipitation to fall as snow
DDFS Degree-day factor for snow mm/◦ C/d 3 [69]
Degree-day factor for debris-free
DDFCI mm/◦ C/d 6 [69]
glaciers
δgw Groundwater recharge delay time d 3
αgw Baseflow recession factor - 0.01 [70]
BFtresh Threshold for baseflow to occur mm 0
kx Flow recession coefficient - 0.7

3. Results
3.1. Regional Climate Change and Runoff Evolution
The USR is a typical alpine and cold basin, which is extremely sensitive to climate
change and global warming. Under climate change and global warming, the USR has
experienced significant warming, while increases of precipitation and runoff have fluctu-
ated in the past few decades. Figure 3 shows the anomalies of annual mean temperature,
annual precipitation and annual runoff in the USR from 1961 to 2020 based on the statistics
of observed temperature, precipitation and runoff. The mean annual temperature in the
USR from 1961 to 2020 was approximately −1.62 ◦ C, and the temperature had a significant
increasing trend with an increase rate of 0.35 ◦ C/10a. However, there was an abrupt change
in the temperature in the USR. As shown in Figure 3a, the abrupt change in temperature
occurred around 1996, with a mean annual temperature of −2.12 ◦ C from 1961 to 1996
before the abrupt change and −0.88 ◦ C from 1997 to 2020 after the abrupt change. The
abrupt change of temperature in the USR was consistent with an abrupt change of temper-
 of observed temperature, precipitation and runoff. The mean annual temperature in the
USR from 1961 to 2020 was approximately −1.62 °C, and the temperature had a significant
increasing trend with an increase rate of 0.35 °C/10a. However, there was an abrupt
change in the temperature in the USR. As shown in Figure 3a, the abrupt change in tem-
Water 2021, 13, 3357 perature occurred around 1996, with a mean annual temperature of −2.12 °C from 1961 8 of to
21

1996 before the abrupt change and −0.88 °C from 1997 to 2020 after the abrupt change. The
abrupt change of temperature in the USR was consistent with an abrupt change of tem-
perature in Northwest
ature in Northwest China
China [33].[33].
The The annual
annual precipitation
precipitation in theinUSR
the USR also tended
also tended to in-
to increase
crease with the increase of temperature; the increase rate of the precipitation from
with the increase of temperature; the increase rate of the precipitation from 1961 to 2020 1961 to
2020 was 14.1 mm/10a. The mean annual precipitation from 1961 to 2020
was 14.1 mm/10a. The mean annual precipitation from 1961 to 2020 was approximately was approxi-
mately 235.2
235.2 mm, mm,
with with afluctuation
a large large fluctuation of annual
of annual precipitation
precipitation and a standard
and a standard devia-
deviation of
tion of 46.4
46.4 mm. mm.
The The change
abrupt abrupt of
change of precipitation
precipitation occurredoccurred
later thanlater
thatthan that of tempera-
of temperature, with
ture, with the
the abrupt abruptpoint
change change point at
at about aboutThe
2000. 2000. The annual
mean mean annual precipitation
precipitation was 218.4
was 218.4 mm
mm from 1961 to 2000 before the abrupt change and 268.8 mm from
from 1961 to 2000 before the abrupt change and 268.8 mm from 2001 to 2020 after2001 to 2020 after the
the
abrupt change.

Interannual variation of annual

Figure 3. Interannual annual mean
mean temperature,
temperature, annual precipitation and annual runoff
in the USR from 1961 to 2020.
2020. (a) Temperature anomaly; (b) Percentage of precipitation anomaly;
(c) Percentage
(c) Percentage of
of runoff
runoff anomaly.
anomaly.

8 3
The mean annual runoff of of the
the USR
USR from
from1961
1961to
to2020
2020was 10.54×× 10
was10.54 108 mm3, and the
the
8 3
increase rate
increase rate was 1.27 ×
was 1.27 × 10
108 mm3/10a.
/10a.The
Theabrupt
abruptchange
changein in runoff
runoff was
was more
more pronounced
pronounced
than that
than that of
of temperature
temperature and
and precipitation.
precipitation. AsAs shown
shown in
in Figure
Figure 3c,
3c, the abrupt change
the abrupt change point
point
of runoff occurred around 2000. The mean annual runoff from 1961–2000 before
of runoff occurred around 2000. The mean annual runoff from 1961–2000 before the ab- the abrupt
change
rupt was was
change × 10×8 10
8.90 8.90 m83 m
/10a, while
3/10a, whilethe
themean
meanannual
annualrunoff
runofffrom
from2001–2020
2001–2020 after
after the
the
abrupt change was 13.83 × 108 m3 /10a. The abrupt change characteristics of runoff were
more consistent with the abrupt change of precipitation.
The monthly variation of temperature and precipitation in the USR had characteristics
of an alpine climate and temperate continental climate. The seasonality of temperature and
precipitation were very obvious, with the maximum temperature in July and the minimum
temperature in January, while 91.4% of precipitation was mainly concentrated in the warm
season from May to September. The significant seasonal characteristics of temperature
and precipitation in the USR led to significant seasonal characteristics of runoff as well,
 abrupt change was 13.83 × 108 m3/10a. The abrupt change characteristics of runoff were
more consistent with the abrupt change of precipitation.
The monthly variation of temperature and precipitation in the USR had characteris-
tics of an alpine climate and temperate continental climate. The seasonality of temperature
Water 2021, 13, 3357
and precipitation were very obvious, with the maximum temperature in July and the9min- of 21
imum temperature in January, while 91.4% of precipitation was mainly concentrated in
the warm season from May to September. The significant seasonal characteristics of tem-
perature and precipitation in the USR led to significant seasonal characteristics of runoff
as
as shown
well, asinshown
Figurein4c.Figure
The maximum monthly monthly
4c. The maximum runoff occurred in July and
runoff occurred August,
in July and and
Au-
68.4% of the
gust, and annual
68.4% runoff
of the was concentrated
annual in June to September.
runoff was concentrated in June toThe maximumThe
September. monthly
maxi-
precipitation
mum monthly of precipitation
the USR was of in the
July,USR
while
wastheinmaximum
July, whilemonthly runoff was
the maximum in July
monthly and
runoff
August. This phenomenon indicates that the runoff of the USR has a hysteresis
was in July and August. This phenomenon indicates that the runoff of the USR has a hys- to the
precipitation, which reflects the water conservation capacity of the USR.
teresis to the precipitation, which reflects the water conservation capacity of the USR.

Figure 4.
Figure 4. Seasonal characteristics
characteristics of multi-year
multi-year monthly
monthly mean
mean temperature,
temperature, precipitation
precipitationand
andmean
mean
runoff in the USR. (a) Monthly mean temperature; (b) Monthly precipitation; (c) Monthly
runoff in the USR. (a) Monthly mean temperature; (b) Monthly precipitation; (c) Monthly mean mean
runoff.
runoff.

3.2.
3.2. Performance of the
the SPHY
SPHY Model
Model
Using
Using the
the same
same method
method asas that
that with
with which
which many
many hydrologic
hydrologic models
models arearecalibrated
calibrated
and
and validated,
validated, this
this study
study calibrated
calibrated andand validated
validated the
the SPHY
SPHY model
model with
with the
the observed
observed
monthly
monthlyrunoff
runoff[46,61].
[46,61].Due
Duetotothe
thelimitation of of
limitation thethe
time series
time of glacier
series inventory
of glacier data,data,
inventory the
simulation period of runoff and runoff components in the USR was set from
the simulation period of runoff and runoff components in the USR was set from 1971 to 1971 to 2020.
2020.Table 3 presents statistics on the evaluation metrics of the monthly runoff simulations
for the calibration period from 1971 to 1990, the validation period from 1991 to 2020 and
the entire simulation period from 1971 to 2020. As can be seen from Table 3, the NSE,
R2 and RE for the calibration period were 0.89, 0.90 and −3.82%, respectively, while the
NSE, R2 and RE for the validation period were 0.92, 0.92 and −11.51%, respectively. The
NSE and R2 of the SPHY model for the validation period were higher than those for the
calibration period, which may be related to improvements in the accuracy of temperature,
precipitation and glacier inventories over the validation period with advancements of
observation technology. However, the RE in the validation period was larger, which was
 bration period, which may be related to improvements in the accuracy of temp
for the calibration period from 1971 to 1990, the validation period from 1991 to 2020 and
precipitation
the entire simulation and1971
period from glacier inventories
to 2020. As can be over the validation
seen from Table 3, the period
NSE, R2 with advancemen
and RE for theservation
calibration technology.
period were 0.89,However, the RErespectively,
0.90 and −3.82%, in the validation period was larger, wh
while the NSE,
R2 and RE forduethe validation period were 0.92, 0.92 and −11.51%, respectively.
to the low simulation of monthly runoff in the cold season The NSE(Figure 5). The com
Water 2021, 13, 3357 and R2 of the SPHY model for the validation period were higher than those for the cali- 10 of 21
of the simulated and observed monthly runoff for the entire simulation period in
bration period, which may be related to improvements in the accuracy of temperature,
6 demonstrates that the simulated and observed monthly runoff were basically th
precipitation and glacier inventories over the validation period with advancements of ob-
with anHowever,
servation technology. NSE of the 0.91,RER
2 of 0.92 and RE of −6.44% indicating that the performanc
due to the low simulation ofinmonthly
the validation period
runoff in wasseason
the cold larger, which5).was
(Figure The comparison
due to the lowSPHY model
simulation of was
monthly acceptable.
runoff in the cold season (Figure 5). The comparison
of the simulated and observed monthly runoff for the entire simulation period in Figure
of the simulated and observedthat
6 demonstrates monthly runoff for
the simulated andthe entire simulation
observed periodwere
monthly runoff in Figure
basically the same,
6 demonstrates thatan
Table
with the simulated
3.NSE
Monthly and
R ofobserved
2
of 0.91,efficiencies formonthly
0.92 and runoff
calibration
RE of −6.44% were from
period basically
indicating 1971the
that tosame,
the 1990, validation
performance of theperiod fr
with an NSE of to 0.91, Rand
2020model
SPHY
2 of 0.92 and RE of −6.44% indicating that the performance of the
entire simulation period from 1971 to 2020.
was acceptable.
SPHY model was acceptable.
Evaluation
Table 3. Monthly Indicators
efficiencies forCalibration Period
calibration period from 1971 toValidation Period
1990, validation period fromSimulation
1991 to P
Table 3. Monthly efficiencies
2020 and entire for calibration
simulation period
period from
from 1971
1971 to to 1990,
2020. validation period from 1991
to 2020 and entire simulationNSEperiod from 1971 to 2020. 0.89 0.92 0.91
EvaluationR Indicators
2 0.90
Calibration Period 0.92 Simulation Period 0.92
Validation Period
Evaluation Indicators Calibration Period Validation Period Simulation Period
RE
NSE −3.82%
0.89 0.92 −11.51% 0.91 −6.44%
NSE 0.89 0.92 0.91
R2 0.90 0.92 0.92
R 2
RE 0.90 −3.82%0.92 −11.51%0.92 −6.44%
RE −3.82% −11.51% −6.44%

Figure 5. Simulated
Figure 5.and Figure
observed
Simulated and 5.observed
Simulated
monthly runoff and
for
monthly theobserved
for themonthly
calibration
runoff period runoff
from
calibration 1971 tofor
period the
1990
from calibration
and
1971 the1990
to andperiod
validation from
period 1971 to
from 1991 tothe validation period
2020. from 1991 to
the validation 2020. from 1991 to 2020.
period

Figure 6. Comparison of the simulated and observed monthly amounts of runoff for the entire sim-
ulation period from 1971 to 2020.

Figure
Figure 7 shows 6.6.Comparison
Figure the Comparison of the
monthly variation
of the ofsimulated andand
the simulated
simulated observed
and observed monthlymonthly
observed monthlyamounts
amounts ofrunoff of the
runoff for runoff for the en
entire
during the calibration
simulation
ulation and
periodvalidation
period periods.
from 1971
from 1971 toto
2020.It can be seen from Figure 7 that the SPHY
2020.
model could simulate the monthly variation of runoff better in both the calibration and
Figure 7 shows the monthly variation of the simulated and observed monthly runoff
duringFigure 7 showsand
the calibration thevalidation
monthly variation
periods. It canofbethe simulated
seen from Figure and observed
7 that the SPHYmonthly
during the calibration
model could and validation
simulate the monthly variationperiods. It caninbeboth
of runoff better seenthefrom Figure
calibration and 7 that th
validation
model periods.
could However,
simulate themonthly
the simulated variation
runoff in August was lower
of runoff thaninthe
better measured
both the calibrat
runoff in both the calibration and validation periods, which may have been caused by an
error of the SPHY model in simulating the recession process of the USR.
 Water 2021, 13, x FOR PEER REVIEW 11 of 21
021, 13, x FOR PEER REVIEW 11 of 21

validation periods. However, the simulated runoff in August was lower than the meas-
validation
Water 2021, 13, 3357 periods.
uredHowever, the simulated
runoff in both runoffand
the calibration in August wasperiods,
validation lower than
whichthemay
meas-
have been caused
11 of 21
ured runoff in both
by anthe calibration
error and model
of the SPHY validation periods, which
in simulating may have
the recession been of
process caused
the USR.
by an error of the SPHY model in simulating the recession process of the USR.

Figure 7. Comparison of the simulated and observed monthly mean amounts of runoff for the cali-
Figure 7. Comparison
Figure 7. Comparison of thebrationof period
simulatedtheand
simulated
from and
1971
observed toobserved
1990 and
monthly monthly mean
validation
mean amounts amounts
period from
of runoff of
therunoff
for1991 for(a)
to 2020. theperiod
calibration cali- from period;
Calibration 1971 (b)
to 1990 andbration period
validation from
period 19711991
Validation
from to 1990 and validation
period.
to 2020. period
(a) Calibration from 1991
period; to 2020. (a)
(b) Validation Calibration period; (b)
period.
Validation period.
AA comparison
comparison of of simulated
simulated and and observed
observed annual
annual runoff
runoffin in the
the USR
USR from from1971
1971 toto 2020
2020
A comparison
is of simulated
is shown in Figure and
Figure 8. observed
8.On Onthe annual
theannual runoff
annualscale,
scale, in
thethe the USR from
correlation
correlation 1971
between
between to 2020
thethe simulated
simulated andandob-
is shown in Figure
served 8. On
observed the annual
annual
annual runoff
runoff scale,
isisgood,thewith
good, withananRR2 2ofof
correlation between the simulated
0.89,which
0.89, which isis considered
consideredand toob-be acceptable.
acceptable. ItIt
served annual runoff
can
can be is good,
be seen
seen from
fromwith
Figure
Figurean R 8b2 of
8b that0.89,
the which is considered
variations
variations to be acceptable.
of the simulated
simulated and observed
observed It annual
annual runoff
runoff
can be seen fromwereFigure
were 8b thatthe
basically
basically thesame.
the variations
same. The of the simulated
The increase
increase rate
rate of and
of the
the observedannual
simulated
simulated annualrunoff
annual runoffwas
runoff was 1.3,
1.3, which
which
were basically was
the slightly
same. Thesmaller
increase than the
rate ofincrease
the rate
simulated of the observed
annual runoff
was slightly smaller than the increase rate of the observed annual runoff of 1.5.annual
was runoff
1.3, which of 1.5. In general,
it
it is
was slightly smalleris acceptable to
to use
than the increase
acceptable use therateSPHY
the of themodel
SPHY observed
model to
to simulate runoff
runoffofevolution
annual runoff
simulate in
in the
1.5. In general,
evolution the USR,
USR, with
with the
the
model
it is acceptable model accurately
to use the SPHY model
accurately reflecting
reflecting the
to simulatecharacteristics
runoff evolution
the characteristics of runoff evolution.
in the
of runoff USR, with
evolution. With
Withthe the help of
the help of the
the
model accurately SPHY
SPHY model,
reflecting
model,thewe can
can better
better understand
we characteristics of runoff
understand the variation
evolution.
the of
of runoff
variation With runoffthe and
and runoff
help of thecomponents
runoff components in in
SPHY model, we the
thecanUSR,
USR, which
better
which can
can help
understandhelp with
the planning,
withvariation
planning, ofmanagement
runoff and runoff
management and
and utilization
components
utilization of
of water
waterresources.
in resources.
the USR, which can help with planning, management and utilization of water resources.

Figure 8. Comparison
Figure 8. Comparison of the simulated
simulated and
and observed
observed annual
annual runoff
runofffrom
from1971
1971to
to2020.
2020. (a)
(a) Simulated
Simulated
Figure 8. Comparison
annual
annualof the simulated
runoff versus and observed
Observed annual
annual runoff
runoff; (b) from 1971
Interannual to 2020. (a)
variation Simulated
of observed runoff
versus Observed annual runoff; (b) Interannual variation of observed runoff and and
sim-
annual runoff versus Observed
ulated runoff.
simulated runoff.annual runoff; (b) Interannual variation of observed runoff and sim-
ulated runoff.
3.3. Glacier Shrinking in the USR
Glaciers are indicators of climate change and global warming. There is a significant
increase in temperature in the USR under climate change and global warming, which
directly leads to rapid melting and retreat of glaciers in the basin. Figure 9a compares
Water 2021, 13, x FOR PEER REVIEW 12 of 21

3.3. Glacier Shrinking in the USR

Water 2021, 13, 3357 12 of 21
Glaciers are indicators of climate change and global warming. There is a significant
increase in temperature in the USR under climate change and global warming, which di-
rectly leads to rapid melting and retreat of glaciers in the basin. Figure 9a compares the
simulated
the glacier
simulated areaarea
glacier of the
of SPHY model
the SPHY and the
model andextracted (observed)
the extracted glacier
(observed) area of
glacier the
area
USR
of thefrom
USR the
fromFirst
theand
FirstSecond Chinese
and Second GlacierGlacier
Chinese Inventories and the
Inventories multi-period
and Qilian
the multi-period
Mountains
Qilian GlacierGlacier
Mountains Inventory. It canItbe
Inventory. seen
can that the
be seen thatsimulated glacier
the simulated areaarea
glacier andand
the the
ob-
served glacier
observed area
glacier correlated
area correlatedwell, with
well, anan
with R2 Rof2 0.98, though
of 0.98, thoughthethe
SPHY
SPHY model more
model ac-
more
curately simulated
accurately simulated glacier
glaciershrinking
shrinkingininthe
theUSR.
USR.

Figure 9.
Figure 9. Comparison
Comparison of of the
the simulated
simulated glacier
glacier area
area and
and the
the extracted
extracted (observed)
(observed)glacier
glacierarea
areafrom
fromaa
multi-period glacier inventory, and the change of the simulated glacier area in the USR in
multi-period glacier inventory, and the change of the simulated glacier area in the USR in the past the past
50 years. (a) Simulated glacier area versus Observed glacier area; (b) Simulated glacier
50 years. (a) Simulated glacier area versus Observed glacier area; (b) Simulated glacier area from area from
1971 to 2020.
1971 to 2020.

Figure 9b
Figure 9b shows
shows thethe changes
changes of of glacier
glacier area
area inin the
the USR
USR from
from1971
1971toto2020
2020simulated
simulated
by the SPHY model. The glacier area in the USR decreased
by the SPHY model. The glacier area in the USR decreased from 473 km in 1971 tofrom 473 km 2 in2 1971 to 370

km2km
370 in22020, andand
in 2020, retreated
retreatedbyby103 km
103 km2, 2
with a aretreat
, with retreatratio
ratioofof21.8%
21.8%and
and aa retreat rate of
retreat rate of
2.06 km
2.06 km /a.
22 /a. With
With the
the increasing
increasing temperature
temperature in in the
the USR,
USR, the
the glacier
glacier accelerated
accelerated itsits melting
melting
and retreated,
and retreated,which
whichdirectly
directlyled
ledfirst
firstto
toananincrease
increaseand andthen
thenaadecrease
decreaseofofGR,
GR,weakening
weakening
the recharge
the recharge andand regulating
regulatingeffect
effectofofthe
theglaciers
glacierson onriver
riverrunoff.
runoff.IfIfthe
theglaciers
glaciersin inthe
theUSR
USR
continue to retreat at a rate of 2.06 km2 2/a, the glacier area in the basin will shrink to 43%
continue to retreat at a rate of 2.06 km /a, the glacier area in the basin will shrink to 43% of
of that
that in the
in the 1970s
1970s by by
thethe
endend of this
of this century.
century.

3.4.
3.4. Annual
Annual Variations
Variations of
of Runoff
Runoff Components
Components
According
According to the structure of
to the structure of the
the SPHY
SPHY model,
model, the
the total
total runoff
runoff isiscomposed
composed of of four
four
runoff components: GR, snow runoff, rainfall runoff and baseflow.
runoff components: GR, snow runoff, rainfall runoff and baseflow. Among them, Among them, the base-
the
flow is mainly
baseflow generated
is mainly by recharge
generated of groundwater
by recharge from rainfall
of groundwater and a and
from rainfall smalla portion
small por- of
glacier
tion of meltwater. The fourThe
glacier meltwater. runoff
fourcomponents have different
runoff components responses
have different to climate
responses to change
climate
and global
change andwarming, resulting resulting
global warming, in large interannual variation ofvariation
in large interannual each runoff component.
of each runoff com-The
interannual variations of the four runoff components are shown in Figure
ponent. The interannual variations of the four runoff components are shown in Figure 10. 10. According to
this, it can be seen that from 1971 to 2020, the GR, rainfall runoff and
According to this, it can be seen that from 1971 to 2020, the GR, rainfall runoff and baseflow showed a
significant increasing trend. The increase rate of rainfall runoff was the largest, followed
by the GR and baseflow, with increase rates of 0.56 × 108 m3 /10a, 0.43 × 108 m3 /10a and
0.21 × 108 m3 /10a, respectively. The decreasing trend of snow runoff was not significant
due to global warming, which led to a decrease in snowfall, resulting in a decrease in
snow meltwater. In the last two decades, snowfall also increased to some extent due to
 precipitation was more in the form of rainfall, and the increase in precipitation did not
slow down the decrease in snowfall due to global warming.
It can also be seen from Figure 10 that there were large interannual fluctuations of
GR and rainfall runoff, while the interannual fluctuations of the baseflow were very small,
indicating that the baseflow was the most stable component of the total runoff. The large
Water 2021, 13, 3357 13 of 21
interannual fluctuations of GR were related to the interannual fluctuations of temperature
and precipitation. In the years with more precipitation, the rainfall runoff was high and
the precipitation had a certain inhibitory effect on glacier melting, while in the years with
the the
less precipitation, increase
high in precipitation.
temperatures andHowever,
increasedthe trends
glacier of snow
melting ledrunoff, rainfall runoff and the
to an increase
of GR, which hadinterannual
a rechargevariation of precipitation
and regulating in Figure
effect on total 10 Interannual
runoff. show that the increase in precipitation
fluctuations
was
of rainfall runoff more
were in the related
directly form oftorainfall, and thedue
precipitation increase in fluctuations
to large precipitationindid not slow down the
annual
decrease in snowfall
precipitation in Northwest China. due to global warming.

Figure 10. Interannual

Figure 10. Interannual variation
variation of simulatedof simulated runoff components.
runoff components. (a)runoff
(a) Glacier Glacier runoff
(GR); (b)(GR);
Snow(b) Snow(c) Rainfall
runoff;
runoff; (c) Rainfall runoff; (d) Baseflow.
runoff; (d) Baseflow.

3.5. Monthly VariationsIt can also be

of Runoff seen from Figure 10 that there were large interannual fluctuations of
Components
GR and rainfall runoff, while the interannual
Due to significant seasonal characteristics fluctuations
of temperature of the baseflow
and precipitation were very small,
in the
indicating that the baseflow was the most stable component
USR and different responses of the four runoff components to temperature and precipita- of the total runoff. The large
interannual fluctuations of GR were related to the interannual
tion, the monthly variations among the four runoff components were very different. It canfluctuations of temperature
and precipitation.
be seen from Figure In the years
11 that GR mainly withfrom
occurred moreJune
precipitation, the rainfall
to September, runoff
with peaks in was high and
July and August, the which
precipitation had a certain
were consistent withinhibitory effect
the highest on glacier melting,
temperatures while in
and amounts of the years with
less precipitation, the high temperatures and increased glacier
precipitation in July and August. Snow runoff mainly occurred from April to June, with melting led to an increase of
GR, which had a recharge and regulating effect on total runoff. Interannual fluctuations
May being the largest month for snow runoff. This seasonal characteristic of snow runoff
of rainfall runoff were directly related to precipitation due to large fluctuations in annual
was strongly influenced by temperature changes during the year, with April being the
precipitation in Northwest China.
time when the temperature gradually rose from below 0 °C in the USR, resulting in snow
cover starting to melt
3.5. during
Monthly this time.
Variations of Snow
Runoffrunoff could occur between March and Oc-
Components
tober, which was related to the topography of the USR. The average elevation of the USR
Due to significant seasonal characteristics of temperature and precipitation in the USR
is close to 4000 m, and the absolute height difference between the mountains on both sides
and different responses of the four runoff components to temperature and precipitation, the
monthly variations among the four runoff components were very different. It can be seen
from Figure 11 that GR mainly occurred from June to September, with peaks in July and
August, which were consistent with the highest temperatures and amounts of precipitation
in July and August. Snow runoff mainly occurred from April to June, with May being the
largest month for snow runoff. This seasonal characteristic of snow runoff was strongly
influenced by temperature changes during the year, with April being the time when the
temperature gradually rose from below 0 ◦ C in the USR, resulting in snow cover starting
to melt during this time. Snow runoff could occur between March and October, which
was related to the topography of the USR. The average elevation of the USR is close to
4000 m, and the absolute height difference between the mountains on both sides of the
basin is large, so even in the highest temperature season, snowfall can still occur at high
elevations, which is why snow runoff can occur from March to October. The seasonal
characteristics of rainfall runoff were clearly consistent with the seasonal characteristics of
precipitation, and the months with the highest precipitation were also the months with the
021, 13, x FOR PEER REVIEW 14 of 21

of the basin is large, so even in the highest temperature season, snowfall can still occur at
Water 2021, 13, 3357 14 of 21
high elevations, which is why snow runoff can occur from March to October. The seasonal
characteristics of rainfall runoff were clearly consistent with the seasonal characteristics
of precipitation, and the months with the highest precipitation were also the months with
the highest rainfall runoff.
highest Therunoff.
rainfall monthly variations
The monthlyof the baseflow
variations were
of the very smooth,
baseflow were veryandsmooth, and the
the baseflow only increased
baseflow only somewhat in the month
increased somewhat of month
in the maximum rainfall runoff,
of maximum rainfallwhich
runoff, which gives
gives a good indication that the baseflow
a good indication was the most
that the baseflow basic
was the andbasic
most stable component
and of the of the inland
stable component
inland rivers inrivers
Northwest China. China.
in Northwest

Figure 11. 11.

Figure Monthly variation
Monthly of simulated
variation runoff
of simulated components.
runoff (a) 1971~2000;
components. (b) 2001~2020.
(a) 1971~2000; (b) 2001~2020.

Figure 11a shows Figurethe 11a showsvariations

monthly the monthly variations
of the four runoffof the four runoffbefore
components components
the before the
abrupt change abrupt
of totalchange of total
runoff from runoff
1971 fromand
to 2000, 1971Figure
to 2000,
11band Figure
shows the11b showsvari-
monthly the monthly varia-
tions
ations of the four of the
runoff four runoffafter
components components
the abruptafter the abrupt
change of totalchange
runoff of total
from runoff
2001 to from 2001 to
2020. It can be seen from these two figures that the monthly mean runoff in the warm in the warm
2020. It can be seen from these two figures that the monthly mean runoff
season increasedseason increased
significantly significantly
after the abruptafter the abrupt
change of totalchange
runoff, of total
and therunoff,
increaseandof the increase of
the monthly meanthe monthly
runoff inmean runoffseason
the warm in thewas
warm season
mainly was mainly
caused caused by
by the increase ofthe
GRincrease of GR
and rainfall
and rainfall runoff, runoff, which
which indicates that theindicates
increasethat the runoff
of total increasewasofmainly
total runoff
due towas
the mainly due to
increase of GR and rainfall runoff. It can also be seen from these two figures that the oc-figures that the
the increase of GR and rainfall runoff. It can also be seen from these two
currence of snowoccurrence
runoff had of snow runoff
a certain hadadvance,
timing a certainindicating
timing advance, indicating
the objective thein-
fact of objective fact of
creasing temperature. That advance may have caused snowmelt floods, which can lead towhich can lead
increasing temperature. That advance may have caused snowmelt floods,
flood disasters to
andflood disasters
personal andand personal
property and property damage.
damage.
4. Discussion
4. Discussion
Water 2021, 13, x FOR PEER REVIEW 15 of 21
4.1. Contributions of Runoff Components to Total Runoff
4.1. Contributions of Runoff Components to Total Runoff
Figure 12 shows the mean annual contributions of the four runoff components to
Figure 12 shows therunoff
the total mean annual contributions
from 1971 to 2020. Itof can
the four runoff
be seen components
from Figure 12tothat
the the GR, snow
total runoff from
Yet, 1971
runoff, to 2020.
withrainfall
climate It can
change
runoff be
andandseen from
global
baseflow Figure 12
warming that 16%,
leading
contributed 28%, the GR,
18%snow
to rapid meltingrunoff,
and 38% and retreat
to the totalofrunoff,
glac-
rainfall runoffiers
and baseflow
and contributed
reduction of snow 28%, 16%,
cover, the18% and
recharge 38%
and to the total
regulation runoff,
effect ofre-
river
respectively. Baseflow was the runoff component that contributed the most to total runoff, runoff by
spectively. Baseflow
glacierswas
and the runoff
snow component
cover meltwater that
is contributed
weakening.
followed by GR, rainfall runoff and snow runoff. the most to total runoff,
followed by GR, rainfall runoff and snow runoff.
According to Figure 3c, a significant abrupt change in the runoff of the USR occurred
at the end of the 20th century. Figure 13 compares the contributions of each runoff com-
ponent to the total runoff in the two periods, i.e., before and after the abrupt change of
runoff in the USR. As can be seen from Figure 13 that the contributions of snow runoff
and the baseflow to the total runoff decreased significantly, while the contribution of rain-
fall runoff to the total runoff increased significantly. Meanwhile, the contribution of GR
to the total runoff did not decrease significantly, but the contribution range of GR to the
total runoff decreased significantly in the period after the abrupt change of runoff from
2001 to 2020. The total mean annual contribution of GR and snow runoff to the total runoff
was 44%, indicating that the recharge of river runoff in the USR was not dominated by
glaciers and snow cover meltwater, but the role of glaciers and snow cover meltwater in
recharging and regulating river runoff was very important. In the years with less precip-
itation, the recharge and regulation of river runoff by glacier meltwater was prominent.

Figure
Figure12.
12.Mean
Meanannual
annualcontributions
contributionsof
ofrunoff
runoffcomponents
componentsto
tothe
thetotal
totalrunoff.
runoff.

According to Figure 3c, a significant abrupt change in the runoff of the USR occurred at
the end of the 20th century. Figure 13 compares the contributions of each runoff component
Water 2021, 13, x FOR PEER REVIEW 15 of 21

Yet, with climate change and global warming leading to rapid melting and retreat of glac-
Water 2021, 13, 3357 iers and reduction of snow cover, the recharge and regulation effect of river runoff
15 ofby
21
glaciers and snow cover meltwater is weakening.

to the total runoff in the two periods, i.e., before and after the abrupt change of runoff
in the USR. As can be seen from Figure 13 that the contributions of snow runoff and the
baseflow to the total runoff decreased significantly, while the contribution of rainfall runoff
to the total runoff increased significantly. Meanwhile, the contribution of GR to the total
runoff did not decrease significantly, but the contribution range of GR to the total runoff
decreased significantly in the period after the abrupt change of runoff from 2001 to 2020.
The total mean annual contribution of GR and snow runoff to the total runoff was 44%,
indicating that the recharge of river runoff in the USR was not dominated by glaciers and
snow cover meltwater, but the role of glaciers and snow cover meltwater in recharging
and regulating river runoff was very important. In the years with less precipitation, the
recharge and regulation of river runoff by glacier meltwater was prominent. Yet, with
climate change and global warming leading to rapid melting and retreat of glaciers and
reduction of snow cover, the recharge and regulation effect of river runoff by glaciers and
snow cover
Figure meltwater
12. Mean is weakening.
annual contributions of runoff components to the total runoff.

Figure 13. Comparison

Comparison ofof contributions
contributions of
of runoff components
components to total runoff before and after the
abrupt
abrupt change
change of
of the total runoff.
the total runoff.

4.2. Increasing
4.2. Increasing Total
Total Runoff
Runoff Dominated
Dominated by by Increasing
Increasing Rainfall
Rainfall
Figure 14 shows the interannual variation of thethe
Figure 14 shows the interannual variation of contributions
contributions offour
of the the runoff
four runoff
com-
components to the total runoff and the five-year fast Fourier filter of the contributions.
ponents to the total runoff and the five-year fast Fourier filter of the contributions. As can As
can be seen from the figure, the contribution of GR to the total runoff
be seen from the figure, the contribution of GR to the total runoff showed a trend of firstshowed a trend of
first increasing and then decreasing. Although the interannual variation
increasing and then decreasing. Although the interannual variation of snow runoff did of snow runoff
did show
not not show a significant
a significant decreasing
decreasing trend,trend, the contribution
the contribution of snowof runoff
snow runoff to therunoff
to the total total
runoff did. The decreasing trend of the contribution of GR and snow
did. The decreasing trend of the contribution of GR and snow runoff to the total runoff runoff to the total
runoff indicates
indicates that thethat the alimentation
alimentation of theof the glaciers
glaciers and snow
and snow covercover
to thetoriver
the river
runoffrunoff
waswas
de-
decreasing. The five-year fast Fourier filter of the contribution of rainfall
creasing. The five-year fast Fourier filter of the contribution of rainfall runoff indicates runoff indicates
that the
that the contribution
contribution of of rainfall
rainfall runoff
runoff toto the
the total
total runoff
runoff showed
showed aa significant
significant increasing
increasing
trend, consistent with the increasing trend of rainfall runoff. Although it can be seen from
trend, consistent with the increasing trend of rainfall runoff. Although it can be seen from
Figure 10d that the baseflow had a significant increasing trend, the contribution of the
Figure 10d that the baseflow had a significant increasing trend, the contribution of the
baseflow to the total runoff showed a significant decreasing trend, meaning the trends of the
baseflow to the total runoff showed a significant decreasing trend, meaning the trends of
baseflow and contribution of the baseflow were the opposite. Moreover, according to the
the baseflow and contribution of the baseflow were the opposite. Moreover, according to
structure of the SPHY model, the water source of the baseflow is mainly the percolation of
the structure of the SPHY model, the water source of the baseflow is mainly the percola-
rainfall and a small part made up by glacier meltwater, meaning an increase of the baseflow
tion of rainfall and a small part made up by glacier meltwater, meaning an increase of the
actually depends on an increase of rainfall. The decreasing trends of the contributions of
baseflow actually depends on an increase of rainfall. The decreasing trends of the contri-
GR, snow runoff and the baseflow to the total runoff, and the increasing trend of rainfall
butions of GR, snow runoff and the baseflow to the total runoff, and the increasing trend
runoff in contributing to the total runoff fully indicate that with climate change and global
of rainfall runoff in contributing to the total runoff fully indicate that with climate change
warming, the rapid melting of glaciers in the past decades was an important reason for the
and global warming, the rapid melting of glaciers in the past decades was an important
increase in runoff. However, with the rapid retreat of glaciers, the contribution of GR to
the increase in the total runoff peaked and gradually decreased, and the contribution of
increased rainfall to the increase in the total runoff exceeded the contribution of GR to the
increase in runoff. This means that a future increase in the total runoff will depend on an
increase in rainfall in the USR.
 reason for the increase in runoff. However, with the rapid retreat of glaciers, the contribu-
tion of GR to the increase in the total runoff peaked and gradually decreased, and the
Water 2021, 13, 3357contribution of increased rainfall to the increase in the total runoff exceeded the contribu- 16 of 21
tion of GR to the increase in runoff. This means that a future increase in the total runoff
will depend on an increase in rainfall in the USR.

Figure 14. Interannual

Figure 14. Interannual variation ofvariation of contributions
contributions of runoff components
of runoff components to total
to total runoff. runoff. (a) Contribu-
(a) Contribution of GR; (b) Contribu-
tion of snowtion of GR;
runoff; (c)(b) Contribution
Contribution of snowrunoff;
of rainfall runoff;(d)
(c)Contribution
Contributionofofbaseflow.
rainfall runoff; (d) Contribution of
baseflow.
However, increases in the rainfall in Northwest China are mainly dominated by
However,short-duration
increases in the rainfall
heavy in Northwest
rainfall episodes, suchChina are mainly
as extreme dominated
rainfall by
and short-term convective
short-duration rainfall,
heavy rainfall
whichepisodes,
may leadsuch as extreme
to disasters suchrainfall and short-term
as flooding and causeconvective
serious damage to flood
rainfall, which control
may lead to disasters
facilities. such asthe
In addition, flooding
rainfalland cause serious
in Northwest damage
China occurstomostly
flood in the form of
control facilities. In addition,
pulses and canthe rainfall in Northwest
be extremely China
irregular. Years occurs
with less mostly
rainfallin thelead
can form
to of
droughts, which
pulses and can in beturn, can have
extremely important
irregular. Years impacts on the
with less lives can
rainfall of residents and industrial
lead to droughts, which and agricultural
in turn, can haveproduction.
importantInimpacts
conclusion,
on theunder
livesclimate change
of residents andand global warming,
industrial the role of glaciers
and agricul-
tural production.andInsnow cover inunder
conclusion, recharging
climate and regulating
change the runoff
and global of thethe
warming, USR is of
role weakening, and
glaciers and snowthe cover
influence of rainfalland
in recharging on the evolution
regulating the of runoff
runoff is becoming
of the increasingly prominent.
USR is weakening,
and the influenceHowever, dueon
of rainfall to the significant
evolution of pulsating
runoff ischaracteristics of regionalpromi-
becoming increasingly rainfall, uncertainty
around future changes in runoff is increasing, along with
nent. However, due to the significant pulsating characteristics of regional rainfall, uncer-the impact of disasters such as
tainty around future changes in runoff is increasing, along with the impact of disasters management
floods and droughts. This creates new requirements for scientific planning,
and sustainable
such as floods and droughts. This usecreates
of the limited water resources
new requirements in Northwest
for scientific China.
planning, man-
agement and sustainable use of the limited water resources in Northwest China.
4.3. Uncertainty in the Simulation of Runoff and Runoff Components
Uncertainty
4.3. Uncertainty in the Simulation in the simulation
of Runoff of runoff
and Runoff and runoff components in the USR using the
Components
SPHY model surrounds the accuracy of the forcing data obtained and the rationality of its
Uncertainty in the simulation of runoff and runoff components in the USR using the
reprocessing, along with the setting of the sensitivity parameters of the SPHY model.
SPHY model surrounds the accuracy of the forcing data obtained and the rationality of its
The SPHY model is a spatially distributed model, and the temperature and precipita-
reprocessing, along with the setting of the sensitivity parameters of the SPHY model.
tion data required to drive it are also distributed [47], whereas the observed temperature
The SPHY model is a spatially distributed model, and the temperature and precipi-
and precipitation are point-based, so the temperature and precipitation need to be interpo-
tation data required to drive it are also distributed [47], whereas the observed temperature
lated across the basin. The interpolation method for temperature and precipitation is based
and precipitation are point-based, so the temperature and precipitation need to be inter-
on bilinear interpolation, which relies on the altitude gradient of temperature and precip-
polated across the basin. The interpolation method for temperature and precipitation is
itation, thus inevitably increasing the uncertainty in the simulation of runoff and runoff
based on bilinear interpolation, which relies on the altitude gradient of temperature and
components. Improving the accuracy of the interpolation of temperature and precipitation
precipitation, thus inevitably increasing the uncertainty in the simulation of runoff and
presents a major difficulty in hydrological modeling, which should be pursued in future
runoff components. Improving the accuracy of the interpolation of temperature and pre-
work to find more accurate methods for the interpolation of temperature and precipitation.
cipitation presents a majorstudies
However, difficulty in hydrological
based modeling,
on temperature which should
and precipitation be pursued show that the
observations
temperature and precipitation in the Qilian Mountains are altitude-dependent [71,72], with
a linear increase in temperature and precipitation with increasing altitude [73], and the
altitude-dependency of precipitation reduces the uncertainty of precipitation interpolation.
Other data such as DEM and land use are widely used and recognized data with a high
level of accuracy. The simulation of glacier thickness in this study was carried out using
Water 2021, 13, 3357 17 of 21

the Glabtop2 method, which is based on DEM and glacier outlines and has been used
with a high accuracy in the simulation of glacier thickness in other regions [56–58]. The
applicability of the Glabtop2 method in the Upper Shule River basin needs to be validated
due to the lack of glacier thickness data, which is dependent on the detection of glacier
thicknesses in this basin.
The parameters of the SPHY model are also one of the reasons for uncertainty. Some
of the sensitive parameters can be obtained from the relevant literature, while others need
to be calibrated. The sensitive parameters of the SPHY model listed in Table 2 actually have
a large spatial variation, and the use of averages is a common way of setting parameters
in hydrological simulations when a significant number of targeted parameters are not
investigated. Previous studies have shown that there is a large spatial variability in the
critical temperature of rain and snow separation in the Qilian Mountains [12]. The critical
temperature of rain and snow separation fluctuates widely, and there is no relatively fixed
value of the critical temperature. Similarly, the degree-day factor of glaciers and snow
cover also has large spatial variability [74], and the spatial variability of parameters is
difficult to obtain in the high mountainous basin, which increases the uncertainty of the
simulation to a certain extent. Admittedly, the method of setting parameters involves great
uncertainty, which should be investigated in future work with a large number of sensitive
parameters to improve the accuracy of runoff and runoff component simulations.
In addition, monthly mean runoff data and multi-period glacier inventory data with
measurement errors were used to validate the performance of the SPHY model, which may
have increased the uncertainty surrounding the simulation. The runoff was measured at the
National Hydrological Station of Changmabao and the measurement errors in runoff were
within specified error limits. Although the construction of reservoirs near the outlet of the
USR in recent years has affected the natural runoff involution, the impact on the total water
resources is not significant, and the errors in measured runoff are considered acceptable.
There is a slight difference in albedo between glaciers and snow cover, and the presence of
snow inevitably affects glacier inventorying [17]. Therefore, this study corrected the glacier
inventory data based on remote-sensing images without snow cover in non-glacierized
areas in the ablation season, which ensured the accuracy of the validation data.

5. Conclusions
Based on the observed hydrometeorological data and multi-period glacier inventory
data, the contributions of runoff components to the total runoff, along with the essential
driving factors of runoff evolution, were quantified in the Upper Shule River. To do so, we
used a distributed cryospheric-hydrological SPHY model coupled with a mass-conserving
glacier module. Comparisons between the monthly observed and simulated runoff and
between changes in the observed and simulated areas of glacier show indicated that the
performance of the SPHY model was reasonable.
The glaciers continually retreated with the increasing of temperature, and the area of
glacier coverage decreased from 473 km2 to 370 km2 between 1971 and 2020, a decrease
of approximately 21.8%, with a reduction rate of 2.06 km2 /a. The accelerated melting of
glaciers led to an increasing trend of glacier runoff, and the increasing precipitation led to
an increasing trend of rainfall runoff and baseflow. Meanwhile, the decreasing trend of
snow runoff was not significant, indicating that the consistent increase in the total runoff
was due to increasing rainfall and glacier melt over the past few decades.
The baseflow is the most stable and greatest contributing component to the total
runoff, with a mean annual contribution of 38%, while the glacier runoff, rainfall runoff
and snow runoff have greater interannual fluctuations, with mean annual contributions of
28%, 18% and 16% to total runoff, respectively. With decreasing glacier area and increasing
temperature, the annual contribution of glacier runoff and snow runoff to the total runoff
showed a decreasing trend.
Although the baseflow was increasing, the contribution of the baseflow to the total
runoff was not increasing, indicating that the increases in the total runoff depended more
Water 2021, 13, 3357 18 of 21

on increases in rainfall runoff dominated by surface runoff and lateral flow in shallow soil.
The role of glaciers in regulating and recharging runoff was diminishing, which indicates
that the future evolution of runoff in the Upper Shule River will depend on rainfall with
pulsating characteristics. This places new demands on the construction of flood- and
drought-control facilities in the basin.

Author Contributions: J.W. set up the SPHY model, simulated the runoff, analyzed the data and
wrote the manuscript; H.L. prepared the driven data and reviewed the manuscript; J.Z., S.T. and
X.W. conducted data collection. All authors have read and agreed to the published version of the
manuscript.
Funding: This study was supported by the National Natural Science Foundation of China (NSFC)
(41771084, 41730751).
Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.
Data Availability Statement: Data sharing not applicable.
Acknowledgments: The authors would like to thank Qiudong Zhao and Jia Qin from the Key
Laboratory of Ecohydrology of the Inland River Basin, Northwest Institute of Eco-Environment and
Resources, Chinese Academy of Sciences, Lanzhou, China for their help with setting up the SPHY
model and collecting data.
Conflicts of Interest: The authors declare no conflict of interest.

References
1. Li, K.; Chen, R.; Liu, G. Cryosphere water resources simulation and service function evaluation in the Shiyang River Basin of
Northwest China. Water 2021, 13, 114. [CrossRef]
2. Ding, Y.; Ye, B.; Han, T.; Shen, Y.; Liu, S. Regional difference of annual precipitation and discharge variation over west China
during the last 50 years. Sci. China Ser. D Earth Sci. 2007, 50, 936–945. [CrossRef]
3. Yang, Z. Glacier meltwater runoff in China and its nourishment to river. Chin. Geogr. Sci. 1995, 5, 66–76. [CrossRef]
4. Yang, Z. Mountain stream types in Northwest China. J. Glaciol. Geocryol. 1981, 3, 24–31.
5. Yang, Z. Glacier water resources of Qilian Mountains. J. Glaciol. Geocryol. 1988, 10, 36–46.
6. Chen, Y.; Yang, Q.; Luo, Y.; Shen, Y.; Pan, X.; Li, L.; Li, Z. Ponder on the issues of water resources in the arid region of northwest
China. Arid. Land Geogr. 2012, 35, 1–9.
7. Zhang, S.; Ye, B.; Liu, S.; Zhang, X.; Hagemann, S. A modified monthly degree-day model for evaluating glacier runoff changes in
China. Part I: Model development. Hydrol. Process. 2012, 26, 1686–1696. [CrossRef]
8. Piao, S.; Ciais, P.; Huang, Y.; Shen, Z.; Peng, S.; Li, J.; Zhou, L.; Liu, H.; Ma, Y.; Ding, Y.; et al. The impacts of climate change on
water resources and agriculture in China. Nature 2010, 467, 43–51. [CrossRef] [PubMed]
9. Liu, S.; Ding, Y.; Li, J.; Shangguan, D.; Zhang, Y. Glaciers in response to recent climate warming in western China. Quat. Sci. 2006,
26, 762–771.
10. Cai, X.; Li, Z.; Zhang, H.; Xu, C. Vulnerability of glacier change in the Tianshan Mountains region of China. J. Geogr. Sci. 2021, 31,
1469–1489. [CrossRef]
11. Duan, J.; Wang, L.; Ren, J.; Li, L. Progress in glacier variations in China and its sensitivity to climatic change during the past
century. Prog. Geogr. 2009, 28, 231–237.
12. Liu, G.; Chen, K.; Li, K. Glacial change and its hydrological response in three inland river basins in the Qilian Mountains, Western
China. Water 2021, 13, 2213.
13. Ding, Y.; Liu, S.; Li, J.; Shangguan, D. The retreat of glaciers in response to recent climate warming in western China. Ann. Glaciol.
2006, 43, 97–105. [CrossRef]
14. Zhang, X.; Chen, R.; Liu, G.; Yang, Y.; Feng, T. Economic value of freshwater provisioning services of the cryosphere in the
Urumqi River, Northwest China. Adv. Clim. Chang. Res. 2021, in press. Available online: https://www.sciencedirect.com/
science/article/pii/S1674927821001295 (accessed on 20 November 2021). [CrossRef]
15. Ding, Y.; Zhao, Q.; Wu, J.; Zhang, S. The future changes of Chinese cryospheric hydrology and their impacts on water security in
arid areas. J. Glaciol. Geocryol. 2020, 42, 23–32.
16. Sun, M.; Liu, S.; Yao, X.; Guo, W.; Xu, J. Glacier changes in the Qilian Mountains in the past half-century: Based on the revised
first and second Chinese glacier inventory. J. Geogr. Sci. 2018, 28, 206–220. [CrossRef]
17. Wang, H.; Jia, B.; Wang, Y.; Hou, S. Glacier changes in Shulenan Mountain from 1973 to 2018. J. Arid. Land Resour. Environ. 2021,
35, 60–65.
Water 2021, 13, 3357 19 of 21

18. Wang, Y.; Li, J.; Wu, L.; Guo, L.; Li, J. Using remote sensing images to monitor the glacier changes in Qilian Mountains during
1987–2018 and analyzing the impact factors. J. Glaciol. Geocryol. 2020, 42, 344–356.
19. Shi, Y.; Cao, X.; Wu, J.; Filippo, G.; Dong, W. Simulation of the changes in snow cover over China under global warming by a high
resolution RCM. J. Glaciol. Geocryol. 2010, 32, 215–222.
20. Che, T.; Hao, X.; Dai, L.; Li, H.; Huang, X.; Xiao, L. Snow cover variation and its impacts over the Qinghai-Tibet Plateau. Bull.
Chin. Acad. Sci. 2019, 34, 1247–1253.
21. Liu, X.; Chen, R.; Liu, J.; Wang, X.; Wang, L. Effects of snow-depth change on spring runoff in cryosphere areas of China. Hydrol.
Sci. J. 2019, 64, 789–797. [CrossRef]
22. Ma, R. Characteristics and Causes of Snow Change in Arid Northwest China from 1979 to 2016; Northwest Normal University: Xi’an,
China, 2018.
23. Wang, C.; Zhang, S.; Li, K.; Zhang, F.; Yang, K. Change characteristics of precipitation in Northwest China from 1961 to 2018.
Chin. J. Atmos. Sci. 2021, 45, 713–724.
24. Li, M.; Sun, H.; Su, Z. Research progress in dry/wet climate variation in Northwest China. Geogr. Res. 2021, 40, 1180–1194.
25. Wang, L.; Chen, R.; Han, C.; Wang, X.; Liu, G.; Song, Y.; Yang, Y.; Liu, J.; Liu, Z.; Liu, X.; et al. Change characteristics of
precipitation and temperature in the Qilian Mountains and Hexi Oasis, Northwestern China. Environ. Earth Sci. 2019, 78, 284.
[CrossRef]
26. Li, B.; Chen, Y.; Shi, X.; Chen, Z.; Li, W. Temperature and precipitation changes in different environments in the arid region of
northwest China. Theor. Appl. Climatol. 2012, 112, 589–596. [CrossRef]
27. Wang, C.; Zhang, S.; Zhang, F.; Li, K.; Yang, K. On the increase of precipitation in the Northwestern China under the global
warming. Adv. Earth Sci. 2021, 36, 980–989.
28. Zhang, Q.; Zhu, B.; Yang, J.; Ma, P.; Liu, X.; Lu, G.; Wang, Y.; Yu, H.; Liu, W.; Wang, D. New characteristics about the climate
humidification trendin Northwest China. Chin. Sci. Bull. 2021, 66, 3757–3771. [CrossRef]
29. Ma, S.; Zhou, T.; Dai, A.; Han, Z. Observed changes in the distributions of daily precipitation frequency and amount over China
from 1960 to 2013. J. Clim. 2015, 28, 6960–6978. [CrossRef]
30. Chang, M.; Liu, B.; Martinez-Villalobos, C.; Ren, G.; Li, S.; Zhou, T. Changes in extreme precipitation accumulations during the
warm season over Continental China. J. Clim. 2020, 33, 10799–10811. [CrossRef]
31. Zhang, W.; Zhou, T. Significant increases in extreme precipitation and the associations with global warming over the global land
Monsoon regions. J. Clim. 2019, 32, 8465–8488. [CrossRef]
32. Zhou, B.; Xu, Y.; Wu, J.; Dong, S.; Shi, Y. Changes in temperature and precipitation extreme indices over China: Analysis of a
high-resolution grid dataset. Int. J. Climatol. 2016, 36, 1051–1066. [CrossRef]
33. Chen, R.; Shen, Y.; Mao, W.; Zhang, S.; Lv, H.; Liu, Y.; Liu, Z.; Fang, S.; Zhang, W.; Chen, C.; et al. Progress and issues on key
technologies in forecasting of snowmelt flood disaster in arid areas Northwest China. Adv. Earth Sci. 2021, 36, 233–244.
34. Li, H.; Zhao, Q.; Wu, J.; Ding, Y. Quantitative simulation of the runoff components and its variation characteristics in the upstream
of the Shule River. J. Glaciol. Geocryol. 2019, 41, 907–917.
35. Li, H. Study on the Change of Runoff and Its Components in an Alpine Watershed Based on SPHY Model; University of Chinese Academy
of Sciences: Beijing, China, 2019.
36. Zhou, J. Multi-Methods Study on Runoff Components in Alpine Regions: A Case Study of the Upper Reaches of Shule River Basin;
University of Chinese Academy of Sciences: Beijing, China, 2019.
37. Zhou, J.; Ding, Y.; Wu, J.; Liu, F.; Wang, S. Streamflow generation in semi-arid, glacier-covered, montane catchments in the upper
shule river, Qilian Mountains, Northeastern Tibetan Plateau. Hydrol. Process. 2021, 35, e14276.
38. Wu, J.-K.; Wu, X.-P.; Hou, D.-J.; Liu, S.-W.; Zhang, X.-Y.; Qin, X. Streamwater hydrograph separation in an alpine glacier area in
the Qilian Mountains, northwestern China. Int. Assoc. Sci. Hydrol. Bull. 2016, 61, 2399–2410. [CrossRef]
39. Zhou, J.; Wu, J.; Liu, S.; Zeng, G.; Qin, J.; Wang, X.; Zhao, Q. Hydrograph separation in the headwaters of the Shule river basin:
Combining water chemistry and stable isotopes. Adv. Meteorol. 2015, 2015, 830306.
40. Zhou, J.; Ding, Y.; Wu, J.; Wang, J. Applied comparison analysis of base-flow separation methods in upper Shule River basin. J.
Glaciol. Geocryol. 2019, 41, 1456–1466.
41. Zhang, Z.; Deng, S.; Zhao, Q.; Zhang, S.; Zhang, X. Projected glacier meltwater and river run-off changes in the Upper Reach of
the Shule River Basin, north-eastern edge of the Tibetan Plateau. Hydrol. Process. 2019, 33, 1059–1074. [CrossRef]
42. Zhao, Q.; Ding, Y.; Wang, J.; Gao, H.; Zhang, S.; Zhao, C.; Xu, J.; Han, H.; Shangguan, D. Projecting climate change impacts on
hydrological processes on the Tibetan Plateau with model calibration against the glacier inventory data and observed streamflow.
J. Hydrol. 2019, 573, 60–81. [CrossRef]
43. Jin, Z.; Zhao, Q.; Qin, X.; Zhang, J.; Zhang, H.; Qin, J.; Qin, Y.; Li, H.; Chen, J.; Liu, Y.; et al. Quantifying the impact of landscape
changes on hydrological variables in the alpine and cold region using hydrological model and remote sensing data. Hydrol.
Process. 2021, 35, e14392. [CrossRef]
44. He, X. Runoff Simulation and Response to Climate Change in the Upper Reaches of the HeiHe River Basin Based on SWAT Model; Northwest
Normal University: Lanzhou, China, 2013.
Water 2021, 13, 3357 20 of 21

45. Zhang, G.; Xie, H.; Yao, T.; Li, H.; Duan, S. Quantitative water resources assessment of Qinghai Lake basin using Snowmelt
Runoff Model (SRM). J. Hydrol. 2014, 519, 976–987. [CrossRef]
46. 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. [CrossRef]
47. Terink, W.; Lutz, A.F.; Simons, G.W.H.; Immerzeel, W.W.; Droogers, P. SPHY v2.0: Spatial processes in hydrology. Geosci. Model
Dev. 2015, 8, 2009–2034. [CrossRef]
48. Gao, M. The Analyzing and Simulating of Hydrologic Processes in the Upstream of the Shule River Basin; University of Chinese Academy
of Sciences: Beijing, China, 2012.
49. Liu, S.; Guo, W.; Xu, J. The Second Glacier Inventory Dataset of China (Version 1.0). National Cryosphere Desert Data Center.
2012. Available online: https://data.tpdc.ac.cn/en/data/f92a4346-a33f-497d-9470-2b357ccb4246/ (accessed on 20 November
2021).
50. Liu, S.; Yao, X.; Wan, W.; Xu, J.; Shangguan, D.; Wei, J.; Bao, W.; Wu, L. The contemporary glaciers in China based on the Second
Chinese Glacier Inventory. Acta Geogr. Sin. 2015, 70, 3–16.
51. 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. 2017, 61, 357–372. [CrossRef]
52. Ling, L.; Chao, L. The Chinese Soil Dataset Based on Harmonized World Soil Database Version (HWSD) (v1.1). National
Cryosphere Desert Data Center. 2019. Available online: https://data.tpdc.ac.cn/zh-hans/data/611f7d50-b419-4d14-b4dd-4a944
b141175/ (accessed on 20 November 2021).
53. Radoux, J.; Bontemps, S.; Defourny, P.; Eric, V.; Lamarche, C.; Frédéric, A.; Mayaux, P.; Martin, B.; Brockmann, C.; Kirches, G.;
et al. Consistent Global Land Cover Maps for Climate Modeling Communities: Current Achievements of the ESA’s Land Cover
CCI. In Proceedings of the ESA Living Planet Symposium, Edinburgh, UK, 9–13 September 2013.
54. Ouellet Dallaire, C.; Lehner, B.; Sayre, R.; Thieme, M. A multidisciplinary framework to derive global river reach classifications at
high spatial resolution. Environ. Res. Lett. 2018, 14, 024003. [CrossRef]
55. Li, J.; Wang, Y.; Li, J. Glacier Outlines over the Qilian Mountain Area (1980-2015). National Tibetan Plateau Data Center. 2019.
Available online: https://data.tpdc.ac.cn/en/data/890e5ad3-56e4-427b-b4b4-f1b4cfc21b2c/ (accessed on 20 November 2021).
56. Frey, H.; Machguth, H.; Huss, M.; Huggel, C.; Bajracharya, S.; Bolch, T.; Kulkarni, A.; Linsbauer, A.; Salzmann, N.; Stoffel, M.
Estimating the volume of glaciers in the Himalayan–Karakoram region using different methods. Cryosphere 2014, 8, 2313–2333.
[CrossRef]
57. Zou, X.; Gao, H.; Zhang, Y.; Ma, N.; Wu, J.; Farhan, S.B. Quantifying ice storage in upper Indus river basin using ground-
penetrating radar measurements and glacier bed topography model version 2. Hydrol. Process. 2021, 35, e14145. [CrossRef]
58. Linsbauer, A.; Frey, H.; Haeberli, W.; Machguth, H.; Azam, M.F.; Allen, S. Modelling glacier-bed overdeepenings and possible
future lakes for the glaciers in the Himalaya—Karakoram region. Ann. Glaciol. 2016, 57, 119–130. [CrossRef]
59. Haeberli, W.; Hoelzle, M. Application of inventory data for estimating characteristics of and regional climate-change effects on
mountain glaciers: A pilot study with the European Alps. Ann. Glaciol. 1995, 21, 206–212. [CrossRef]
60. Paterson, W. The Physics of Glaciers, 3rd ed.; Pergamon Press: Oxford, UK, 1994.
61. Karssenberg, D. The value of environmental modeling languages for building distributed hydrological models. Hydrol. Process.
2002, 16, 2751–2766. [CrossRef]
62. Karssenberg, D.; Schmitz, O.; Salamon, P.; de Jong, K.; Bierkens, M.F.P. A software framework for construction of process-based
stochastic spatio-temporal models nd data assimilation. Environ. Model. Softw. 2010, 25, 489–502. [CrossRef]
63. Schmitz, O.; Karssenberg, D.; de Jong, K.; Kok, J.-L.; de Jong, S. Map algebra and model algebra for integrated model building.
Environ. Model. Softw. 2013, 48, 113–128. [CrossRef]
64. Liang, X.; Wood, E.F.; Lettenmaier, D.P. Surface soil moisture parameterization of the VIC-2L model: Evaluation and modification.
Glob. Planet. Chang. 1996, 13, 195–206. [CrossRef]
65. Liang, X.; Lettenmaier, D.P.; Wood, E.; Burges, S.J. A simple hydrologically based model of land surface water and energy fluxes
for general circulation models. J. Geophys. Res. Atmos. 1994, 99, 14415–14428. [CrossRef]
66. Terink, W.; Lutz, A.F.; Immerzeel, W.; Nepal, S.; Khanal, S.; Shrestha, A.B. Improvement of the SPHY model glacier module and
its application in the Tamakoshi River Basin, Nepal. In Proceedings of the Agu Fall Meeting, San Francisco, CA, USA, 12–16
December 2016.
67. Santhi, C.; Arnold, J.; Williams, J.R.; Dugas, W.A.; Srinivasan, R.; Hauck, L.M. Validation of the SWAT model on a large river
basin with point and nonpoint sources. JAWRA J. Am. Water Resour. Assoc. 2001, 37, 1169–1188. [CrossRef]
68. Chen, R.; Kang, E.; Ding, Y. Some knowledge on and parameters of China’s alpine hydrology. Adv. Water Sci. 2014, 25, 307–317.
69. Zhang, Y.; Liu, S.; Ding, Y. Spatial variation of degree-day factors on the observed glaciers in Western China. Acta Geogr. Sin. 2006,
61, 89–98. [CrossRef]
70. Seibert, J. Estimation of parameter uncertainty in the HBV model. Hydrol. Res. 1997, 28, 247–262. [CrossRef]
71. Chen, R.; Liu, J.; Kang, E.; Yang, Y.; Han, C.; Liu, Z.; Song, Y.; Qing, W.; Zhu, P. Precipitation measurement intercomparison in the
Qilian Mountains, north-eastern Tibetan Plateau. Cryosphere 2015, 9, 1995–2008. [CrossRef]
72. Wang, L.; Chen, R.; Song, Y.; Yang, Y.; Liu, J.; Han, C.; Liu, Z. Precipitation–altitude relationships on different timescales and at
different precipitation magnitudes in the Qilian Mountains. Theor. Appl. Climatol. 2017, 134, 875–884. [CrossRef]
Water 2021, 13, 3357 21 of 21

73. Wang, L. Observation and Experimental Study on the Spatial-Temporal Variation of Precipitation in the Qilian Mountains; University of
Chinese Academy of Sciences: Beijing, China, 2020.
74. Zhang, Y.; Liu, S.; Ding, Y. Observed degree-day factors and their spatial variation on glaciers in western China. Ann. Glaciol.
2017, 43, 301–306. [CrossRef]