Abstract

The Water and Energy transfer Processes in Cold Regions (WEP-COR) model is an improved version of the Water and Energy transfer Processes in Large basins (WEP-L) model that integrates a multi-layer frozen soil model to simulate the hydrological processes in cold regions and the heat fluxes at different depths of frozen soil. The temperature, water content, freezing depth of the soil, and daily discharge were simulated and compared with observations. The simulated and observed data were used to analyze the runoff flow components. The results showed that the WEP-COR model can effectively simulate the distributions of the soil temperature and water content. The average root mean squared errors of the temperature, unfrozen water content, total water content, and freezing depth of the soil were 1.21 °C, 0.035 cm3/cm3, 0.034 cm3/cm3, and 17.6 cm, respectively. The mean Nash–Sutcliffe efficiency and relative error of the daily discharge were 0.64 and 6.58%, respectively. Compared with the WEP-L model, the WEP-COR model simulated the discharge with higher accuracy, especially during the soil thawing period. This improvement was mainly due to the addition of the frozen soil mechanism. The WEP-COR model can provide support for agricultural and water resources management in cold regions.

INTRODUCTION

The Qinghai-Tibet plateau, northwestern alpine area, and Northeastern China are the three major cold regions of China (Chen et al. 2006). They not only cover about 43.5% of the total land area but are also the primary headsprings of the water supply for the arid and semiarid regions of China (Kang et al. 1999; Chen et al. 2008a; Liu & Yao 2016). Due to their high latitude and high altitude, these regions have a wide distribution of glaciers, snow cover, permafrost, and seasonal frozen ground (Jiang et al. 2009; Homan et al. 2011; Wang et al. 2015). The impermeability of frozen soil changes its capacity for surface storage, infiltration, and evaporation; this influences the hydrological cycle processes of the land surface (Gusev & Nasonova 2002; Yamazaki et al. 2006; Tian et al. 2009). Meanwhile, global warming has led to shrinking glaciers and the degeneration of permafrost (Liu et al. 2003; Smith et al. 2004; Cheng & Wu 2007; Niu et al. 2016). These changes alter the runoff generation processes and runoff amounts of these regions, which in turn affect their water supply capacity (Barnett et al. 2005; Han et al. 2014; Zhang et al. 2015). Northeastern China is an important base for producing grain as a commodity. The freezing and thawing (FT) processes of seasonal frozen soil in this region affect changes in the water phase and are vital to crop cultivation (Hejduk & Kasprzak 2010). Therefore, exploring the hydrological cycle under frozen soil conditions is very important for water resources management and agricultural production in cold regions.

In the past decades, various models have been developed to simulate the transfer of water and heat fluxes in frozen soil and improve the modeling parameterization for frozen soil (Harlan 1973; Flerchinger et al. 1998; Jansson & Moon 2001; Li et al. 2010). However, these models only treat one-dimensional water flows (Takata 2002). They only consider the vertical migration of water without considering the lateral flow, so they cannot calculate the runoff and overland flow in cold regions. Therefore, some researchers have tried coupling frozen soil mechanisms to land surface models, such as the Variable Infiltration Capacity (VIC) model (Liang et al. 1994; Cherkauer et al. 2003; Zhao et al. 2013) and Community Land Model (CLM) (Niu & Yang 2006; Wang et al. 2014). In addition, the soil FT processes have been incorporated into distributed hydrological models, such as the Cold Regions Hydrological Model (CHRM) platform (Pomeroy et al. 2007; Zhou et al. 2014), Soil and Water Assessment Tool (SWAT), Distributed Water–Heat Coupled (DWHC) model (Chen et al. 2008a, 2008b), Water and Energy Budget-based Distributed Hydrological Model (WEB-DHM) (Shrestha et al. 2010), and Geomorphology Based Hydrological Model (GBHM) (Zhang et al. 2013). However, with the exception of CRHM and the DWHC model, most frozen soil modules in hydrological models are simplified and empirical. Although CRHM contains various modules for hydrological processes, it cannot simulate the distributions of soil temperature and water content at different depths and does not consider the land use type. The DWHC model was set up by referencing SWAT and Coup Model (Jansson & Moon 2001). Because the calculation of some thermal parameters depends on the results of other parameters, this correlation increases the uncertainty of the model. In addition, the use of many parameters increases the difficulty of modeling hydrological processes in a cold region, which lacks detailed hydrological and soil data. Therefore, a DWHC model that considers the land use type and frozen soil with some conceptual parameters to compensate for the lack of detailed data should be developed. The Water and Energy Transfer processes in Large basins (WEP-L) model is a physicallybased and distributed hydrological model. Each calculation unit of the model considers land use heterogeneity by adopting the mosaic method (Jia et al. 2006). However, the WEP-L model does not contain a module to simulate soil FT processes. Thus, the focus of this study was to couple a frozen soil model to the WEP-L model and improve the simulation performance for basins in cold regions. The improved Water and Energy Processes in Cold Regions (WEP-COR) model can consider the impact of soil FT processes on hydrological processes in cold regions. The WEP-COR model was applied to simulating the soil heat and water fluxes, soil temperature, freezing depth, and runoff of the Second Songhua River (SSR) basin.

METHODOLOGY

Description of the WEP-L model

