Change in frozen soils and its effect on regional hydrology, upper Heihe basin, northeastern Qinghai–Tibetan Plateau

Abstract

Frozen ground has an important role in regional hydrological cycles and ecosystems, particularly on the Qinghai–Tibetan Plateau (QTP), which is characterized by high elevations and a dry climate. This study modified a distributed, physically based hydrological model and applied it to simulate long-term (1971–2013) changes in frozen ground its the effects on hydrology in the upper Heihe basin, northeastern QTP. The model was validated against data obtained from multiple ground-based observations. Based on model simulations, we analyzed spatio-temporal changes in frozen soils and their effects on hydrology. Our results show that the area with permafrost shrank by 8.8 % (approximately 500 km${}^{\mathrm{2}})$, predominantly in areas with elevations between 3500 and 3900 m. The maximum depth of seasonally frozen ground decreased at a rate of approximately 0.032 m decade${}^{-\mathrm{1}}$, and the active layer thickness over the permafrost increased by approximately 0.043 m decade${}^{-\mathrm{1}}$. Runoff increased significantly during the cold season (November–March) due to an increase in liquid soil moisture caused by rising soil temperatures. Areas in which permafrost changed into seasonally frozen ground at high elevations showed especially large increases in runoff. Annual runoff increased due to increased precipitation, the base flow increased due to changes in frozen soils, and the actual evapotranspiration increased significantly due to increased precipitation and soil warming. The groundwater storage showed an increasing trend, indicating that a reduction in permafrost extent enhanced the groundwater recharge.

Introduction

Global warming has led to significant changes in frozen soils, including both permafrost and seasonally frozen ground at high latitudes and high elevations (Hinzman et al., 2013; Cheng and Wu, 2007). Changes in frozen soils can greatly affect land–atmosphere interactions and the energy and water balances of the land surface (Subin et al., 2013; Schuur et al., 2015), altering soil moisture, water flow pathways, and stream flow regimes (Walvoord and Kurylyk, 2016). Understanding the changes in frozen soils and their impacts on regional hydrology is important for water resources management and ecosystem protection in cold regions.

Previous studies based on either experimental observations or long-term meteorological or hydrological observations have examined changes in frozen soils and their impacts on hydrology. Several studies reported that permafrost thawing might enhance base flow in the arctic and the subarctic (Walvoord and Striegl, 2007; St. Jacques and Sauchyn, 2009; Ye et al., 2009), as well as in northeastern China (Liu et al., 2003; Duan et al., 2017). A few studies reported that permafrost thawing might reduce river runoff (here, runoff is defined as all liquid water flowing out of the study area), especially on the Qinghai–Tibetan Plateau (QTP) (e.g., Qiu, 2012; Jin et al., 2009). Those studies that include intensive field observations of frozen soils have typically been performed at small spatial scales over short periods (e.g., Cheng and Wu, 2007; Wu et al., 2010). Consequently, regional patterns and long-term trends are not typically captured. Long-term meteorological and hydrological observations are available in many areas, but they do not provide information on soil freezing and thawing processes (McClelland et al., 2004; Liu et al., 2003; Niu et al., 2011). Therefore, previous observation-based studies have not provided a sufficient understanding of the long-term changes in frozen soils and their impact on regional hydrology (Woo et al., 2008).

As an alternative strategy, hydrological models have been coupled with soil freezing–thawing schemes to simulate impacts of the changes in frozen soils on catchment hydrology. Several hydrological models (Rawlins et al., 2003; Chen et al., 2008) used simple freezing–thawing schemes, but these cannot simulate the vertical soil temperature profiles. The modified Variable Infiltration Capacity (VIC) model (Cherkauer and Lettenmaier, 1999) and the Community Land Model (CLM) (Oleson et al., 2010) simulate vertical soil freezing–thawing processes, but they simplify the flow routing using linear schemes. Subin et al. (2013) and Lawrence et al. (2015) used the CLM to simulate global changes in permafrost. Cuo et al. (2015) used the VIC to simulate frozen soil changes and their hydrological impacts on the plot scale in the headwaters of the Yellow River. The GEOtop model (Endrizzi et al., 2014) simulates three-dimensional water flux and vertical heat transfer in soil, but it is difficult to apply to regional investigations. Wang et al. (2010) and Y. L. Zhang et al. (2013) incorporated frozen soil schemes in a distributed hydrological model and showed improved performance in a small mountainous catchment. More regional studies are necessary to better understand the frozen soil changes and their impacts on regional hydrologic processes and water resources.

The QTP is known as Asia's water tower, and runoff changes on the plateau have significant impacts on water security in downstream regions (Immerzeel et al., 2010); hence, such changes have attracted considerable attention in recent years (Cuo et al., 2014). The QTP is characterized by high elevations and a cold climate. Consequently, cryospheric processes have great impacts on its hydrological processes (Cheng and Jin, 2013; Cuo et al., 2014). The thickness of permafrost on the QTP varies from 1 to 130 m, and the temperature ranges between $-\mathrm{0.5}$ and $-\mathrm{3.5}$ ${}^{\circ }$C (Yang et al., 2010). Compared with arctic and subarctic soils, the frozen soils on the QTP are more sensitive to increased air temperature (Yang et al., 2010), and changes in the frozen soils may have more significant impacts on the regional hydrology.

Clear increases in annual and seasonal air temperatures have been observed on the QTP (Li et al., 2005; Liu and Chen, 2000; Zhao et al., 2004). Several studies have shown changes in frozen soils based on long-term observations. For example, Cheng and Wu (2007) analyzed soil temperature profiles from boreholes on the QTP and found that the active layer thickness of frozen soils increased by 0.15–0.50 m during the period of 1996–2001. Zhao et al. (2004) observed a decreasing trend of freezing depth in the seasonally frozen soils at 50 stations. Several studies have analyzed the relationship between changes in frozen soils and river discharge using observational data (Zhang et al., 2003; Jin et al., 2009; Niu et al., 2011). However, the spatio-temporal characteristics of the long-term changes in frozen soils are not sufficiently clear. Based on comprehensive field experiments (Cheng et al., 2014), a hydrological model coupling cryospheric processes and hydrological processes has been developed (Yang et al., 2015; Gao et al., 2016). This model provides a basis upon which to analyze the spatio-temporal changes in frozen soils and their impacts on the regional hydrology in the upper Heihe basin, northeastern QTP. Specifically, this study aims to (1) explore the spatial and temporal changes in frozen soils using a distributed hydrological model with comprehensive validation and (2) analyze the hydrological responses to the changes in frozen soils during the past 40 years in the upper Heihe basin.

The study area, hydrological stations, borehole observation and flux tower stations.

[Figure omitted. See PDF]

Study area and data

