Development of WEP-COR model to simulate land surface water and energy budgets in a cold region

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.

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.

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

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.

The soil freezing status is divided into three types. (1) For the unfrozen type, the soil temperature (T_s) is above 0 °C, and there is no solid water in the soil layer. (2) For the frozen type, T_s is lower than T_f (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 cm³/cm³ (Wu et al. 2013). Here, the T_f value was revised to −10 °C. (3) For the partially frozen type, T_s is higher than T_f but lower than 0 °C. Liquid water and solid water can coexist in the soil.

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. T_t+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.

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 × 10⁴ km² (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 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.

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 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 cm³/cm³, respectively.

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.

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

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

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 cm³/cm³, 0.034 cm³/cm³, 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.

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