The WEP-L model has been successfully applied to several basins in Japan, Korea, and China (Jia & Tamai 1998; Jia et al. 2009). It combines the merits of physically based spatially distributed (PBSD) models (Refsgaard et al. 1996) and the soil–vegetation–atmosphere transfer (SVAT) scheme to represent spatially variable water and energy processes in soils with complex land covers. To consider the impact of the topography and land cover on the water cycle in large basins, the spatial calculation unit of the WEP-L model is contour bands inside small sub-basins. Each unit in the vertical direction includes five layers from top to bottom: the interception layer, depression layer, soil layer, transition layer, and aquifer (Figure 1). In addition, each unit in the horizontal direction is classified into five classes with the mosaic method: water body, soil–vegetation, irrigated farmland, non-irrigated farmland, and impervious area (Avissar & Pielke 1989). The soil layers of the soil–vegetation and farmland classes are further divided into three layers to describe soil evaporation and the water uptake of vegetation roots. The average water and heat flux is obtained by areally averaging those from each land use in a contour band.

Figure 1

Vertical structure of the WEP-COR model.

Figure 1

Vertical structure of the WEP-COR model.

The time step of the WEP-L model is generally 1 day. Evapotranspiration from the water body and soil is calculated by using the Penman equation, while evapotranspiration with vegetation canopies is calculated by using the Penman–Monteith equation (Monteith 1973). The canopy resistance (Noilhan & Planton 1989) is related to the soil moisture condition. The infiltration and surface runoff during rainfall greater than 10 mm/d are calculated with a generalized Green–Ampt model (Jia & Tamai 1998). The soil moisture movement in unsaturated soils is calculated with the Richards model. The air temperature is used to adjust the calculation of the saturated hydraulic conductivity in frozen soil (Jia et al. 2009). The snow melt is calculated with the temperature index approach. The subsurface runoff is generated according to the land slope and soil hydraulic conductivity. The ground water flow is calculated by using the Boussinesq equations (Zaradny 1993). The groundwater outflow is calculated according to the hydraulic conductivity of the riverbed material and the difference between the river water stage and groundwater level. The overland flow and river flow are calculated by using the kinematic wave method (Jia et al. 2001). The model approaches for water and energy processes are described in Jia et al. (2001).

Model improvement: frozen soil scheme

Seasonal frozen soil commonly occurs in Northeast China. Phase changes of the soil water have a considerable influence on the infiltration and evaporation (Konrad & Duquennoi 1993). In addition, the soil water flux and runoff generation process in a cold region are strongly dependent on the soil freezing depth (Hayashi et al. 2007). The heat fluxes at different soil depths (i.e., soil layers) can change the soil water phase and soil temperature, while the soil moisture content and temperature difference between soil layers determine the heat flux. Therefore, the soil temperature and soil solid water need to be represented in a model through an energy balance. Although the WEP-L model can be used to simulate the soil surface temperature and heat flux of soil, the simulation results cannot reflect the heat flux transfer in soil layers. Besides, the WEP-L model divides the upper soil (2 m) into three layers, so only three average values of the soil moisture content can be shown at three depths. The distributions of the soil moisture and temperature in a soil group cannot be revealed either. To overcome these deficiencies, we added the calculation of the heat transfer and water phase change for different soil layers. The improved model can be used to simulate the soil temperature, soil solid water content, and freezing depth. The depth of a soil layer is one of the parameters for calculating the heat flux. Therefore, in consideration of the computation time, the WEP-COR model divides the upper soil (2 m) into 11 layers. The top two layers are set to a depth of 10 cm because the surface soil is sensitive to climate changes, and the other layers are each set to a standard depth of 20 cm. The number of layers can be adjusted if the soil thickness is less than 2 m. The calculations for the soil heat transfer and soil FT processes were added to all three land use classes (i.e., soil–vegetation, irrigated farmland, and non-irrigated farmland). Figure 1 shows the vertical structure of the WEP-COR model. The methods for the coupled soil water–heat processes were mainly derived from the Coup model and Shang's model (Shang et al. 1997; Wang et al. 2014). These models can clearly represent the water and heat flux transfer of soil FT processes and parameters for different soil status. Thus, they are coupled in the WEP-COR model. The water–heat continuous equation of frozen soil is solved numerically based on the soil freezing status and empirical formulas.

Soil freezing status

The soil freezing status is divided into three types. (1) For the unfrozen type, the soil temperature (Ts) is above 0 °C, and there is no solid water in the soil layer. (2) For the frozen type, Ts is lower than Tf (i.e., the threshold temperature value). The soil is assumed to be completely frozen with a residual unfrozen amount. According to measurements, when the soil temperature is less than −10 °C, the liquid water content of the soil is stable at about 0.09 cm3/cm3 (Wu et al. 2013). Here, the Tf value was revised to −10 °C. (3) For the partially frozen type, Ts is higher than Tf but lower than 0 °C. Liquid water and solid water can coexist in the soil.

Boundary condition and heat flux into soil

The WEP-COR model assumes that the upper boundary of the soil group is the atmosphere, which controls the input and output of the energy. The upper boundary energy transfer can be calculated from meteorological variables including the temperature, wind speed, hours of sunshine, and relative humidity. The bottom boundary of the soil group is the transition layer or aquifer, which has a constant temperature. The energy balance equation on land surface is expressed as (Jia et al. 2001): 
formula
(1)
where RN (J/m2/d) is the net radiation, Ae (J/m2/d) is the anthropogenic energy source, lE (J/m2/d) is the latent heat flux, H (J/m2/d) is the sensible heat flux, and G (J/m2/d) is the heat conduction into soil. The force-restore method (FRD) (Hu & Islam 1995) is used to solve G and the surface temperature of different land covers.

Soil heat flux transfer