The Heihe River is one of the major inland basins in northwestern China. As shown in Fig. 1, the upper reaches of the Heihe River, representing a drainage area of 10 009 km${}^{\mathrm{2}}$, are located on the northeastern QTP at elevations ranging from 2200 to 5000 m. The upper reaches of this river provide the majority of the water supplied to the middle and lower reaches (Cheng et al., 2014). The annual precipitation in the upper Heihe basin ranges from 200 to 700 mm, and the mean annual air temperature ranges from $-\mathrm{9}$ to 5 ${}^{\circ }$C. Permafrost dominates the high-elevation region above 3700 m (Wang et al., 2013), and seasonal frozen ground covers the remaining portion of the study area. Glaciers are found at elevations above 4000 m and cover approximately 0.8 % of the upper Heihe basin. The upper Heihe basin contains two tributaries, each with a hydrological station; one at Qilian (on the eastern tributary) and the other at Zhamashike (on the western tributary). The outlet of the upper Heihe basin also features a hydrological station, namely Yingluoxia (Fig. 1).

The spatial data used in this study includes: atmospheric forcing data, land surface data, and actual evapotranspiration data based on remote sensing. The atmospheric forcing data include a 1 km gridded dataset of daily precipitation, air temperature, sunshine hours, wind speed, and relative humidity. The gridded daily precipitation was interpolated from observations at meteorological stations (Fig. 1) provided by the China Meteorological Administration (CMA) using the method developed by Wang et al. (2017). The other atmospheric forcing data were interpolated from observations at meteorological stations using the inverse distance weighted method. The interpolation of air temperature considers the elevation-dependent temperature gradient provided by the Heihe Watershed Allied Telemetry Experimental Research (HiWATER) experiment (Li et al., 2013).

The land surface data used to run the model include land use, topography, leaf area index (LAI), and soil parameters. The topography data were obtained from the Shuttle Radar Topography Mission (SRTM) dataset (Jarvis et al., 2008) with a spatial resolution of 90 m. The land use and land cover data were provided by the Institute of Botany, Chinese Academy of Sciences (Zhou and Zheng, 2014). The LAI data with 1 km resolution were developed by Fan (2014). The soil parameters were developed by Song et al. (2016) and include the saturated hydraulic conductivity, residual soil moisture content, saturated soil moisture content, soil sand content, soil clay content, and soil organic content. Monthly actual evapotranspiration data with 1 km resolution during the period of 2002–2012 were estimated based on remote sensing data (Wu et al., 2012; Wu, 2013).

The field observation data used in this study include river discharge, soil temperature, frozen depth, soil moisture, and borehole observations. Daily river discharge data were obtained from the Hydrology and Water Resources Bureau of Gansu Province. The CMA provided daily soil temperature data collected at the Qilian station from 1 January 2004 to 31 December 2013, and daily frozen depth data collected at the Qilian and Yeniugou stations from 1 January 2002 to 31 December 2013. We obtained ground temperature observations from six boreholes, the locations of which are shown in Fig. 1, from Wang et al. (2013). We used the observations at specific dates instead of annual averages due to a lack of continuous measurements. The borehole depths are 100 m for T1, 69 m for T2, 50 m for T3, 90 m for T4, and 20 m for T5 and T7. The HiWATER experiment (Li et al., 2013; Liu et al., 2011) provided the soil moisture data from 1 January to 31 December 2014 at the A'rou Sunny Slope station (100.52${}^{\circ }$ E, 38.09${}^{\circ }$ N).

Model structure and vertical discretization of soil column.

[Figure omitted. See PDF]

Methodology

Brief introduction of the hydrological model

This study used the distributed Geomorphology-Based Ecohydrological Model (GBEHM), which was developed by Yang et al. (2015) and Gao et al. (2016). The GBEHM is a spatially distributed model for large-scale river basins. It employs geomorphologic properties to reduce the lateral two dimensions into one dimension for flow routing within a sub-catchment, which greatly improves computational efficiency while retaining the spatial heterogeneity in water flow paths at the basin scale. As shown in Fig. 2, the GBEHM used a 1 km grid system to discretize the study catchment, which was divided into 251 sub-catchments. Each sub-catchment was further divided into flow intervals along its main stream. To capture the sub-grid topography, each 1 km grid was represented by a number of hillslopes with an average length and gradient, but different aspect, which were estimated from the 90 m digital elevation models (DEM). Additional hillslope properties include soil and vegetation types (Yang et al., 2015).

The hillslope is the basic unit in the hydrological simulation of the water and heat transfers (both conduction and convection) in the vegetation canopy, snow/glacier, and soil layers. The canopy interception, radiation transfer in the canopy and the energy balance of the land surface are described using the methods of Simple Biosphere Model 2 (SIB2) (Sellers et al., 1985, 1996). The surface runoff on the hillslope is solved using the kinematic wave equation. The groundwater aquifer is considered as individual storage unit corresponding to each grid. Exchange between the groundwater and the river water is calculated using Darcy's law (Yang et al., 1998, 2002; Cong et al., 2009).

The model runs with a time step of 1 h. Runoff generated from the grid is the lateral inflow into the river over the same flow interval in the corresponding sub-catchment. Flow routing in the river network is calculated using the kinematic wave equation following the sequence determined by the Horton–Strahler scheme (Strahler, 1957). The model is driven by the atmosphere forcing data and land surface data introduced in Sect. 2.

Simulation of cryospheric processes

The simulation of cryospheric processes in the GBEHM includes glacier ablation, snow melting, and soil freezing and thawing.

Glacier ablation

Glacier ablation is simulated using the following energy balance model (Oerlemans, 2001): $$Q_{\mathrm{M}}=\mathrm{SW}(\mathrm{1}-\mathit{\alpha })+{\mathrm{LW}}_{\mathrm{in}}-{\mathrm{LW}}_{\mathrm{out}}-Q_{\mathrm{H}}-Q_{\mathrm{L}}-Q_{\mathrm{G}}+Q_{\mathrm{R}},$$ where $Q_{\mathrm{M}}$ is the net energy absorbed by the surface of the glacier (W m${}^{-\mathrm{2}})$; SW is the incoming shortwave radiation (W m${}^{-\mathrm{2}})$; $\mathit{\alpha }$ is the surface albedo; LW${}_{\mathrm{in}}$ is the incoming longwave radiation (W m${}^{-\mathrm{2}})$; LW${}_{\mathrm{out}}$ is the outgoing longwave radiation (W m${}^{-\mathrm{2}})$; $Q_{\mathrm{H}}$ is the sensible heat flux (W m${}^{-\mathrm{2}})$; $Q_{\mathrm{L}}$ is the latent heat flux (W m${}^{-\mathrm{2}})$; $Q_{\mathrm{R}}$ is the energy from rainfall (W m${}^{-\mathrm{2}})$; and $Q_{\mathrm{G}}$ is the penetrating shortwave radiation (W m${}^{-\mathrm{2}})$. The surface albedo is calculated as follows (Oerlemans and Knap, 1998): $$\mathit{\alpha }={\mathit{\alpha }}_{\mathrm{snow}}+({\mathit{\alpha }}_{\mathrm{ice}}-{\mathit{\alpha }}_{\mathrm{snow}})e^{-h/d^{*}},$$ where ${\mathit{\alpha }}_{\mathrm{snow}}$ is the albedo of snow on the glacier surface; ${\mathit{\alpha }}_{\mathrm{ice}}$ is the albedo of the ice surface; $h$ is the snow depth on the glacier surface (m); $d^{*}$ is a parameter describing the snow depth effect on the albedo (m).

