This study simulated the processes of soil freezing and thawing and coupled water-heat movement based on Harlan’s model. Comparing water movement in frozen soil to that in unsaturated soil, Harlan firstly developed a coupled model to analyze water and heat movement in frozen soil under some hypotheses (Harlan, 1971; Harlan, 1973). Subsequently, a numerical model that adopted the finite difference method was derived by Kung and Steenhuis (1986) to simulate soil water and heat movement. This model overcame limitations in description of water flux for a partly frozen soil in Harlan’s model. Based on Harlan’s theories, Taylor, Sheppard, Fukuda, et al. (Sheppard et al., 1978; Taylor and Luthin, 1978; Fukuda, 1980; Fukuda et al., 1987) analyzed and improved Harlan’s model or developed their own models.

Soil water movement for the processes of soil freezing and thawing was described by Richards equation as follows: ∂ θ u ∂ t = ∂ ∂ z [ D ( θ u ) ∂ θ u ∂ z − K ( θ u ) ] − ρ i ρ w ∂ θ i ∂ t , (1) where θ u is the volumetric liquid content ( m 3 ∙ m − 3 ); θ i is the ice content ( m 3 ∙ m − 3 ); ρ i is the density of ice ( 900 k g ∙ m − 3 ); ρ w is the density of water ( 1000 k g ∙ m − 3 ); K ( θ u ) is the unsaturated hydraulic conductivity of soil ( m ∙ s − 1 ); D ( θ u ) is the diffusivity of frozen soil ( m 2 ∙ s − 1 ); t is time (s) and z is depth (m).

In this paper, due to lack of available soil data in TRSR, the unsaturated hydraulic conductivity and the diffusivity of frozen soil are given as the following empirical formulas respectively (Lei et al., 1988; Hao et al., 2009): D ( θ u ) = D 1 ( θ u / θ s ) D 2 , (2) K ( θ u ) = K 1 ( θ u / θ s ) K 2 , (3) where D1, D2, K1 and K2 are parameters, D1 = 226.4, D2 = 8.4, K1 = 0.9, K2 = 10.87 (Lei et al., 1988; Hao et al., 2009); θs is saturated water content ( m 3 ∙ m − 3 ).

For partially frozen soil, hydraulic conductivity and diffusivity are adjusted by the following formulas respectively (Taylor and Luthin, 1978; Lei et al., 1988): D i ( θ u ) = D ( θ u ) / 10 10 θ i   , (4) K i ( θ u ) = K ( θ u ) / 10 10 θ i ,   (5) where D i ( θ u ) is the adjusted diffusivity for partially frozen soil ( m 2 ∙ s − 1 ); K i ( θ u ) is the adjusted hydraulic conductivity in the partially frozen soil ( m ∙ s − 1 ).

The one-dimensional heat flux equation of frozen soil is given by: C V S ∂ T ∂ t = ∂ ∂ z [ λ ∂ T ∂ z ] + L ρ i ∂ θ i ∂ t , (6) where C v s is volumetric heat capacity of soil ( J ∙ m − 3 ∙ K − 1 ); λ is thermal conductivity of soil ( W ∙ m − 1 ∙ K − 1 ); L is latent heat of fusion ( J ∙ k g − 1 ); T is soil temperature (°C); other symbols are the same as before.

Volumetric heat capacity, C v s , of soil is from De Vries (De Vries and Van Wijk, 1963): C v s = c s ( 1 − θ s a t ) + θ l i q c l i q + θ i c e c i c e , (7) where θ l i q ( m 3 ∙ m − 3 ) and θ i c e ( m 3 ∙ m − 3 ) are liquid water content and ice content, respectively; θ s a t is the saturated volumetric water content ( m 3 ∙ m − 3 ) ; C s is the volumetric heat capacity of the soil solids ( J ∙ m − 3 ∙ K − 1 ); c l i q and c i c e are the specific heat capacities of liquid water and ice, respectively ( 4.213 J ∙ m − 3 ∙ K − 1 ; 1.94 J ∙ m − 3 ∙ K − 1 ). So Equation 7 can be written as follows: C v s = c + 4.213 × θ l i q + 1.94 × θ i c e , (8) where c is a parameter about the soil properties and varies with the sand and clay content, which is mainly influenced by the regional climate including long term average temperature and precipitation and land cover. It also exhibits an obvious vertical distribution which could be characterized by latitude-longitude and elevation. As an important parameter in frozen soil water-heat movement model, this paper established a multiple regression model between c and factors that significantly correlated with it (i.e. climate characteristics, vegetation characteristics and geographical characteristics). After that, the multiple regression model was applied to calculate the spatial distribution of c in TRSR for simulating the processes of soil freezing and thawing. The calculation were conducted at a resolution of 1 km by 1 km to match the vegetation coverage and MODIS LST data.