According to the law of energy conservation, the equation of the soil vertical heat flux transfer can be written as (Shang et al. 1997; Wang et al. 2014): 
formula
(2)
where z is the soil depth (m) that represents each soil layer, λs is the soil thermal conductivity (W/(m•°C)), Ts is the soil temperature (°C), CV is the composite soil heat capacity (J/(m3•°C)), t is the time (s), ρi is the ice density (920 kg/m3), Li is the latent heat of fusion (3.35 × 106 J/kg), and θi is the soil volumetric ice content. Equation (2) describes the heat transfer among soil layers at different depths and the changes in the soil temperature and water phase. We can use the numerical iterative method to solve Equation (2). Then, the finite difference scheme can be written as: 
formula
(3)
 
formula
(4)
 
formula
(5)
 
formula
(6)
 
formula
(7)
where j is the number of soil layers and k is the time. If the time step is 1 day, k + 1 represents the day after. The other variables are the same as those defined previously.
The composite soil heat capacity CV is expressed as (Jansson & Moon 2001): 
formula
(8)
where θs is the soil saturated water content, θl is the soil liquid water content, and Cs, Cl, and Ci are the heat capacities of soil, water, and ice, respectively.
The thermal conductivity is a complex function of the soil moisture and constituents. The Coup model and other models calculate the thermal conductivity by using different equations depending on the soil freezing status; however, this approach requires the determination of several parameters. In consideration of the parameter determination, here the soil thermal conductivity was calculated by using the IBIS model (Foley et al. 1996): 
formula
 
formula
(10)
where λst is the dry soil thermal conductivity (W/(m•°C)) and ωsand, ωsilt, and ωclay represent the volumetric fractions of sand, silt, and clay, respectively.

Soil temperature

The energy flux drives the changes in the soil temperature and water phase, while the soil temperatures of different soil layers impact the soil sensitive heat. The temperature of the top soil layer is calculated with the FRD method. For other soil layers, the temperature at the middle part of a soil layer is used to represent its average temperature. The soil temperature and heat flux between adjacent soil layers can be calculated as follows (Chen et al. 2008b): 
formula
(11)
 
formula
(12)
where Hi,i+1 is the sensible heat flow between soil layers i and i + 1. The initial Ts of each soil layer is the input of the model. The soil temperature and moisture content are simulated by numerical iteration.

Soil water flux transfer

During soil FT periods, migration only occurs in liquid water. The soil vertical heat flux transfer can be written as (Shang et al. 1997; Wang et al. 2014): 
formula
(13)
where D(θl) and K(θl) are the hydraulic diffusivity and hydraulic conductivity, respectively, for unsaturated soil and ρl is the density of water (kg/m3). K(θl) is closely related to the saturated hydraulic conductivity Ks, which is corrected by the soil temperature (Jansson & Moon 2001): 
formula
(14)
 
formula
(15)
where θr is the soil residual moisture content, n is the Mualem constant, Ks is the saturated hydraulic conductivity of the soil temperature correction (cm/s), and K0 is the initial saturated hydraulic conductivity (cm/s).

Modeling procedure of soil heat and water transfer

The WEP-COR model couples the calculation of the soil heat and water transfer processes with the WEP-L model. Figure 2 shows the simulation flowchart of the soil heat and water transfer. Both vertical and lateral flows are calculated. In these charts, R is the soil lateral flow, E is the evaporation of vegetation and bare soil, Q is the soil gravity drainage, and QD is the water flow between adjacent soil layers. The WEP-COR model iteratively calculates the finite difference scheme to compute the soil temperature and soil water flow. The time index represents the time, and the change interval of the time index is equal to the time step. Tt+1 represents T the next day. The cycle index represents the iterations. There are two parts to the iterative computations. First, the water and heat equations of the frozen soil are solved following Equations (3)–(7). Then, the soil moisture migration caused by evaporation, infiltration, and gravity drainage is calculated. n is the iterations of the former, and m is the iterations of the latter. An error of within 0.001 is acceptable; the iterative calculation does not stop until the error value is acceptable.

Figure 2

Flowchart of the soil heat and water transfer.

Figure 2

Flowchart of the soil heat and water transfer.

STUDY SITES AND DATA

Study area

The WEP-COR model was applied to the SSR basin. The SSR is a tributary of the Songhua River in Northeast China and covers an area of 7.4 × 104 km2 (Figure 3). Observations from 1971 to 1995 indicated that the long-term average annual air temperature is approximately 4.2 °C and the average annual precipitation is 700 mm for the SSR basin. The basin is located in a typical cold region where seasonal frozen soil is common. The soil FT period of this basin is usually from November to May. The water and heat flux simulations of the WEP-COR model were evaluated by taking measurements from the Qianguo irrigation experimental station (124°30′30″E, 45°14′00″N) and a meteorological station (station 54266, 125°37′59″E, 42°31′59″N). To consider the integrity of the data, the Wudaogou station (126°37′43″E, 42°53′51″N) was selected to evaluate the modeling performance for the daily discharge. This station is located upstream of the SSR.

Figure 3

Locations of the Second Songhua River basin and stations.

Figure 3

Locations of the Second Songhua River basin and stations.

Data description

Figure 3 shows the distribution of the rivers and main hydro-meteorological stations in the SSR basin. The WEP-COR model requires data on the geography (elevation, vegetation), land use, meteorology, and parameters (soil type, hydraulic conductivity, and soil moisture characteristic). Elevation data were obtained from SRTM90. Monthly vegetation information (LAI and area fraction of vegetation) was obtained from the NOAA-AVHRR data. The soil data (soil type and corresponding characteristic parameters) were acquired from the National Second Soil Survey Data and Soil Types of China. Land use data were obtained from the Landsat TM data and statistical data in the yearbooks of administrative districts. Meteorological daily data including the precipitation, temperature, wind speed, sunshine, and humidity can be downloaded from the website http://cdc.cma.gov.cn. Based on the DEM data and observed river information, the SSR basin was divided into 1,305 sub-basins, each of which was assigned by using the Pfafstetter code (Verdin & Verdin 1999). The meteorological daily data were interpolated to each sub-basin by using the inverse distance weighted method.