The amount of melt water is calculated as follows (Oerlemans, 2001): $$M={\displaystyle \frac{Q_{\mathrm{M}}}{L_{\mathrm{f}}}}\mathrm{d}t,$$ where d$t$ is the time step used in the model (s) and $L_{\mathrm{f}}$ is the latent heat of fusion (J kg${}^{-\mathrm{1}})$.

Snowmelt

A multilayer snow cover model is used to describe the mass and energy balance of snow cover. The snow parametrization is based on Jordan (1991), and two constituents, namely, ice and liquid water, are used to describe each snow layer. For each snow layer, temperature is solved using an energy balance approach (Bartelt and Lehnin, 2002): $$C_{\mathrm{s}}{\displaystyle \frac{\partial T_{\mathrm{s}}}{\partial t}}-L_{\mathrm{f}}{\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{i}}{\mathit{\theta }}_{\mathrm{i}}}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left(K_{\mathrm{s}},{\displaystyle \frac{\partial T}{\partial z}}\right)+{\displaystyle \frac{\partial I_{\mathrm{R}}}{\partial z}}+Q_{\mathrm{R}},$$ where $C_{\mathrm{s}}$ is the heat capacity of snow (J m${}^{-\mathrm{3}}$ K${}^{-\mathrm{1}})$; $T_{\mathrm{s}}$ is the temperature of the snow layer (K); ${\mathit{\rho }}_{\mathrm{i}}$ is the density of ice (kg m${}^{-\mathrm{3}})$; ${\mathit{\theta }}_{\mathrm{i}}$ is the volumetric ice content; $K_{\mathrm{s}}$ is the thermal conductivity of snow (W m${}^{-\mathrm{1}}$ K${}^{-\mathrm{1}})$; $L_{\mathrm{f}}$ is the latent heat of ice fusion (J kg${}^{-\mathrm{1}})$ ; $I_{\mathrm{R}}$ is the radiation transferred into the snow layer (W m${}^{-\mathrm{2}})$; and $Q_{\mathrm{R}}$ is the energy delivered by rainfall (W m${}^{-\mathrm{2}})$, which is only considered for the top snow layer. The solar radiation transfer in the snow layers and the snow albedo are simulated using the SNow, ICe, and Aerosol Radiative (SNICAR) model, which is solved using the method developed by Toon et al. (1989). Equation (4) is solved using an implicit centered finite difference method, and a Crank–Nicholson scheme is employed.

The mass balance of the snow layer is described as follows (Bartelt and Lehnin, 2002): $$\begin{array}{}& {\displaystyle }& {\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{i}}{\mathit{\theta }}_{\mathrm{i}}}{\partial t}}+M_{\mathrm{iv}}+M_{\mathrm{il}}=\mathrm{0}& {\displaystyle }& {\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{l}}{\mathit{\theta }}_{\mathrm{l}}}{\partial t}}+{\displaystyle \frac{\partial U_{\mathrm{l}}}{\partial z}}+M_{\mathrm{lv}}-M_{\mathrm{il}}=\mathrm{0},\end{array}$$ where ${\mathit{\rho }}_{\mathrm{l}}$ is the density of the liquid water (kg m${}^{-\mathrm{3}})$; ${\mathit{\theta }}_{\mathrm{l}}$ is the volumetric liquid water content; $U_{\mathrm{l}}$ is the liquid water flux (kg m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}})$; $M_{\mathrm{iv}}$ is the mass of ice that changes into vapor within a time step (kg m${}^{-\mathrm{3}}$ s${}^{-\mathrm{1}})$; $M_{\mathrm{il}}$ is the mass of ice that changes into liquid water within a time step (kg m${}^{-\mathrm{3}}$ s${}^{-\mathrm{1}})$; and $M_{\mathrm{lv}}$ is the mass of liquid water that changes into vapor within a time step (kg m${}^{-\mathrm{3}}$ s${}^{-\mathrm{1}})$. The liquid water flux of the snow layer is calculated as follows (Jordan, 1991): $$U_{\mathrm{l}}=-{\displaystyle \frac{k}{{\mathit{\mu }}_{\mathrm{l}}}}{\mathit{\rho }}_{\mathrm{l}}^{\mathrm{2}}g,$$ where $k$ is the hydraulic permeability (m${}^{\mathrm{2}})$, ${\mathit{\mu }}_{\mathrm{l}}$ is dynamic viscosity of water at 0 ${}^{\circ }$C ($\mathrm{1.787}\times {\mathrm{10}}^{-\mathrm{3}}$ N s m${}^{-\mathrm{2}})$, ${\mathit{\rho }}_{\mathrm{l}}$ is the density of liquid water (kg m${}^{-\mathrm{3}})$ and $g$ is gravitational acceleration (m s${}^{-\mathrm{2}})$. The water flux of the bottom snow layer is considered snowmelt runoff.

Soil freezing and thawing