Soil thermal conductivity is highly complex because it is not only related to the proportion of the components, but also related to the structure and shape of soil components and other factors. Considering this, in this paper, an empirical formula was used for estimation of soil thermal conductivity (Shi et al., 1998). λ = λ m − ( λ m − λ T ) θ u θ u + ρ i ρ w ∙ θ i   , (9) where λ m is thermal conductivity of frozen soil ( W ∙ m − 1 ∙ K − 1 ); λ T is thermal conductivity of unfrozen soil ( W ∙ m − 1 ∙ K − 1 ). Based on the measured data, Pang et al. (2006) calculated the values of λ m and λ T in this region and λ m = 1.84 W ∙ m − 1 ∙ K − 1 and λ T = 1.57 W ∙ m − 1 ∙ K − 1 .

In the processes of soil freezing and thawing, the connection between water and heat movement of frozen soil is mainly reflected in the dynamic equilibrium of water content and soil subzero temperature, namely a one-one correspondence exists between the soil subzero temperature and water content. Therefore, the connecting equation is written as (Shang et al., 1997): θ u = θ max ( T ) , (10) where θ max ( T ) is the largest unfrozen water content in the condition of the corresponding soil subzero temperature. Thus, the soil water content coupled with heat transfer model comprises Eqss. 1, 6, 10.