Parameter sensitivity analysis of the WEP-L model was previously done by Jia et al. (2006). The conductivity of riverbed material, soil layer thickness, maximum soil moisture content, and groundwater aquifer hydraulic conductivity were identified as parameters with high sensitivity. Most parameters in the model do not need to be calibrated, but the high-sensitivity parameters were adjusted by comparing the simulated discharge with observed values during the selected calibration period. The soil physical properties presented in Table 1 were measured at the Qianguo experiment station. Based on the texture information, the soils were reclassified into four categories: sand, loam, clay loam, and clay. The main soil moisture characteristics in Table 2 were taken from Jia et al. (2006). The conductivity of the riverbed was set to 1.728 m/d. The groundwater aquifer hydraulic conductivity and specific yield in the basin were deduced from groundwater simulation and geological exploration data. The hydraulic conductivity was set to 1.056 m/d, and the specific yield was set to 0.05 m/d.

Table 1

Some soil physical properties and the permeability coefficient

Soil depth (cm) Particle size distribution (%)
 
Bulk weight (g·cm−3Permeability coefficient (10−4 cm·s−1
<2 μm 2–50 μm >50 μm 
0–15 28.0 41.0 31.0 1.2 3.24 
15–28 30.5 35.8 33.7 1.4 1.25 
28–100 18.5 28.9 52.6 1.5 2.81 
Soil depth (cm) Particle size distribution (%)
 
Bulk weight (g·cm−3Permeability coefficient (10−4 cm·s−1
<2 μm 2–50 μm >50 μm 
0–15 28.0 41.0 31.0 1.2 3.24 
15–28 30.5 35.8 33.7 1.4 1.25 
28–100 18.5 28.9 52.6 1.5 2.81 
Table 2

Soil moisture characteristics

Parameters Sand Loam Clay loam Clay 
θs 0.4 0.466 0.475 0.479 
θf 0.174 0.278 0.365 0.387 
θr 0.077 0.120 0.170 0.250 
SW 6.1 8.9 12.5 17.5 
n 3.37 3.97 3.97 4.38 
Parameters Sand Loam Clay loam Clay 
θs 0.4 0.466 0.475 0.479 
θf 0.174 0.278 0.365 0.387 
θr 0.077 0.120 0.170 0.250 
SW 6.1 8.9 12.5 17.5 
n 3.37 3.97 3.97 4.38 

θf: field capacity; SW: suction at wet front of soils; n: constant parameter of Mualem.

Methods of evaluation

Experimental data from October 2011 to May 2012 were used to evaluate the simulation results of the soil FT processes at the Qianguo experiment station. The daily freezing depth of station 54266 and daily discharge data of the Wudaogou station were split into two parts: data from 1971 to 1985 for calibration and data from 1986 to 1995 for validation. The modeling performance was statistically evaluated by a qualitative assessment though graphs first and then a quantitative assessment with statistical measures. The root mean squared error (RMSE), Nash–Sutcliffe coefficient (NSE), and relative error (RE) were used for the quantitative evaluation. The RMSE, NSE, and RE can be calculated as follows: 
formula
(16)
 
formula
(17)
 
formula
(18)
where n is the number of observations, Oi is the observed value, is the mean observed value, and Si is the simulated value. For the NSE, the best value is 1, and a negative value means that the model is not credible.

RESULTS AND DISCUSSION

Soil temperature simulation

Figure 4 shows the simulated and observed daily mean soil temperatures at different depths from 24 October 2011 to 11 May 2012 at the Qianguo experiment station. The temperature of the upper layers (0–20 cm) fluctuated wildly, and the soil temperature and its fluctuation diminished with depth. This result agrees with other studies (Guo et al. 2011; Xiang et al. 2013). The simulated soil temperature of the upper soil layers matched the observations well, but the temperatures of the middle and lower layers were slightly underestimated during the thawing period. This underestimation may have been due to the homogenized soil hydraulic conductivity and thermal parameters in the experiment (Table 1); the parameters were the same from the third to 11th layers. In addition, the minimum simulated soil temperature was lower than that observed. This previous error may have resulted in the later underestimation during the thawing period. The soil thermal conductivity may also be greater than the simulation, so the simulated soil temperature rose slower than that observed during the thawing period. Table 3 presents the NSE and RMSE of the soil temperature at different depths. The layer at 120 cm had missing data, which led to a higher NSE value compared to other layers because only the freezing period was represented and not the FT periods. The mean NSE was 0.92, and the mean RMSE was 1.21 °C. Overall, the soil temperature simulated with the WEP-COR model was similar to the observed data.

Figure 4

Simulated and observed soil temperatures at different depths.

Figure 4

Simulated and observed soil temperatures at different depths.

Table 3

Statistical values from the soil temperature simulation

Depth (cm) 10 20 35 60 90 100 110 120 
NSE 0.95 0.97 0.90 0.86 0.91 0.92 0.90 0.98 
RMSE 1.47 0.94 1.17 1.23 1.14 1.27 1.48 1.04 
Depth (cm) 10 20 35 60 90 100 110 120 
NSE 0.95 0.97 0.90 0.86 0.91 0.92 0.90 0.98 
RMSE 1.47 0.94 1.17 1.23 1.14 1.27 1.48 1.04 

Soil moisture simulation

Figure 5 illustrates the simulated and measured volumetric liquid water and total water (both unfrozen and frozen moisture) of the soil at different depths during the FT periods from 2011 to 2012. The liquid water content was equal to the total water content on 2 November, which means that there was no solid water at all and the soil was unfrozen. As the air temperature decreased, the soil froze from top to bottom (Figure 5(b) and 5(c)) (Slater et al. 1998). At the early stage of soil freezing, solid water appeared in the upper soil layers. As the liquid water shifted from unfrozen soil to frozen soil, the liquid water content of the top layer decreased, and its total water content increased (Iwata & Hirota 2005; Sheshukov & Nieber 2011). Figure 5(b) and 5(c) show that the simulated total water content was less than that measured, which may have resulted from the higher simulated values of evaporation. The simulated liquid water content was higher than that measured, which was probably caused by the differences in soil temperature. On 22 February, the soil water distribution showed a ‘V’ shape, which indicates that the soil frozen depth was almost at the maximum (Cheng & Wu 2007; Hayashi et al. 2007). Figure 5(e) shows an ‘O’ shape for the water distribution, which means that both the upper and bottom soil layers were thawing. However, the simulated thaw rate was less than that measured. The heat flux of the middle layer was overestimated, so the liquid water content was higher than that measured. Figure 5(f) shows that the total water content was equal to the liquid water content, which indicates that the soil thawing process was completed. Table 4 presents the statistical values, and the mean RMSEs for the unfrozen water content and total water content were 0.035 and 0.034 cm3/cm3, respectively.

Figure 5

Simulated and observed soil moisture contents at different depths.

Figure 5

Simulated and observed soil moisture contents at different depths.

Table 4

Statistical values from the soil water content simulation

Depth (cm) 0–10 10–20 20–40 40–60 60–80 80–100 100–120 120–140 140–160 
Unfrozen 0.030 0.048 0.031 0.029 0.028 0.024 0.016 0.041 0.069 
Total 0.031 0.025 0.027 0.048 0.026 0.017 0.014 0.043 0.076 
Depth (cm) 0–10 10–20 20–40 40–60 60–80 80–100 100–120 120–140 140–160 
Unfrozen 0.030 0.048 0.031 0.029 0.028 0.024 0.016 0.041 0.069 
Total 0.031 0.025 0.027 0.048 0.026 0.017 0.014 0.043 0.076 

Unfrozen: unfrozen water content; total: sum of unfrozen and frozen water contents.

Soil freezing depth simulation

Figure 6 shows the simulated and observed soil freezing depths during the FT periods from 2011 to 2012. As the temperature decreased, the soil began to freeze in mid-November and reached the maximum frozen depth in early March. When heat from the atmosphere was less than the energy required to keep the basal freezing state, the frozen soil began to thaw on the surface and at the bottom (Cherkauer & Lettenmaier 1999; Woo et al. 2004), and the thawing process finished on 25 April. The simulated freezing depth matched the measurement during the freezing period but was a little greater than that measured during the thawing period. The deepest simulated frozen depth was 8 cm deeper than the measurement. The simulated FT periods were 10 d shorter than that measured. The average RMSE of the simulated freezing depth was 17.68 cm.

Figure 6

Simulated and observed soil freezing depths during the FT period.

Figure 6

Simulated and observed soil freezing depths during the FT period.

The soil frozen depth affects both the lateral and vertical soil water fluxes in cold regions (Hayashi et al. 2007). The long-term observations of the daily freezing depth at station 54266 were used to evaluate the modeling performance. Figure 7 shows the daily simulated and observed soil freezing depths from 1971 to 1995. The simulated freezing depth generally matched the observed result well. The NSE and RMSE of the calibration were 0.95 and 11.3 cm, respectively. The NSE and RMSE of the validation period were 0.92 and 13.2 cm, respectively.

Figure 7

Daily simulated and observed soil freezing depths at station 54266. (a) Calibration (1971–1985). (b) Validation (1986–1995).

Figure 7

Daily simulated and observed soil freezing depths at station 54266. (a) Calibration (1971–1985). (b) Validation (1986–1995).

Water discharge simulation

Figure 8 presents the simulation results of the WEP-L and WEP-COR models for the daily discharge, and Table 5 presents the statistical test. The graphs indicate that the variation tendencies of the simulations were consistent with the observations. However, the WEP-L model without the frozen soil scheme systematically underestimated the stream flow. The NSE and RE were 0.38 and −53.69%, respectively, for the calibration period. The WEP-L model performed better during the validation period than the calibration period with an NSE and RE of 0.64 and −41.62%, respectively. As frozen soil has a considerable impact on the surface storage capacity, hydrological modeling of a basin in a cold region must include a frozen soil scheme, and the model performance cannot be improved solely by parameter adjustment. The model performance was obviously improved when coupled with a frozen soil scheme. The WEP-COR model achieved an NSE and RE of 0.42 and −0.97%, respectively, for the calibration period. For the validation period, the NSE was similar to that of the WEP-L model, but the RE endpoint showed a 35.04% decrease. Figure 8 shows that the simulated daily discharge during the thawing period (February–May) tended to underestimate the observations, and the simulations were consistent with the highest observations. We inferred that the high RE values were mainly due to the underestimation of the runoff during the thawing period (from February to May). In general, the WEP-COR model demonstrated an acceptable performance for the SSR basin and achieved efficiency coefficients of NSE > 0.6 and RE < 10% for the validation period. The simulated discharge can be used for further analysis.

Figure 8

Observed and simulated daily discharges of the Wudaogou station from 1971 to 1995. (a) Calibration (1971–1985). (b) Validation (1986–1995).

Figure 8

Observed and simulated daily discharges of the Wudaogou station from 1971 to 1995. (a) Calibration (1971–1985). (b) Validation (1986–1995).

Table 5

Statistical values for the daily discharge simulation from 1971 to 1995

 Calibration (NSE/RE) Validation (NSE/RE) 
WEP-L 0.38/−53.69% 0.64/−41.62% 
WEP-COR 0.42/−0.97% 0.64/6.58% 
 Calibration (NSE/RE) Validation (NSE/RE) 
WEP-L 0.38/−53.69% 0.64/−41.62% 
WEP-COR 0.42/−0.97% 0.64/6.58% 

Analysis of flow components

Based on the analysis presented in the last section, we inferred that the improved performance of the WEP-COR model was mainly for the thawing period. To further clarify the improvement, the simulated daily river discharge from February to May was compared with observations for the validation period (Figure 9). Table 6 presents the statistical test. In addition, the flow components from February to May were calculated, and the monthly statistical results are shown in Figure 10. As shown in Figure 9, the simulated daily river discharge of the WEP-L model was lower than that observed, while the simulation result of the WEP-COR model showed good agreement with the observation. The RE of the WEP-L model was −73.33%. In contrast, the RE of the WEP-COR model was −6.26%. Although both models underestimated the daily river discharge, the WEP-COR model clearly showed an improved performance.

Figure 9

Observed and simulated discharges from February to May (1986–1995). (a) WEP-L. (b) WEP-COR.

Figure 9

Observed and simulated discharges from February to May (1986–1995). (a) WEP-L. (b) WEP-COR.

Table 6

Statistical values for the simulated daily discharge from February to May (1986–1995)

 WEP-L WEP-COR 
(NSE/RE) −0.06/−73.33% 0.33/−6.26% 
 WEP-L WEP-COR 
(NSE/RE) −0.06/−73.33% 0.33/−6.26% 
Figure 10

Components of runoff from February to May during the validation period. (a) Snowmelt runoff. (b) Rainfall runoff. (c) Infiltration. (d) Groudwater discharge (e) Recharge groundwater.

Figure 10

Components of runoff from February to May during the validation period. (a) Snowmelt runoff. (b) Rainfall runoff. (c) Infiltration. (d) Groudwater discharge (e) Recharge groundwater.

Figure 10 shows the monthly mean variations in flow components from February to May during the validation period. The flow components represent the sources of the river discharge, which include the surface flow (snowmelt runoff and rainfall runoff), subsurface flow, and base flow (groundwater discharge). As shown in Figure 10(a), snowmelt runoff occurred in March and April as the air temperature rose. There was no distinct difference between the simulations of the two models. As shown in Figure 10(b) and 10(c), there was slight rainfall runoff and infiltration in February as the soil was frozen. Subsequently, the rainfall runoff of the two models tended to increase, but the WEP-L model showed a lower rainfall runoff than the WEP-COR model. The results indicated that the frozen soil decreased the soil infiltration capacity. The differences in the snowmelt runoff, rainfall runoff, and infiltration between the two models were small. The main differences were from the groundwater discharge and recharge. In this study, the groundwater discharge represented the water flux from the groundwater or aquifer layer into the river, and the groundwater recharge represented the water flux from the soil layer to the groundwater or aquifer layer. Due to the impermeability of frozen soil, the soil water flux could not further transfer to the deep layer and formed an aquifer layer above the frozen soil layer. The flow above the frozen soil layer was classified as base flow. Therefore, the calculated groundwater level was higher, so the lateral flow was higher. The WEP-L model without the frozen soil scheme underestimated the base flow. Regardless of the occurrence of frozen soil, the water flux infiltrated into the deep layer, and the groundwater level was lower with hardly any lateral flow. These may be the main reasons for the improved results with the WEP-COR model.

CONCLUSIONS

The WEP-COR model was developed to improve the modeling performance of the WEP-L model by adding a soil heat–water coupled module to help simulate the land surface water and energy budgets in a cold region. In addition, the number of soil layers was increased from three to 11 to simulate the soil temperature and soil moisture content at different depths. The simulated soil thermal and moisture distributions were compared with measurements taken at Qianguo irrigation experimental station. The results showed that the WEP-COR model can be used to simulate the vertical distributions of the soil temperature, soil liquid and solid water contents, and soil freezing depth in a cold region. The mean RMSEs of the soil temperature, liquid moisture content, total water content, and freezing depth were 1.21 °C, 0.035 cm3/cm3, 0.034 cm3/cm3, and 17.6 cm, respectively. The average NSE of the soil temperature was 0.92. The simulated daily freezing depth was compared with observations at station 54266 from 1971 to 1995. The simulated freezing depth generally matched the observed result well. The WEP-COR model significantly improved the predicted discharge for the SSR basin, especially during the soil thawing period. The analysis of the daily discharge and runoff components demonstrated that frozen soil needs to be considered when modeling the hydrological processes in a cold region. The simulated results showed an obvious improvement in the model performance when it was coupled with the frozen soil scheme. When the observed and simulated daily discharges were compared, the NSE of the WEP-COR model was 0.64, and the RE was 6.58%. In conclusion, the developed WEP-COR model contributes to the qualitative evaluation and prediction of the spatial and temporal distributions of frozen soil and change in water resources. It can act as a reference for agricultural and water resource management in a cold region. It can also be used to explore the hydrothermal transfer and hydrological cycle response to climatic change.

ACKNOWLEDGEMENTS

This work is supported by the National Natural Science Foundation of China (51179203), and the National Science and Technology Major Project for Water Pollution Control and prevention (2008ZX07207-006, 2012ZX07201-006).

REFERENCES

REFERENCES
Barnett
T. P.
,
Adam
J. C.
&
Lettenmaier
D. P.
2005
Potential impacts of a warming climate on water availability in snow-dominated regions
.
Nature
438
(
7066
),
303
309
.
Chen
R. S.
,
Kang
E. S.
,
Ji
X. B.
,
Yang
J. P.
&
Yang
Y.
2006
Cold regions in China (SCI)
.
Cold Regions Science & Technology
45
(
2
),
95
102
.
Chen
R. S.
,
Lu
S. H.
,
Kang
E. S.
,
Ji
X. B.
,
Zhang
Z. H.
,
Yang
Y.
&
Qing
W. W.
2008a
A distributed water–heat coupled model for mountainous watershed of an inland river basin of Northwest China (I) model structure and equations
.
Environmental Geology
53
(
6
),
1299
1309
.
Chen
R. S.
,
Kang
E. S.
,
Lu
S. H.
,
Ji
X. B.
,
Zhang
Z. H.
,
Yang
Y.
&
Qing
W. W.
2008b
A distributed water–heat coupled model for mountainous watershed of an inland river basin in Northwest China (II) using meteorological and hydrological data
.
Environmental Geology
55
(
1
),
17
28
.
Cheng
G. D.
&
Wu
T. H.
2007
Responses of permafrost to climate change and their environmental significance, Qinghai-Tibet Plateau
.
Journal of Geophysical Research Atmospheres
112
(
112
),
93
104
.
Cherkauer
K. A.
&
Lettenmaier
D. P.
1999
Hydrologic effects of frozen soils in the upper Mississippi River basin
.
Journal of Geophysical Research Atmospheres
104
(
D16
),
19599
19610
.
Cherkauer
K. A.
,
Bowling
L. C.
&
Lettenmaier
D. P.
2003
Variable infiltration capacity cold land process model updates
.
Global & Planetary Change
38
(
1–2
),
151
159
.
Foley
J. A.
,
Prentice
I. C.
,
Ramankutty
N.
,
Levis
S.
,
Pollard
D.
,
Sitch
S.
&
Haxeltine
A.
1996
An integrated biosphere model of land surface processes, terrestrial carbon balance, and vegetation dynamics
.
Global Biogeochemical Cycles
10
(
4
),
603
628
.
Han
L. J.
,
Tsunekawa
A.
,
Tsubo
M.
,
He
C. Y.
&
Shen
M. G.
2014
Spatial variations in snow cover and seasonally frozen ground over northern China and Mongolia, 1988–2010
.
Global & Planetary Change
116
(
3
),
139
148
.
Harlan
R. L.
1973
Analysis of coupled heat-fluid transport in partially frozen soil
.
Water Resources Research
9
(
5
),
1314
1323
.
Hayashi
M.
,
Goeller
N.
,
Quinton
W. L.
&
Wright
N.
2007
A simple heat-conduction method for simulating the frost-table depth in hydrological models
.
Hydrological Processes
21
(
19
),
2610
2622
.
Homan
J. W.
,
Kane
D. L.
&
Sturm
M.
2011
Arctic snow distribution patterns at the watershed scale
.
Hyrology Research
46
(
4
),
507
520
.
Jia
Y. W.
&
Tamai
N.
1998
Integrated analysis of water and heat balances in Tokyo Metropolis with a distributed model
.
Journal of Japan Society of Hydrology & Water Resources
11
(
11
),
150
163
.
Jia
Y. W.
,
Ni
G. H.
,
Kawahara
Y.
&
Suetsugi
T.
2001
Development of WEP model and its application to an urban watershed
.
Hydrological Processes
15
(
11
),
2175
2194
.
Jia
Y.
,
Ding
X.
,
Qin
C.
&
Wang
H.
2009
Distributed modeling of landsurface water and energy budgets in the inland Heihe river basin of China
.
Hydrology & Earth System Science
13
(
10
),
1849
1866
.
Jiang
F.
,
Liu
S.
,
Liu
J.
&
Wang
X. Y.
2009
Measurement of ice movement in water using electrical capacitance tomography
.
Journal of Thermal Science
18
(
1
),
8
12
.
Konrad
J. M.
&
Duquennoi
C.
1993
A model for water transport and ice lensing in freezing soils
.
Water Resources Research
29
(
3
),
3109
3124
.
Li
Q.
,
Sun
S. F.
&
Xue
Y. K.
2010
Analyses and development of a hierarchy of frozen soil models for cold region study
.
Journal of Geophysical Research Atmospheres
115
(
D3
),
315
317
.
Liang
X.
,
Lettenmaier
D. P.
,
Wood
E. F.
&
Burges
S. J.
1994
A simple hydrologically based model of land surface water and energy fluxes for general circulation models
.
Journal of Geophysical Research Atmospheres
99
(
D7
),
14415
14428
.
Liu
Z.
&
Yao
Z.
2016
Contribution of glacial melt to river runoff as determined by stable isotopes at the source region of the Yangtze River, China
.
Hydrology Research
47
(
2
),
442
453
.
Liu
J. S.
,
Hayakawa
N.
,
Lu
M. J.
,
Dong
S. H.
&
Yuan
J. Y.
2003
Hydrological and geocryological response of winter streamflow to climate warming in Northeast China
.
Cold Region Science & Technology
37
(
1
),
15
24
.
Monteith
J. L.
1973
Principles of Environmental Physics
.
Edward Arnold Publishers
,
London
.
Niu
G. Y.
&
Yang
Z. L.
2006
Assessing a land surface model's improvements with GRACE estimates
.
Geophysical Research Letters
33
(
7
),
L07401
.
Niu
L.
,
Ye
B. S.
,
Ding
Y. J.
,
Li
J.
,
Zhang
Y. S.
,
Sheng
Y.
&
Yue
G. Y.
2016
Response of hydrological processes to permafrost degradation from 1980 to 2009 in the Upper Yellow River Basin, China
.
Hydrology Research
47
(
5
),
1014
1024
.
Noilhan
J.
&
Planton
S.
1989
A simple parameterization of land surface processes for meteorological models
.
Monthly Weather Review
117
(
3
),
536
549
.
Pomeroy
J. W.
,
Gray
D. M.
,
Brown
T.
,
Hedstrom
N. R.
,
Quinton
W. L.
,
Granger
R. J.
&
Carey
S. K.
2007
The cold regions hydrological model: a platform for basing process representation and model structure on physical evidence
.
Hydrological Processes
21
(
19
),
2650
2667
.
Refsgaard
J. C.
,
Storm
B.
,
Abbott
M. B.
1996
Comment on ‘A discussion of distributed hydrological modeling’
. In:
Distributed Hydrological Modeling
(
Abbott
M. B.
&
Refsgaard
J. C.
, eds).
Kluwer Academic Publishers
,
Dordrecht
, pp.
279
287
.
Shang
S. H.
,
Lei
Z. D.
&
Yang
S. X.
1997
Numerical simulation improvement of coupled moisture and heat transfer during soil freezing
.
Journal of Tsinghua University (Sci& Tech)
8
,
62
64
.
Shrestha
M.
,
Wang
L.
,
Koike
T.
,
Xue
Y.
&
Hirabayashi
Y.
2010
Improving the snow physics of WEB-DHM and its point evaluation at the SnowMIPsites
.
Hydrology & Earth System Sciences
14
(
12
),
2577
2594
.
Slater
A. G.
,
Pitman
A. J.
&
Desborough
C. E.
1998
Simulation of freeze–thaw cycles in a general circulation model land surface scheme
.
Journal of Geophysical Research Atmospheres
103
(
D10
),
11303
11312
.
Smith
N. V.
,
Saatchi
S. S.
&
Randerson
J. T.
2004
Trends in high northern latitude soil freeze and thaw cycles from 1988 to 2002
.
Journal of Geophysical Research Atmospheres
109
(
12
),
221
236
.
Tian
K. M.
,
Liu
J. S.
,
Kang
S. C.
,
Campbell
I. B.
,
Zhang
F.
,
Zhang
Q. G.
&
Lu
W.
2009
Hydrothermal pattern of frozen soil in Nam Co lake basin, the Tibetan Plateau
.
Environmental Geology
57
(
8
),
1775
1784
.
Verdin
K. L.
&
Verdin
J. P.
1999
A topological system for delineation and codification of the Earth's river basins
.
Journal of Hydrology
218
(
1–2
),
1
12
.
Wang
A. W.
,
Xie
Z. H.
,
Feng
X. B.
,
Tian
X. J.
&
Qin
P. H.
2014
A soil water and heat transfer model including changes in soil frost and thaw fronts
.
Science China Earth Sciences
57
(
6
),
1325
1339
.
Wang
X. J.
,
Yang
M. X.
,
Pang
G. J.
,
Wan
G. N.
&
Chen
X. L.
2015
Simulation and improvement of land surface processes in Nameqie, Central Tibetan Plateau, using the Community Land Model (CLM3.5)
.
Environmental Earth Sciences
73
(
11
),
7343
7357
.
Woo
M. K.
,
Arain
M. A.
,
Mollinga
M.
&
Yi
S.
2004
A two-directional freeze and thaw algorithm for hydrologic and land surface modelling
.
Geophysical Research Letters
31
(
12
),
261
268
.
Wu
M.
,
Kang
W.
,
Xiao
T.
&
Huang
J.
2013
Water movement in soil freezing and thawing cycles and flux simulation
.
Advances in Water Science
24
(
4
),
543
550
.
Xiang
X. H.
,
Wu
X. L.
,
Wang
C. H.
,
Xi
C.
&
Shao
Q. Q.
2013
Influences of climate variation on thawing–freezing processes in the northeast of Three-River Source Region China
.
Cold Regions Science & Technology
86
(
2
),
86
97
.
Yamazaki
Y.
,
Kubota
J.
,
Ohata
T.
,
Vuglinsky
V.
&
Mizuyama
T.
2006
Seasonal changes in runoff characteristics on a permafrost watershed in the southern mountainous region of eastern Siberia
.
Hydrological Processes
20
(
3
),
453
467
.
Zaradny
H.
1993
Groundwater Flow in Saturated and Unsaturated Soil
.
Balkema Press
,
Rotterdam
.
Zhang
Y. L.
,
Cheng
G. D.
,
Li
X.
,
Han
X. J.
,
Wang
L.
,
Li
H. Y.
,
Chang
X. L.
&
Flerchinger
G. N.
2013
Coupling of a simultaneous heat and water model with a distributed hydrological model and evaluation of the combined model in a cold region watershed
.
Hydrological Processes
27
(
25
),
3762
3776
.
Zhang
D.
,
Cong
Z.
,
Ni
G.
,
Yang
D.
&
Hu
S.
2015
Effects of snow ratio on annual runoff within the Budyko framework
.
Hydrology & Earth System Sciences
19
,
1977
1992
.
Zhao
Q.
,
Ye
B.
,
Ding
Y.
,
Zhang
S.
,
Yi
S.
,
Wang
J.
,
Shangguan
D.
,
Zhao
C.
&
Han
H.
2013
Coupling a glacier melt model to the variable infiltration capacity (VIC) model for hydrological modeling in north-western China
.
Environmental Earth Sciences
68
(
1
),
87
101
.
Zhou
J.
,
Pomeroy
J. W.
,
Zhang
W.
,
Cheng
G.
,
Wang
G.
&
Chen
C.
2014
Simulating cold regions hydrological processes using a modular model in the west of China
.
Journal of Hydrology
509
(
4
),
13
24
.