The energy balance of the soil layer is solved as follows (Flerchinger and Saxton, 1989): $$C_{\mathrm{s}}{\displaystyle \frac{\partial T}{\partial t}}-{\mathit{\rho }}_{\mathrm{i}}{L}_{\mathrm{f}}{\displaystyle \frac{\partial {\mathit{\theta }}_{\mathrm{i}}}{\partial t}}-{\displaystyle \frac{\partial }{\partial z}}\left({\mathit{\lambda }}_{\mathrm{s}},{\displaystyle \frac{\partial T}{\partial z}}\right)+{\mathit{\rho }}_{\mathrm{l}}{c}_{\mathrm{l}}{\displaystyle \frac{\partial q_{\mathrm{l}}T}{\partial z}}=\mathrm{0},$$ where $C_{\mathrm{s}}$ is the volumetric soil heat capacity (J m${}^{-\mathrm{3}}$ K${}^{-\mathrm{1}})$; $T$ is the temperature (K) of the soil layers; $z$ is the vertical depth of the soil (m); ${\mathit{\theta }}_{\mathrm{i}}$ is the volumetric ice content; ${\mathit{\rho }}_{\mathrm{i}}$ is the density of ice (kg m${}^{-\mathrm{3}})$; ${\mathit{\lambda }}_{\mathrm{s}}$ is the thermal conductivity (W m${}^{-\mathrm{1}}$ K${}^{-\mathrm{1}})$; ${\mathit{\rho }}_{\mathrm{l}}$ is the density of liquid water (kg m${}^{-\mathrm{3}})$; and $c_{\mathrm{l}}$ is the specific heat of liquid water (J kg${}^{-\mathrm{1}}$ K${}^{-\mathrm{1}})$. In addition, $q_{\mathrm{l}}$ is the water flux between different soil layers (m s${}^{-\mathrm{1}})$ and is solved using the 1-D vertical Richards equation. The unsaturated soil hydraulic conductivity is calculated using the modified van Genuchten's equation (Wang et al., 2010), as follows: $$\begin{array}{cc}& {\displaystyle }& {\displaystyle }K={\displaystyle }& {\displaystyle }{f}_{\mathrm{ice}}{K}_{\mathrm{sat}}{\left({\displaystyle \frac{{\mathit{\theta }}_{\mathrm{l}}-{\mathit{\theta }}_{\mathrm{r}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}}\right)}^{\mathrm{1}/\mathrm{2}}{\left[\mathrm{1},-,{\left(\mathrm{1},-,{\left({\displaystyle \frac{{\mathit{\theta }}_{\mathrm{l}}-{\mathit{\theta }}_{\mathrm{r}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}}\right)}^{-\mathrm{1}/m}\right)}^m\right]}^{\mathrm{2}},\end{array}$$ where $K$ is the unsaturated soil hydraulic conductivity (m s${}^{-\mathrm{1}})$; $K_{\mathrm{sat}}$ is the saturated soil hydraulic conductivity (m s${}^{-\mathrm{1}})$; ${\mathit{\theta }}_{\mathrm{l}}$ is the volumetric liquid water content; ${\mathit{\theta }}_{\mathrm{s}}$ is the saturated water content; ${\mathit{\theta }}_{\mathrm{r}}$ is the residual water content; $m$ is an empirical parameter in van Genuchten's equation and $f_{\mathrm{ice}}$ is an empirical hydraulic conductivity reduction factor that is calculated using soil temperature as follows (Wang et al., 2010): $$f_{\mathrm{ice}}=\exp [-\mathrm{10}(T_{\mathrm{f}}-T_{\mathrm{soil}}\left)\right],\mathrm{0.05}\le f_{\mathrm{ice}}\le \mathrm{1},$$ where $T_{\mathrm{f}}$ is 273.15 K and $T_{\mathrm{soil}}$ is the soil temperature.

Equation (8) solves the soil temperature with the upper boundary condition as the heat flux into the uppermost soil layer. When the ground is not covered by snow, the heat flux from the atmosphere into the uppermost soil layer is expressed as follows (Oleson et al., 2010): $$h=S_{\mathrm{g}}+L_{\mathrm{g}}-H{}_{\mathrm{g}}-\mathit{\lambda }{E}_{\mathrm{g}}+Q_{\mathrm{R}},$$ where $h$ is the upper boundary heat flux into the soil layer (W m${}^{-\mathrm{2}})$; $S_{\mathrm{g}}$ is the solar radiation absorbed by the uppermost soil layer (W m${}^{-\mathrm{2}})$; $L_{\mathrm{g}}$ is the net long-wave radiation absorbed by the ground (W m${}^{-\mathrm{2}})$, $H_{\mathrm{g}}$ is the sensible heat flux from the ground (W m${}^{-\mathrm{2}})$; $\mathit{\lambda }{E}_{\mathrm{g}}$ is the latent heat flux from the ground (W m${}^{-\mathrm{2}})$; and $Q_{\mathrm{R}}$ is the energy delivered by rainfall (W m${}^{-\mathrm{2}})$. When the ground is covered by snow, the heat flux into the uppermost soil layer is calculated as follows: $$h=I_{\mathrm{p}}+G,$$ where $I_{\mathrm{p}}$ is the radiation that penetrates the snow cover, and $G$ is the heat conduction from the bottom snow layer to the uppermost soil layer. Equation (8) is solved using a finite difference scheme with an hourly time step, similar to the solution of Eq. (4).

There are no available observations of the geothermal heat flux for the northeastern QTP. To simulate permafrost we assume an upward thermal heat flux at the bottom boundary and estimate its value to be 0.14 W m${}^{-\mathrm{2}}$ at a depth of 50 m using the average geothermal gradient from the four boreholes (T1–T4) shown in Fig. 3, which is reasonable based on a comparison with the observations (0.02–0.16 W m${}^{-\mathrm{2}})$ from the interior of the QTP (Wu et al., 2010). The vertical soil column is divided into 39 layers in the model (Fig. 2). The 1.7 m topsoil layer is subdivided into 9 layers. The first soil layer at the surface is 0.05 m, and the layer thicknesses increase with depth linearly from 0.05 to 0.3 m at a depth of 0.8 m. Then thicknesses decrease linearly with depth reaching 0.1 m at a depth of 1.7 m. From 1.7 to 3.0 m, there are 12 soil layers with a constant thickness of 0.1 m to try to capture the maximum freezing depths according to field observations. From the depth of 3 to 50 m, there are 18 layers with thicknesses increasing exponentially from 0.1 to 12 m. The liquid soil moisture, ice content, and soil temperature of each layer is calculated at each time step. The soil heat capacity and soil thermal conductivity are estimated using the method developed by Farouki (1981). Permafrost is defined as ground with a temperature at or below 0${}^{\circ }$C for at least two consecutive years (Woo, 2012). This study differentiated permafrost from seasonally frozen ground based on the simulated vertical ground temperature profile in each grid. For each year in each grid, the frozen ground condition was determined by searching the ground temperature profile within a four-year window from the previous three years to the current year.

Comparison of the simulated and the observed soil temperature at borehole observation sites, with the observed data provided by Wang et al. (2013).

[Figure omitted. See PDF]

Model calibration

To initialize the model, we first estimated the soil temperature profiles based on the assumption that there is a linear relationship between the ground temperature at a given depth below the surface and elevation. This temperature–elevation relationship is estimated from the observed ground temperatures in six boreholes (see Fig. 1). Next, the model had a 500-year spin-up-run to specify the initial values of the hydrological variables (e.g., soil moisture, soil ice content, ground temperature, and groundwater table) by repeating the atmospheric forcing data from 1961 to 1970.

This study used the period of 2002 to 2006 for model calibration and the period of 2008 to 2012 for model validation. The daily ground temperature at the Qilian station and the frozen depths at the Qilian and Yeniugou stations were used to calibrate the ground surface reflectance according to vegetation type. The other parameters, such as groundwater hydraulic conductivity, were calibrated according to the observed baseflow discharge in the winter season at the Qilian, Zhamashike and Yingluoxia stations. The Nash–Sutcliffe efficiency (NSE) and relative error are calculated using observed and simulated discharge to evaluate the model performance. We calibrated the surface retention capacity and surface roughness to match the observed flood peaks, and calibrated the leaf reflectance, leaf transmittance, and maximum Rubsico capacity of the top leaf based on the remote sensing evapotranspiration data. Table 1 shows the major parameters used in the model.

We also simulated the hydrological processes without the frozen soil scheme in order to investigate the impact of frozen soils. In this case, the phase transition of soil water between the solid and the liquid is not considered, although the ground temperature is still simulated. Other processes are simulated in the same manner as in the normal run.