However, due to the different initial water content and the hysteresis effect of soil water retention characteristics, the relationship between water content and subzero temperature is not a one-to-one correspondence single value function. In this paper, the relationship is given as the following empirical formula (An et al., 1987): { θ u = 0.0644 ⁡ exp ( 0.0551 T )   T < − 0.75 θ u = 0.3058 + 0.596 ( T + 0.30 ) − 0.75 ≤ T ≤ − 0.30 θ u = 0.3058   T ≥ − 0.30 , (11)

Due to the nonlinearity of soil water-heat movement basic equations, the anisotropy of soils and the complexity of the initial and boundary conditions, a numerical simulation method instead of analytical methods was generally used to solve these equations. In this paper, the finite difference method was used to numerically simulate the water–heat coupled movements.

Using a central difference scheme, Equations 1, 6 can be transferred respectively as follows: α i ( θ u ) i − 1 k + 1 + β i ( θ u ) i k + 1 + γ i ( θ u ) i + 1 k + 1 = h i   , (12) A i T i − 1 k + 1 + B i T i k + 1 + C i T i + 1 k + 1 = D i , (13) where subscript i is the ith soil layer, superscript k is the number of computation time intervals in a time step. The coefficients of α _i , β _i , γ _i , h _i , A _i , B _i , C _i and D _i can be calculated from initial conditions of soil temperature and water content.

Equations 12 and 13 are all in a set of tridiagonal equations, which can be solved by a chasing method. The time step was set as 1 day and the spatial step was 0.01 m. So, on the basis of certain initial and boundary conditions, the numerical solution of water and heat coupled movement in frozen soil can be obtained.

The data were collected from 19 meteorological stations in and around TRSR. The station information is presented in S1, and the geographical location of the stations is shown in Figure 1. The data includes daily air temperature, precipitation, frozen depth, thawed depth, surface temperature and soil temperatures at 5, 10, 15, 20 and 40 cm depth from 1 January 2006 to 31 December 2008, and the missing data were obtained by linear interpolation.

The index of vegetation coverage was used to analyze the geographical law of parameters c in TRSR. The parameter c in the coupled model relates to the soil properties, and experiments showed that Normalized Difference Vegetation Index (NDVI) is sensitive to the change of soil properties (Baret and Guyot, 1991; Fang and Tian, 1998). NDVI is a comprehensive reflection of the vegetation types, form and growth status in the unit pixel and its value depends on vegetation coverage and leaf area index (LAI) and other factors, thus providing a measure of vegetation coverage (Carlson and Ripley, 1997). The NDVI data were acquired from NASA (https://ladsweb.nascom.nasa.gov).

A dimidiate pixel model, which is a simplified linear spectral unmixing method, was used to calculate fractional vegetation coverage (FVC) due to its simplicity in principle and form (Ivits et al., 2013). Many studies examined this model by theoretical method and applied research (Leprieur et al., 1994; QI et al., 2000; Zribi et al., 2003). An assumption of the model is that a pixel contains only two elements of vegetation and soil, which can be represented as follows: F V C = ( N D V I − N D V I s o i l ) / ( N D V I v e g − N D V I s o i l ) , (14) where NDVI used in this study was organized at 1 km by 1 km spatial resolution with a 16-days time interval from Moderate-resolution Imaging Spectroradiometer (MODIS). N D V I s o i l and N D V I s o i l are two key input parameters that represent the NDVI values of pure pixels of bare soil and vegetation, respectively. For more detailed processes of FVC estimation, readers can refer to Zeng et al. (2000) for further information.

Figure 2 shows the spatial distribution of FVC that used in this study, i.e. the mean value of 12 months from July 2006 to June 2007 in the TRSR, with the vegetation coverage gradually deteriorating from southeast to northwest. The cover in the east of YeSR represent the highest (60%–90%) in TRSR. The vegetation coverage in the LcSR, southeast YaSR and southwest YeSR are relatively high (40%–60%). However, the vegetation coverage was relatively low and generally below 40% in the northwestern YeSR and west of YaSR.

Land surface temperature (LST) from MODIS were used as the input data in the frozen soil water-heat model in this study. It plays an important role in material exchange and energy balance between the surface and the atmosphere, and it is a key variable in the study of land surface physical processes at regional and global scale (Wang et al., 2012). Due to the large area of missing pixels of daily MODIS LST product in TRSR, we used the MYD11A2 product with 1 km spatial resolution. This product provides 8-days LST data recovered by the spilt-window algorithm. The empty pixels of MYD11A2 product were filtered from its neighbors by using the spline interpolation. Then we used the liner interpolation method to obtain the daily LST data in TRSR.

The MODIS/AQUA LST products provide two instantaneous observations every day, and the overpass times are in 1:30a.m. and 1:30p.m. In this paper, the daily land surface temperature was estimated by the average of daytime LST and nighttime LST that are respectively closed to the minimum and maximum temperature in a day (Geerts, 2003).

The model performance is evaluated by correlation coefficient r, Nash-Sutcliffe efficiency coefficient NSE and the root mean square error RMSE. The methods of calculating r, NSE and RMSE are shown as follows: r = ∑ i = 1 N ( x i − x ¯ ) ( y i − y ¯ ) ∑ i = 1 N ( x i − x ¯ ) 2 ∑ i = 1 N ( y i − y ¯ ) 2   , (15) N S E = 1 − ∑ i = 1 N ( x i − y i ) 2 ∑ i = 1 N ( y i − y ¯ ) 2 , (16) R M S E = ∑ i = 1 N ( x i − y i ) 2 N , (17) Where N is the number of sample; x i and y i is the simulation value and observation value respectively; x ¯ and y ¯ are mean of x i and y i respectively.

In TRSR, Jiuzhi and Maduo are the typical stations in the processes of soil freezing and thawing. Figure 3 shows that the processes can be divided into four stages.

1) Unstable freezing phase. In this phase, the soil temperatures periodically decrease to less than 0°C with the air temperatures fluctuate around 0°C, which makes the shallow-layer soil periodically experience freezing and thawing. When the negative accumulated temperature is higher than 0°C and the shallow-layer soil is always frozen in the daytime.

Model calibration and verification were carried out at 19 stations in and around TRSR by comparing simulated and observed freezing depths and soil temperatures. The observation data were separated into two periods, 1 July 2006-30 June 2007 and 1 July 2007-30 June 2008, for model calibration and verification, respectively. The calibration results are shown in Table 1.

Taking Xinghai and Dari stations for examples, NSE values for the soil freezing depths were 0.92 and 0.75 during the calibration period, respectively. The NSE for the daily soil temperatures at the Xinghai station for 5, 10, 15, 20, 40 cm were 0.98, 0.98, 0.98, 0.97 and 0.94, respectively. The NSE for the daily soil temperatures at the Dari station for 5, 10, 15, 20, 40 cm were 0.98, 0.99, 0.99, 0.99 and 0.95, respectively.

During the verification period, the r, NSE and RMSE values for simulation of freezing depths of 19 meteorological stations are presented in Table 2. Table 2 shows that the lowest Nash-Sutcliffe efficiency coefficient NSE and correlation coefficient r in simulated freezing depths at Maduo station (Figure 4B) are 0.4 and 0.73, respectively. The highest NSE and r at Zaduo station (Figure 4A) are 0.94 and 0.97, respectively. The average values of NSE and r are 0.89 and 0.77, respectively. It is proved that the model can capture the characteristics of the soil freezing-thawing processes well.

The simulated and observed frozen depths at Dari and Xinghai stations during the period from 1 July 2006 to 30 June 2008 are shown in Figure 4. The simulation of frozen depths is more accurate than the thawed depths. The simulated and observed daily soil temperatures at Dari station for 5, 10, 15, 20 and 40 cm are shown in Figures 5A–E. The r and NSE values between observed and simulated soil temperatures at Dari station for 5, 10, 15, 20 and 40 cm are all above 0.98 and 0.95, respectively (Figures 5F–J). The RMSE for the five layers is less than 1.52°C (Figures 5F–J). Therefore, the simulated values agreed well with the observed values.