Major parameters of the GBEHM model.

Parameters	Coniferous	Shrub	Steppe	Alpine	Alpine Sparse	Desert
	forest			meadow	vegetation
Surface retention capacity (mm)	30.0	25.0	10.0	15.0	15.0	5.0
Surface roughness (Manning coefficient)	0.5	0.3	0.1	0.1	0.1	1.0
Soil reflectance to visible light	0.20	0.20	0.20	0.28	0.14	0.11
Soil reflectance to near-infrared radiation	0.225	0.225	0.225	0.28	0.225	0.225
Leaf reflectance to visible light	0.105	0.105	0.105	0.105	0.105	–
Leaf reflectance to near-infrared radiation	0.35	0.58	0.58	0.58	0.58	–
Leaf transmittance to visible light	0.05	0.07	0.07	0.07	0.07	–
Leaf transmittance to near-infrared radiation	0.10	0.25	0.25	0.25	0.25	–
Maximum Rubsico capacity of top leaf	6.0	6.0	3.3	3.3	3.0	–
(10${}^{-\mathrm{5}}$ mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}})$
Plant root depth (m)	2.0	1.0	0.40	0.40	0.1	0.0
Intrinsic quantum efficiency (mol mol${}^{-\mathrm{1}})$	0.08	0.08	0.05	0.05	0.05	–
Canopy top height (m)	9.0	1.9	0.3	0.3	0.2	–
Leaf length (m)	0.055	0.055	0.3	0.3	0.04	–
Leaf width (m)	0.001	0.001	0.005	0.005	0.001	–
Stem area index	0.08	0.08	0.05	0.05	0.08	–

Daily soil temperature at the Qilian station: (a) observation; (b) simulation; (c) difference (simulation $-$ observation).

[Figure omitted. See PDF]

Results

Validation of the hydrological model

We conducted a comprehensive validation of the GBEHM using the ground temperature profiles observed from six boreholes, the long-term observations of ground temperature and frozen depths from the Qilian and Zhamashike stations, the soil moisture observations from the A'rou Sunny Slope station, the long-term observations of streamflow from the three hydrological stations shown in Fig. 1, and the monthly actual evapotranspiration estimated from remote sensing data.

Figure 3 shows the comparison of the model-simulated and observed ground temperature profiles at the six boreholes. The model generally captured the vertical distribution of the ground temperature at T1, T2, T3, and T4 in the permafrost area, but the temperatures were overestimated above 20 m depth for T1 and T3. Good agreement between the simulated and observed ground temperature profiles below the depth of 20 m is probably due to the fitting of initial values. Therefore, the deep ground temperatures are stable, which is confirmed by the comparison of temperature profiles in different years, as shown in Fig. S1 in the Supplement. Figure S1 also illustrates that the temperatures above 20 m have shown significant increasing trends over the past 40 years. The errors in simulating the vertical temperature profile near the surface might be caused by simplification of the 3-D topography. At T5, which is located in seasonally frozen ground, the simulated ground temperature profile did not agree well with the observed profile at depths of 4–20 m. This error might also be related to heterogeneity in the ground properties, especially the thermal conductivity and heat capacity, as no such information was available. The model simulation agrees well with the borehole observations at T7, which is located in the transition zone from permafrost to seasonally frozen ground. Therefore, the model can identify the boundary between the permafrost and seasonally frozen ground.

We also validated the model simulation of the freezing–thawing cycles based on long-term observations of ground temperature and frozen depth. Figure 4 compares the simulated ground temperature with the observed temperature at the Qilian station, which is located in seasonally frozen ground (observed daily ground temperature data are available from 2004). Generally, the model simulations accurately captured the seasonal changes in the ground temperature profile. Validation of the ground temperature at different depths (0.05, 0.1, 0.2, 0.4, 0.8, 0.16, and 0.32 m) showed that the root mean square error (RMSE) decreases with increasing depth. The RMSE was approximately 2.5${}^{\circ }$ for the uppermost three depths (0.5, 0.10 and 0.2 m). The RMSE for depths of 0.4 cm and 0.8 m was 1.7 and 1.5${}^{\circ }$, respectively, and the RMSE for a depth of 3.2 m was 0.9${}^{\circ }$. Uncertainties in the simulations may be related to the ground heat capacity and thermal conductivity estimated according to Farouki (1981), and the results are similar to the findings by Ou et al. (2016) using the Northern Ecosystem Soil Temperature (NEST) model. We compared the model-simulated daily frozen depth with in situ observations at the Qilian and Yeniugou stations from 2002 to 2014, as shown in Fig. 5. The model reproduced the daily variations in frozen depth well, although the depth was underestimated by approximately 0.5 m at the Yeniugou station. In general, the validation of ground temperature and frozen depth indicates that the model effectively captured the freezing and thawing processes in the upper Heihe basin.

Comparison of the simulated and observed daily frozen depths during the period of 2002–2014 at (a) the Qilian station, and (b) the Yeniugou station.

[Figure omitted. See PDF]

Comparison of the simulated and the observed daily river discharge at (a) the Yingluoxia gauge, (b) the Qilian gauge, and (c) the Zhamashike gauge. For each gauge, the upper and lower panels show the calibration and validation periods, respectively. NSE and relative error coefficients are indicated.

[Figure omitted. See PDF]

Furthermore, we used the observed hourly liquid soil moisture at the A'rou Sunny Slope station for an additional independent validation. Figure S2 in the Supplement shows the comparison between the simulated and observed liquid soil moisture at different depths from 1 January to 31 December 2014. This comparison demonstrates that the model simulation of liquid soil moisture is reasonable.

Figure 6 compares the model-simulated and the observed daily streamflow discharge at the Yingluoxia, Qilian, and Zhamashike stations. The model simulations agreed well with the observations. The model simulations captured the flood peaks and the magnitude of base flow in both the calibration and validation periods. For the Yingluoxia, Qilian, and Zhamashike stations, the NSE coefficients were 0.64, 0.63, and 0.72, in the calibration period and 0.64, 0.60, and 0.73, in the validation period. The relative error (RE) was within 10 % for both the calibration and validation periods (Fig. 6). Figure S3 shows the comparison of the model-simulated monthly actual evaporation data and the remote-sensing-based evaporation data for the entire calibration and validation periods. The GBEHM simulation showed similar temporal variations in actual evapotranspiration compared with the remote sensing based estimation, and the RMSE of the simulated monthly evapotranspiration was 9.1 mm in the calibration period and 7.1 mm in the validation period.

We also compared the model-simulated river discharges with and without the frozen soil scheme. Table S1 in the Supplement shows that the model with the frozen soil scheme achieves a better simulation of the daily hydrograph than the model without the frozen soil scheme. Figure S4 shows that the model without the frozen soil scheme overestimates the river discharge in the freezing season and underestimates flood peaks in the warming season.

Simulated ground temperature changes in (a) the freezing season (from November to March) and (b) the thawing season (from April to October). Results from linear regressions are indicated.

[Figure omitted. See PDF]

Long-term changes in frozen soils

In the upper Heihe basin, the ground surface starts to freeze in November and begins to thaw in April (Wang et al., 2015a). From November to March, the ground surface temperature is below 0${}^{\circ }$ in both the permafrost and seasonally frozen ground regions, and precipitation mainly falls in the period from April to October. Therefore, to investigate the changes in frozen soils and their hydrological impact, a year is subdivided into two seasons, i.e., the freezing season (November to March) and the thawing season (April to October). Increasing precipitation and air temperature in the study area during both seasons over the past 50 years were reported in a previous study (Wang et al., 2015b). Compared to the decadal mean for 1971 to 1980, the annual mean air temperature for the 2001 to 2010 period was approximately 1.2${}^{\circ }$ higher, with a larger increase in the freezing season (1.4${}^{\circ }$) than in the thawing season (1.1${}^{\circ }$) (Table S2).

Figure 7 shows the changes in the basin-averaged ground temperature in the freezing and thawing seasons. The ground temperature increased in all seasons, especially over the past 30 years. The increasing trend of ground temperature was larger in the freezing season than in the thawing season. In the freezing season (Fig. 7a), the top-layer ground temperature was lower than the deep-layer temperature. The linear trend of the top layer (0–0.5 m) ground temperature was 0.49 ${}^{\circ }$C decade${}^{-\mathrm{1}}$ and the trend of the deep layer (2.5–3 m) temperature was 0.32 ${}^{\circ }$C decade${}^{-\mathrm{1}}$. The ground temperature in the deep layer (2.5–3 m) changed from $-\mathrm{0.7}$ ${}^{\circ }$C in the 1970s to approximately 0 ${}^{\circ }$C in the most recent decade. In the thawing season (Fig. 7b), the increasing trend of the top layer (0–0.5 m) ground temperature (0.29 ${}^{\circ }$C decade${}^{-\mathrm{1}})$ was greater than that of the deep layer (2.5–3 m) temperature (0.22 ${}^{\circ }$C decade${}^{-\mathrm{1}})$. The warming trend was larger in shallow ground layers; this is because the surface heat flux is impeded by the thermal inertia as it penetrates to greater depths.

Change of the frozen soils in the upper Heihe basin: (a) areas of permafrost and basin-average annual air temperature; (b) the basin-average annual maximum depths of seasonally frozen ground, and thaw above permafrost. Results from linear regressions are indicated.

[Figure omitted. See PDF]

Distribution of permafrost and seasonally frozen ground for two periods: (a) 1971–1980 and (b) 2001–2010; (c) area where permafrost degraded to seasonally frozen ground from (a) to (b); percentage of permafrost area with respect to elevation on the (d) sunny and (e) the shaded slopes for the two periods. Note that (d) and (e) share a legend.

[Figure omitted. See PDF]

Spatially averaged monthly ground temperatures simulated from 1971 to 2013 for two elevation intervals: (a) seasonally frozen ground between 3300 and 3500 m; (b) permafrost that degraded to seasonally frozen ground between 3500 and 3700 m.

[Figure omitted. See PDF]

Figure 8 shows the change in permafrost area during 1971–2013. As shown in Fig. 8a, the permafrost areas decreased by approximately 8.8 % (from 5700 km${}^{\mathrm{2}}$ in the 1970s to 5200 km${}^{\mathrm{2}}$ in the 2000s), indicating an evident decrease in the permafrost extent in the upper Heihe basin in the past 40 years.

Figure 8b shows the changes in the basin-averaged maximum frozen depth in the seasonally frozen ground areas and active layer thickness in the permafrost areas. The basin-averaged annual maximum frozen depth showed a significant decreasing trend (0.032 m decade${}^{-\mathrm{1}})$. In addition, the maximum frozen depth had a significantly negative correlation with the annual mean air temperature ($r=-\mathrm{0.71}$). Simulated active layer thickness in the permafrost regions increased (0.043 m decade${}^{-\mathrm{1}}$), and correlated positively with the annual mean air temperature ($p=\mathrm{0.005}$).

Changes in annual basin water balance and runoff components (mm year${}^{-\mathrm{1}})$ in different seasons.

							Runoff components
							Freezing season			Thawing season
					Runoff	Runoff	(from November			(from April
					ratio	ratio	to March)			to October)
Decade	$P$	$E$	Sim $R$	Obs $R$	(Obs)	(Sim)	$T$	$G$	$S$	$T$	$G$	$S$
1971–1980	439.1	282.1	154.1	143.8	0.33	0.35	18.5	0.0	0.0	135.6	3.5	13.8
1981–1990	492.8	300.8	188.5	174.1	0.35	0.38	20.5	0.0	0.0	168.0	3.1	27.8
1991–2000	471.0	307.6	161.9	157.4	0.33	0.34	20.5	0.0	0.0	141.4	3.8	18.4
2001–2010	504.3	319.0	180.6	174.3	0.35	0.36	26.2	0.0	0.0	154.3	3.7	24.1

Note: $P$ means precipitation, $E$ means actual evaporation, $R$ means runoff, $T$ means total runoff, and $G$ means glacier runoff. $S$ means snowmelt runoff, Sim means simulation, and Obs means observation.

Runoff and simulated evapotranspiration in (a) the freezing season and (b) the thawing season for the period of 1971 to 2013. Trend lines are for simulated data and regression results are shown. The upper two panels represent the freezing season and the lower two panels represent the thawing season.

[Figure omitted. See PDF]

Figure 9 shows the frozen soil distributions in the periods of 1971–1980 and 2001–2010. Comparing the frozen soil distributions of the two periods, we observed major changes in the frozen soils on sunny slopes at elevations between 3500 and 3900 m, especially in the west tributary, where large areas of permafrost changed into seasonally frozen ground. Figure S5, illustrating the taliks simulated in the period of 2001–2010, shows that taliks were mainly located on the edge of the permafrost area and that talik development was not significant.

Figure 10 shows the monthly mean ground temperatures for areas with elevations between 3300 and 3500 m and over areas with elevations between 3500 and 3700 m in the upper Heihe basin. In the areas with elevations between 3300 and 3500 m located in the seasonally frozen ground region (Fig. 10a), the frozen depth decreased and the ground temperature in deep layers (depths greater than 2 m) increased. Figure 10b shows that the increase in ground temperature was larger in the area with higher elevation (3500–3700 m). This figure shows that the thickness of the permafrost layer decreased as the ground temperature increased, and the permafrost changed into seasonally frozen ground after 2000. The surface thaw depths changed slowly compared with the depth to the base of permafrost as shown in Fig. 10, which may be primarily due to the geothermal heat flux. Additionally, the faster increase in the air temperature in the freezing season (0.41 ${}^{\circ }$C decade${}^{-\mathrm{1}}$) than in the thawing season (0.26 ${}^{\circ }$C decade${}^{-\mathrm{1}}$) may be another reason.

Changes in the water balance and runoff

Table 2 shows the decadal changes in the annual water balance from 1971 to 2010 based on the model simulation. The annual precipitation, annual runoff, and annual runoff ratio exhibited the same decadal variation. However, the annual evapotranspiration maintained an increasing trend starting in the 1970s that was consistent with the rising air temperature and soil warming. Although the actual evapotranspiration increased, the runoff ratio remained stable during the past 4 decades because of the increased precipitation.

Figure 11 and Table 2 show the changes in runoff (both simulated and observed) in different seasons. The model-simulated and observed runoff both exhibited significant increasing trends in the freezing season and in the thawing season. Therefore, the model simulation effectively reproduced the observed long-term changes. In the freezing season, since there was no glacier or snow melting (see Table 2), the runoff was mainly subsurface flow (groundwater flow and lateral flow from the unsaturated zone). In the thawing season, as shown in Table 2, snowmelt runoff contributed approximately 14 % of the total runoff, whereas glacier runoff contributed only a small fraction of the total runoff (approximately 2.2 %). Rainfall runoff was the major component of the total runoff in the thawing season, and the runoff increase in the thawing season was mainly due to increased precipitation and snowmelt. As shown in Fig. 11, the actual evapotranspiration increased significantly in both seasons due to increased precipitation and ground warming. The increase in actual evapotranspiration was greater in the thawing season than in the freezing season.

Figure 12 shows the changes in the basin-averaged annual liquid soil water storage (0–3 m) and groundwater storage. The annual liquid soil water storage showed a significant increasing trend, especially in the most recent 3 decades. This long-term change in liquid water storage was similar to the runoff change in the freezing season, as shown in Fig. 11a, exhibiting a correlation coefficient of 0.79. The annual ice water storage in the top 0–3 m soil layers showed a significant decreasing trend due to frozen soil changes. Annual groundwater storage showed a significantly increasing trend especially in the most recent three decades, which indicates that the groundwater recharge has increased with permafrost degradation.

Basin-average annual water storage (equivalent water depth) changes simulated over the period of 1971 to 2013 for (a) liquid water in the top layer of the ground (0–3 m), (b) ice in the top layer of the ground (0–3 m), (c) and groundwater. Results from linear regressions are indicated.

[Figure omitted. See PDF]

Discussion

Impact of frozen soil changes on the soil moisture and runoff

We have plotted the long-term changes in the spatially averaged liquid soil moisture in the areas with elevations between 3300 and 3500 m and in the areas with elevations between 3500 and 3700 m in Fig. S6 in the Supplement. In the seasonally frozen ground at elevations of 3300–3500 m, the liquid soil moisture increased slightly due to the decrease in the frozen depth (Fig. 10a). At elevations of 3500–3700 m, the liquid soil moisture in the deep layer increased significantly since the 1990s, due to the change of the permafrost into seasonally frozen ground (Fig. 10b).

In the freezing season, since the surface ground is frozen, runoff is mainly subsurface flow coming from seasonally frozen ground. Runoff has the highest correlation ($r=\mathrm{0.82}$) with the liquid soil moisture in the freezing season, which indicates that the frozen soil changes were the primary cause of the increase in liquid soil moisture, resulting in increased runoff in the freezing season. During the past 40 years, parts of permafrost changed into seasonally frozen ground and the frozen depth of the seasonally frozen ground decreased, leading to increases in the liquid soil moisture in the deep layers during the freezing season. The increase in liquid soil moisture also increased the hydraulic conductivity, which enhanced the subsurface flow. Figure 13c shows the seasonal pattern of runoff from the entire basin. From April to October (the thawing season), runoff in the permafrost area was much larger than in the area of seasonally frozen ground. However, in the freezing season the inverse was true. Figure S7 in the Supplement shows runoff changes from a typical area (with elevations of 3500–3700 m) that featured permafrost during the period of 1971 to 1980, but degraded to seasonally frozen ground during the period of 2001 to 2010. This illustrates that the thawing of permafrost increased runoff in the freezing season and slowed hydrological recession processes in autumn. Figure S4 illustrates the increase in freezing season runoff and the shift in the seasonal flow patterns simulated by the model without the frozen soil scheme.

Figure 13 shows the large difference in runoff variation with elevation between the freezing and thawing seasons. In the freezing season, the runoff change from the 1970s to the 2000s in the areas of seasonally frozen ground (mainly located below 3500 m, see Fig. 9) was relatively small. The areas with elevations of 3500–3900 m showed large changes in runoff. This pattern is due to the shift from permafrost to seasonally frozen ground in some areas in the elevation range of 3500 to 3900 m, as simulated by the model, particularly for sunny hillslopes (see Fig. 9). This finding illustrates that a change from permafrost to seasonally frozen ground has a larger impact on the runoff than a change in frozen depth in areas of seasonally frozen ground. In the thawing season, runoff increased with elevation due to the increase in precipitation with increasing elevation, and the magnitude of the runoff increase was mainly determined by magnitude of the precipitation increase (Gao et al., 2016). Precipitation in the region with elevations below 3100 m was low, and the air temperature was high. Hence, runoff in this region was lower during 2001–2010 than during 1971–1980 because of greater evapotranspiration.

Model-simulated runoff showing changes from the 1971–1980 period to the 2001–2010 period with elevation for (a) the freezing season and (b) the thawing season, and (c) monthly averaged seasonal runoff in permafrost and seasonally frozen ground for the period of 2001–2010.

[Figure omitted. See PDF]

Comparison with the previous similar studies

In this study, the model simulation showed that the thawing of frozen soils led to increased freezing season runoff and base flow in the upper Heihe basin. This result is consistent with previous findings based on observations in high-latitude regions (Walvoord and Striegl, 2007; St. Jacques and Sauchyn, 2009; Ye et al., 2009) and in northeastern China (Liu et al., 2003). However, those studies did not consider spatial variability. This study found that the impact of the frozen soil thawing on runoff varied regionally. In the upper Heihe basin (see Fig. 13), the change in the freezing season runoff was strongly affected by the change from permafrost to seasonally frozen ground in the higher-elevation region and by the evaporation increase in the lower-elevation region due to rising air temperature. However, runoff at the basin scale mainly came from the higher-elevation regions.

This study also showed that the thawing of frozen soils increased the liquid soil moisture in the upper Heihe basin, which is consistent with the finding of Subin et al. (2013) using the CLM to simulate northern high-latitude permafrost regions, and the findings of Cuo et al. (2015) using the VIC model to simulate 13 sites on the QTP. In contrast, Lawrence et al. (2015) found that permafrost thawing reduced soil moisture based on CLM simulations of the global permafrost region. This finding might be related to the uncertainties in the soil water parameters and the high spatial heterogeneity of soil properties, which are difficult to consider in a global-scale model. Subin et al. (2013) and Lawrence et al. (2015) simulated the soil moisture changes in the active layer of permafrost over large areas with coarse spatial resolution. Unlike those studies, this study investigated the spatio-temporal variability in soil moisture using a high spatial resolution and analyzed the impacts of frozen soil changes.

Jin et al. (2009) found decreased soil moisture and runoff due to permafrost degradation based on observations at the plot scale in the source area of the Yellow River basin. These results are different from those in the present study, possibly due to the difference in the hydrogeological structure and soil hydraulic parameters between the source area of the Yellow River and the upper Heihe basin. Wang et al. (2015a) estimated the increasing trend of the maximum frozen depth in the seasonally frozen ground to be 0.04 m decade${}^{-\mathrm{1}}$ during 1972–2006 in the Heihe River basin based on plot observations, which is consistent with the results in this study. The increase in groundwater storage illustrated in this study is also consistent with the findings of Cao et al. (2012) based on Gravity Recovery And Climate Experiment (GRACE) data, which showed that groundwater storage increased during the period of 2003–2008 in the upper Heihe basin.

Uncertainty in simulation of the frozen soils

Estimation of the change in permafrost area is a great challenge due to such complex factors as climatology, vegetation, and geology. Guo and Wang (2013) reported that the permafrost area for the whole QTP decreased from approximately $\mathrm{175.0}\times {\mathrm{10}}^{\mathrm{4}}$ km${}^{\mathrm{2}}$ in 1981 to $\mathrm{151.5}\times {\mathrm{10}}^{\mathrm{4}}$ km${}^{\mathrm{2}}$ in 2010, with a relative change of 13.4 %. Wu et al. (2005) reported that the permafrost area decreased by 12 % from 1975 to 2002 in the Xidatan basin of the QTP based on a ground penetration radar survey. Jin et al. (2006) found an area reduction of 35.6 % in island permafrost in Liangdaohe, which is located along the southern portion of the Qinghai–Tibet Highway, from 1975 to 1996. Compared with the borehole observations by Wang et al. (2013) shown in Fig. 2, our model slightly overestimated the soil temperature in permafrost areas, possibly leading to an overestimation of the rate of permafrost area reduction.

There were two major uncertainties in the frozen soil simulation: uncertainty in the simulation of the land surface energy balance and uncertainty in the simulation of the soil heat–water transfer processes (Wu et al., 2016). Uncertainty in the land surface energy balance simulation might result from uncertainty in the radiation and surface albedo estimates due to the complex topography, vegetation cover, and soil moisture distribution, thereby introducing uncertainties into the estimated ground temperature and soil heat flux. The uncertainty in the simulation of soil heat–water transfer processes might result from the soil water and heat parameters and the bottom boundary conditions of heat flux. For example, the soil depth and the fraction of rock in the soil can greatly affect the ground temperature simulation. Permafrost degradation is closely related to the thermal properties of rocks and soils, the geothermal flow, the initial ground temperature, and soil ice conditions. Sub-grid topography may also affect the frozen soil simulation. For example, active layer thickness is different between the low-elevation valleys and higher-elevation slopes due to the temperature inversion caused by the accumulation of cold air in valleys (Bonnaventure et al., 2012; Zhang et al., 2013; O'Neill et al., 2015). In areas with high groundwater flow rates, laterally advected heat flux may increase the thawing of permafrost (Kurylyk et al., 2016; Sjöberg et al., 2016). Not considering the lateral heat flux may lead to an underestimation of talik development and thawing rates of permafrost. In addition, uncertainties in the input data, particularly solar radiation (which is estimated using interpolated sunshine hour data from a limited number of observational stations) and precipitation (which is also interpolated based on observations at these stations), may also influence the results of the model simulation. Due to the complexity of the distributed model and the large number of model parameters, quantifying the overall simulation uncertainty is challenging, but is part of our ongoing research.

Conclusions

This work carefully validated a distributed hydrological model coupled with cryospheric processes in the upper Heihe River basin using available observations of soil moisture, soil temperature, frozen depth, actual evaporation, and streamflow discharge. Based on the model simulations from 1971 to 2013 in the upper Heihe River, the long-term changes in frozen soils were investigated, and the effects of the frozen soil changes on the hydrological processes were explored. Based on these analyses, we have reached the following conclusions:

• The model simulation suggests that 8.8 % of the permafrost areas degraded into seasonally frozen ground in the upper Heihe River basin during 1971–2013, predominantly between elevations of 3500 and 3900 m. The results indicate that the decreasing trend of the annual maximum frozen depth of the seasonally frozen ground is 0.032 m decade${}^{-\mathrm{1}}$, which is consistent with previous observation-based studies at the plot scale. Additionally, our work indicates a trend of increasing active layer thickness in the permafrost regions of 0.043 m decade${}^{-\mathrm{1}}$.

• The model-simulated runoff trends agree with the observed trends. In the freezing season (November–March), based on the model simulation, runoff was mainly sourced from subsurface flow, which increased significantly in the higher elevation regions where significant frozen soil changes occurred. This finding implies that the runoff increase in the freezing season is primarily caused by frozen soil changes (permafrost degradation and reduced seasonally frozen depth). In the thawing season (April–October), the model simulation indicates that runoff was mainly sourced from rainfall and showed an increasing trend at higher elevations, which can be explained by the increase in precipitation. In both the freezing and thawing seasons, the model-simulated runoff decreased in the lower-elevation regions, which can be explained by increased evaporation due to rising air temperatures.

• The model-simulated changes in soil moisture and ground temperature indicate that the annual storage of liquid water increased, especially in the most recent three decades, due to frozen soil changes. The annual ice water storage in the top 0–3 m of soil showed a significant decreasing trend due to soil warming. The model-simulated annual groundwater storage had an increasing trend, which is consistent with the changes observed by the GRACE satellite. Therefore, groundwater recharge in the upper Heihe basin has increased in recent decades.

• The model simulation indicated that regions where permafrost changed into seasonally frozen ground had larger changes in runoff and soil moisture than the areas covered by seasonally frozen ground throughout the study period.

For a better understanding of the changes in frozen soils and their impact on ecohydrology, the interactions among soil freezing–thawing processes, vegetation dynamics, and hydrological processes need to be investigated in future studies. There are uncertainties in simulations of frozen soils and hydrological processes that also warrant further investigation.

All data cited in this paper are available from the references. The model code with a working example is freely available online (https://github.com/gb03/GBEHM, Gao, 2017) or upon request from the corresponding author ([email protected]).

The Supplement related to this article is available online at https://doi.org/10.5194/tc-12-657-2018-supplement.

The authors declare that they have no conflict of interest.

Acknowledgements

This research was supported by the major plan of “Integrated Research on the Ecohydrological Processes of the Heihe basin” (project nos. 91225302 and 91425303) funded by the National Natural Science Foundation of China (NSFC). The authors would like to thank the editor and reviewers for their constructive suggestions, which have substantially improved the paper. Edited by: Peter Morse Reviewed by: two anonymous referees